Entropy Multiparticle Correlation Expansion for a Crystal

Abstract

As first shown by H. S. Green in 1952, the entropy of a classical fluid of identical particles can be written as a sum of many-particle contributions, each of them being a distinctive functional of all spatial distribution functions up to a given order. By revisiting the combinatorial derivation of the entropy formula, we argue that a similar correlation expansion holds for the entropy of a crystalline system. We discuss how one- and two-body entropies scale with the size of the crystal, and provide fresh numerical data to check the expectation, grounded in theoretical arguments, that both entropies are extensive quantities.

1. Introduction

The entropy multiparticle correlation expansion (MPCE) is an elegant statistical-mechanical formula that entails the possibility of reconstructing the total entropy of a many-particle system term by term, including at each step of summation the integrated contribution from spatial correlations between a specified number of particles.

The original derivation of the entropy MPCE is found in a book by H. S. Green (1952) [1]. Green’s expansion applies for the canonical ensemble (CE). In 1958, Nettleton and M. S. Green derived an apparently different expansion valid in the grand-canonical ensemble (GCE) [2]. It took the ingenuity of Baranyai and Evans to realize, in 1989, that the CE expansion can indeed be reshuffled in such a way as to become formally equivalent to the GCE expansion [3].

A decisive step forward was eventually taken by Schlijper [4] and An [5], who have highlighted the similarity of the entropy formula to a cumulant expansion, and the close relationship with the cluster variation method (see, e.g., [6]). Other papers wherein in various ways the combinatorial content of the entropy MPCE is emphasized are references [7,8,9,10].

Since the very beginning it has been clear that the successive terms in the entropy expansion for a homogeneous fluid are not all of equal importance. In particular, the contributions from correlations between more than two particles are only sizable at moderate and higher densities. However, while the two-body entropy is easily accessed in a simulation, computing the higher-order entropy terms is a prohibitive task (see, however, reference [11]). Hence, the only viable method to compute the total entropy in a simulation remains thermodynamic integration (see e.g., [12]). The practical interest for the entropy expansion has thus shifted towards the residual multiparticle entropy (RMPE), defined as the difference between excess entropy and two-body entropy. The RMPE is a measure of the impact of non-pair multiparticle correlations on the entropy of the fluid. For hard spheres, Giaquinta and Giunta have observed that the RMPE changes sign from negative to positive very close to freezing [13]. At low densities the RMPE is negative, reflecting a global reduction (largely driven by two-body correlations) of the phase space available to the system as compared to the ideal gas. The change of sign of the RMPE close to freezing indicates that fluid particles, which at high enough densities are forced by more stringent packing constraints, start exploring, this time in a cooperative way, a different structural condition on a local scale, preluding to crystallization on a global scale. Since the original observation in [13], a clear correspondence between the RMPE zero and the ultimate threshold for spatial homogeneity in the system has been found in many simple and complex fluids [14,15,16,17,18,19,20,21,22,23,24], thereby leading to the belief that the vanishing of the RMPE is a signature of an impending structural or thermodynamic transition of the system from a less ordered to a more spatially organized condition (freezing is just an example of many). Albeit empirical, this entropic criterion is a valid alternative to the far more demanding exact free-energy methods when a rough estimate of the transition point is deemed sufficient. For a simple discussion of the interplay between entropy and ordering, the reader is referred to reference [25]; see instead references [26,27] for general considerations about the entropy of disordered solids.

A pertinent question to ask is, what happens to the RMPE on the solid side of the phase boundary, considering that an entropy expansion also holds for the crystal? This is precisely the problem addressed in this paper. Can the scope of the entropic criterion be extended in such a way that it also applies for melting? As it turns out, we can offer no definite answer to this question, since theory alone does not go far enough and we ran into a serious computational bottleneck: while the formulae are clear and the numerical procedure is straightforward, it is extremely hard to obtain reliable data for the two-body entropy of a three-dimensional crystal. We have only carried out a limited test on a triangular crystal of hard disks, but our results are affected by finite-size artifacts that make them inconclusive. Nevertheless, a few firm points have been established: (1) the approximate entropy expressions obtained by truncating the MPCE at a given order can all be derived from an explicit functional of the correlation functions up to that order; (2) the one-body entropy for a crystal is an extensive quantity (the same is held to be true for the two-body entropy, but our arguments are not sufficient for a proof); (3) the peaks present in the crystal one-body density have a nearly Gaussian shape; (4) we have also clarified the role of lattice symmetries in dictating the structure of the two-body density, which is explicitly determined at zero temperature.

This paper is organized as follows. In Section 2 we resume the formalism of the entropy expansion for homogeneous fluids and provide the basic tools needed for its extension to crystals. Then, in Section 3 we exploit the symmetries of one- and two-body density functions to predict the scaling of one- and two-body entropies with the size of the crystal. The final Section 4 is reserved to concluding remarks.

2. Derivation of the Entropy MPCE

In this Section, we collect a number of well-established results on the entropy MPCE, with the only purpose of setting the language and notation for the rest of the paper. First, we recall the derivation of the entropy formula for a one-component system of classical particles in the canonical ensemble. Such an ensemble choice is by no means restrictive, since, as we show next, it is always possible to take advantage of the sum rules obeyed by the canonical correlation functions to arrange the entropy MPCE in an ensemble-invariant form. Then, in the following Section we present an application of the formalism to crystals.

The canonical partition function of a system of N classical particles of mass m at temperature T is $Z_N=Z_N^{\mathrm{id}}{Z}_N^{\mathrm{exc}}$, where the ideal and excess parts are given by

$$Z_N^{\mathrm{id}}=\frac{1}{N!}{\left(\frac{V}{{\Lambda }^3}\right)}^N\mathrm{and}{Z}_N^{\mathrm{exc}}=\frac{1}{V^N}\int {\mathrm{d}}^3{R}_1\cdots {\mathrm{d}}^3{R}_N{e}^{-\beta U\left({\mathbf{R}}^N\right)}.$$

(1)

In Equation (1), V is the system volume, $\beta =1/\left(k_{B}T\right)$, $\Lambda =h/\sqrt{2\pi m{k}_{B}T}$ is the thermal wavelength, and $U\left({\mathbf{R}}^N\right)$ is an arbitrary potential energy. As the particles are identical, for each $n=1,2,\dots ,N$ the cumulative sum of all n-body terms in U is invariant under permutations of particle coordinates (we can also say that U is $S_N$-invariant, $S_N$ being the symmetric group of the permutations on N symbols). The CE average of a function f of coordinates reads

$$〈f\left({\mathbf{R}}^N\right)〉\equiv \frac{1}{V^N}\int {\mathrm{d}}^3{R}_1\cdots {\mathrm{d}}^3{R}_{N}f\left({\mathbf{R}}^N\right){\pi }_{\mathrm{can}}\left({\mathbf{R}}^N\right)\mathrm{with}{\pi }_{\mathrm{can}}\left({\mathbf{R}}^N\right)=\frac{e^{-\beta U\left({\mathbf{R}}^N\right)}}{Z_N^{\mathrm{exc}}},$$

(2)

where ${\pi }_{\mathrm{can}}\left({\mathbf{R}}^N\right)$ is the configurational part of the canonical density function. Finally, the excess entropy $S_N^{\mathrm{exc}}\equiv S_N-S_N^{\mathrm{id}}$ reads

$$\frac{S_N^{\mathrm{exc}}}{k_B}=-\frac{1}{V^N}\int {\mathrm{d}}^3{R}_1\cdots {\mathrm{d}}^3{R}_N{\pi }_{\mathrm{can}}\left({\mathbf{R}}^N\right)ln{\pi }_{\mathrm{can}}\left({\mathbf{R}}^N\right)=-〈ln{\pi }_{\mathrm{can}}\left({\mathbf{R}}^N\right)〉.$$

(3)

We define a set of marginal density functions (MDFs) by

$$\begin{array}{ccc}\hfill P^{\left(N\right)}\left({\mathbf{R}}^N\right)& =& {\pi }_{\mathrm{can}}\left({\mathbf{R}}^N\right);\hfill \\ \hfill P^{\left(n\right)}\left({\mathbf{R}}^n\right)& =& \frac{1}{V^{N-n}}\int {\mathrm{d}}^3{R}_{n+1}\cdots {\mathrm{d}}^3{R}_N{\pi }_{\mathrm{can}}\left({\mathbf{R}}^N\right)(n=1,\dots ,N-1).\hfill \end{array}$$

(4)

Owing to $S_N$-invariance of ${\pi }_{\mathrm{can}}$, it makes no difference which vector radii are integrated out in Equation (4); hence, $P^{\left(n\right)}\left({\mathbf{r}}^n\right)$ is $S_n$-invariant (for example, $P^{\left(2\right)}(\mathbf{r},{\mathbf{r}}^{\prime })=P^{\left(2\right)}({\mathbf{r}}^{\prime },\mathbf{r})$). The following properties are obvious:

$$\frac{1}{V^n}\int {\mathrm{d}}^3{R}_1\cdots {\mathrm{d}}^3{R}_n{P}^{\left(n\right)}\left({\mathbf{R}}^n\right)=1\mathrm{and}\frac{1}{V}\int {\mathrm{d}}^3{R}_{n+1}{P}^{(n+1)}\left({\mathbf{R}}^{n+1}\right)=P^{\left(n\right)}\left({\mathbf{R}}^n\right).$$

(5)

Then, the n-body density functions (DFs), for $n=1,\dots ,N$, can be expressed as

$${\rho }^{\left(n\right)}\left({\mathbf{r}}^n\right)\equiv 〈\underset{i_1\dots i_n}{{\sum }^{\prime }}{\delta }^3({\mathbf{R}}_{i_1}-{\mathbf{r}}_1)\cdots {\delta }^3({\mathbf{R}}_{i_n}-{\mathbf{r}}_n)〉=\frac{N!}{(N-n)!}\frac{P^{\left(n\right)}\left({\mathbf{r}}^n\right)}{V^n},$$

(6)

where the sum in (6) is carried out over all n-tuples of distinct particles (for example, the sum for $n=2$ contains $N(N-1)$ terms). We note that $P^{\left(1\right)}=1$ and ${\rho }^{\left(1\right)}=N/V\equiv \rho$ if no one-body term is present in U, i.e., if no external potential acts on the particles (then U is translationally invariant). $P^{\left(1\right)}\left(\mathbf{r}\right)/V$ is the probability density of finding a particle in $\mathbf{r}$; hence, ${\rho }^{\left(1\right)}\left(\mathbf{r}\right)=N{P}^{\left(1\right)}\left(\mathbf{r}\right)/V$ is the number density at $\mathbf{r}$. Similarly, $P^{\left(2\right)}(\mathbf{r},{\mathbf{r}}^{\prime })/V^2$ is the probability density of finding one particle in $\mathbf{r}$ and another particle in ${\mathbf{r}}^{\prime }$; hence, ${\rho }^{\left(2\right)}(\mathbf{r},{\mathbf{r}}^{\prime })=N(N-1)P^{\left(2\right)}(\mathbf{r},{\mathbf{r}}^{\prime })/V^2$ is the density of the number of particle pairs at $(\mathbf{r},{\mathbf{r}}^{\prime })$. As ${\mathbf{r}}^{\prime }$ increasingly departs from $\mathbf{r}$, the positions of two particles become less and less correlated, until $P^{\left(2\right)}(\mathbf{r},{\mathbf{r}}^{\prime })=P^{\left(1\right)}\left(\mathbf{r}\right)P^{\left(1\right)}\left({\mathbf{r}}^{\prime }\right)$ at infinite distance. We stress that this cluster property holds in full generality, even for a broken-symmetry phase.

The n-body reduced density functions, for $n=2,\dots ,N$, read

$$\begin{array}{ccc}& & g^{\left(n\right)}\left({\mathbf{r}}^n\right)\equiv \frac{{\rho }^{\left(n\right)}\left({\mathbf{r}}^n\right)}{{\rho }^{\left(1\right)}\left({\mathbf{r}}_1\right)\cdots {\rho }^{\left(1\right)}\left({\mathbf{r}}_n\right)}=\left(1-\frac{1}{N}\right)\cdots \left(1-\frac{n-1}{N}\right)Q^{\left(n\right)}\left({\mathbf{r}}^n\right)\hfill \\ & & \mathrm{with}{Q}^{\left(n\right)}\left({\mathbf{r}}^n\right)=\frac{P^{\left(n\right)}\left({\mathbf{r}}^n\right)}{P^{\left(1\right)}\left({\mathbf{r}}_1\right)\cdots P^{\left(1\right)}\left({\mathbf{r}}_n\right)}.\hfill \end{array}$$

(7)

These functions fulfill the property

$$\frac{1}{V}\int {\mathrm{d}}^3{R}_{n+1}{P}^{\left(1\right)}\left({\mathbf{R}}_{n+1}\right)g^{(n+1)}\left({\mathbf{R}}^{n+1}\right)=\left(1-\frac{n}{N}\right)g^{\left(n\right)}\left({\mathbf{R}}^n\right),$$

(8)

which also holds for $n=1$ if we define $g^{\left(1\right)}\equiv 1$. For a homogeneous fluid, $g^{\left(2\right)}(\mathbf{r},{\mathbf{r}}^{\prime })=g\left(\right|\mathbf{r}-{\mathbf{r}}^{\prime }\left|\right)$. From now on, we adopt the shorthand notation $P_{12\dots n}=P^{\left(n\right)}\left({\mathbf{R}}^n\right)$ and $Q_{12\dots n}=Q^{\left(n\right)}\left({\mathbf{R}}^n\right)$. Moreover, any integral of the kind $V^{-n}\int {\mathrm{d}}^3{R}_1\cdots {\mathrm{d}}^3{R}_n(\cdots )$ is hereafter denoted as $\int (\cdots )$. For example, Equations (3) and (4) indicate that $S_N^{\mathrm{exc}}/k_B=-\int P_{12\dots N}ln{P}_{12\dots N}$.

To build up the CE expansion term by term, our strategy is to consider a progressively larger number of particles. For a one-particle system, the excess entropy in units of the Boltzmann constant is $S_1^{\mathrm{exc}}/k_B=-\int P_{1}ln{P}_1$, leading to a first-order approximation to the excess entropy of a N-particle system in the form $S_N^{\mathrm{exc}}/k_B\approx S_N^{\left(1\right)}/k_B\equiv -N\int P_{1}ln{P}_1$ (that is, each particle contributes to the entropy independently of the other particles). For a two-particle system, the excess entropy is $S_2^{\left(1\right)}$ plus a remainder $k_B{R}_2$, given by:

$$R_2\equiv \frac{S_2^{\mathrm{exc}}-S_2^{\left(1\right)}}{k_B}=-\int P_{12}ln{P}_{12}+2\int P_{1}ln{P}_1=-\int P_{12}ln{Q}_{12}.$$

(9)

Equation (9) suggests a second-order approximation for $S_N^{\mathrm{exc}}$, where each distinct pair of particles contributes the same two-body residual term to the entropy:

$$\frac{S_N^{\left(2\right)}}{k_B}=-N\int P_{1}ln{P}_1-\left(\genfrac{}{}{0pt}{}{N}{2}\right)\int P_{12}ln{Q}_{12}.$$

(10)

Notice that Equation (10) is exact for $N=2$, i.e., $S_2^{\left(2\right)}=S_2^{\mathrm{exc}}$. Similarly, for a three-particle system the excess entropy is $S_3^{\left(2\right)}$ plus a remainder $k_B{R}_3$:

$$\begin{array}{ccc}\hfill R_3\equiv \frac{S_3^{\mathrm{exc}}-S_3^{\left(2\right)}}{k_B}& =& -\int P_{123}ln{P}_{123}+3\int P_{1}ln{P}_1+\left(\genfrac{}{}{0pt}{}{3}{2}\right)\int P_{12}ln{Q}_{12}\hfill \\ & =& -\int P_{123}ln{Q}_{123}+\left(\genfrac{}{}{0pt}{}{3}{2}\right)\int P_{12}ln{Q}_{12}.\hfill \end{array}$$

(11)

Hence, a third-order approximation follows for $S_N^{\mathrm{exc}}$ in the form

$$\frac{S_N^{\left(3\right)}}{k_B}=-N\int P_{1}ln{P}_1-\left(\genfrac{}{}{0pt}{}{N}{2}\right)\int P_{12}ln{Q}_{12}-\left(\genfrac{}{}{0pt}{}{N}{3}\right)\left[\int P_{123}ln{Q}_{123}-\left(\genfrac{}{}{0pt}{}{3}{2}\right)\int P_{12}ln{Q}_{12}\right].$$

(12)

Again, $S_3^{\left(3\right)}=S_3^{\mathrm{exc}}$. Equation (12) reproduces the first three terms in the rhs of Equation (5.9) of reference [8], and one may legitimately expect that the further terms in the entropy expansion are similarly obtained by arguing for $N=4,5,\dots$ like we did for $N=1,2,3$ (see the proof in [8]).

The general entropy formula finally reads:

$$\begin{array}{ccc}\hfill \frac{S_N^{\mathrm{exc}}}{k_B}& =& -\int P_{12\dots N}ln{P}_{12\dots N}=-N\int P_{1}ln{P}_1-\int P_{12\dots N}ln{Q}_{12\dots N}\hfill \\ & =& -N\int P_{1}ln{P}_1-\sum _{n=2}^N\left(\genfrac{}{}{0pt}{}{N}{n}\right)\sum _{a=2}^n{(-1)}^{n-a}\left(\genfrac{}{}{0pt}{}{n}{a}\right)\int P_{1\dots a}ln{Q}_{1\dots a}.\hfill \end{array}$$

(13)

This equation is trivially correct since, for any finite sequence $\left\{c_a\right\}$ of numbers,

$$c_N=\sum _{n=2}^N\left(\genfrac{}{}{0pt}{}{N}{n}\right)\sum _{a=2}^n{(-1)}^{n-a}\left(\genfrac{}{}{0pt}{}{n}{a}\right)c_a.$$

(14)

To prove (14), it is sufficient to observe that, for each fixed $k=2,\dots ,N$, the coefficient of $c_k$ in the above sum is

$$\begin{array}{ccc}\hfill \sum _{n=k}^N{(-1)}^{n-k}\left(\genfrac{}{}{0pt}{}{N}{n}\right)\left(\genfrac{}{}{0pt}{}{n}{k}\right)& =& \sum _{n=0}^{N-k}{(-1)}^n\left(\genfrac{}{}{0pt}{}{N}{n+k}\right)\left(\genfrac{}{}{0pt}{}{n+k}{k}\right)=\left(\genfrac{}{}{0pt}{}{N}{k}\right)\sum _{n=0}^{N-k}{(-1)}^n\left(\genfrac{}{}{0pt}{}{N-k}{n}\right)\hfill \\ & =& \left\{\begin{array}{cc}0,& \mathrm{for}2\le k<N\hfill \\ 1,& \mathrm{for}k=N.\hfill \end{array}\right.\hfill \end{array}$$

(15)

A more compact entropy formula is

$$\frac{S_N^{\mathrm{exc}}}{k_B}=-\sum _{n=1}^N\left(\genfrac{}{}{0pt}{}{N}{n}\right)\sum _{a=1}^n{(-1)}^{n-a}\left(\genfrac{}{}{0pt}{}{n}{a}\right)\int P_{1\dots a}ln{P}_{1\dots a},$$

(16)

which follows from

$$c_N=\sum _{n=1}^N\left(\genfrac{}{}{0pt}{}{N}{n}\right)\sum _{a=1}^n{(-1)}^{n-a}\left(\genfrac{}{}{0pt}{}{n}{a}\right)c_a.$$

(17)

The entropy expansion, (13) or (16), is only valid in the CE. Eliminating $Q_{1\dots a}$ in favor of $g_{1\dots a}$ by Equation (7), an overall constant comes out of the integral in Equation (13), namely,

$$\sum _{n=2}^N\left(\genfrac{}{}{0pt}{}{N}{n}\right)\sum _{a=2}^n{(-1)}^{n-a}\left(\genfrac{}{}{0pt}{}{n}{a}\right)ln\frac{(N-1)\cdots (N-a+1)}{N^{a-1}},$$

(18)

which, by Equation (15), equals $ln\left(N!/N^N\right)$; this term exactly cancels an identical term present in the ideal-gas entropy. In the end, a modified entropy MPCE emerges:

$$\begin{array}{c}\hfill \frac{S_N}{k_B}=N\left[\frac{3}{2}-ln\left(\rho {\Lambda }^3\right)\right]-N\int P_{1}ln{P}_1-\sum _{n=2}^N\left(\genfrac{}{}{0pt}{}{N}{n}\right)\sum _{a=2}^n{(-1)}^{n-a}\left(\genfrac{}{}{0pt}{}{n}{a}\right)\int P_{1\dots a}ln{g}_{1\dots a}.\end{array}$$

(19)

Notice that the first term in the rhs differs by N from the ideal-gas entropy expression in the thermodynamic limit. In order that Equation (19) conforms to the GCE entropy expansion, for each n a suitable fluctuation integral of value $-N/\left[n\right(n-1\left)\right]$ should be summed to (and subtracted from) the n-th term in the expansion. For example, using Equations (7) and (8) the second-order term in (13) can be rewritten as

$$\begin{array}{ccc}\hfill -\left(\genfrac{}{}{0pt}{}{N}{2}\right)\int \frac{{\mathrm{d}}^3{r}_1{\mathrm{d}}^3{r}_2}{V^2}{P}_{12}ln\frac{P_{12}}{P_1{P}_2}& =& \left(\genfrac{}{}{0pt}{}{N}{2}\right)ln\frac{N-1}{N}-\left(\genfrac{}{}{0pt}{}{N}{2}\right)\int \frac{{\mathrm{d}}^3{r}_1{\mathrm{d}}^3{r}_2}{V^2}{P}_1{P}_2{g}_{12}ln{g}_{12}\hfill \\ & =& \left(\genfrac{}{}{0pt}{}{N}{2}\right)ln\frac{N-1}{N}+\frac{N}{2}-\frac{1}{2}{\rho }^2\int {\mathrm{d}}^3{r}_1{\mathrm{d}}^3{r}_2{P}_1{P}_2\left(g_{12}ln{g}_{12}-g_{12}+1\right).\hfill \end{array}$$

(20)

Overall, the extra constants appearing in each term of the entropy formula (for example, the quantity $N/2$ in Equation (20)) add to N. By absorbing such a N in the first term of (19) we recover the ideal-gas entropy in the thermodynamic limit, and the CE expansion becomes formally identical to the grand-canonical MPCE [9].

In Appendix A we present another derivation of the entropy formula in the CE, which is closer in spirit to the one given by H. S. Green. In parallel, we show that the approximation obtained by truncating the MPCE at a given order can be derived from a modified $P_{12\dots N}$ distribution, which is an explicit functional of the spatial correlation functions up to that order.

3. The First Few Terms in the Expansion of Crystal Entropy

The entropy expansion in the CE is formally identical for a fluid system and a crystal, since the origin of (13) is purely combinatorial. However, the DFs of the two phases are radically different: most notably, while $P_1=1$ and ${\rho }^{\left(1\right)}=\rho$ for a homogeneous fluid, the one-body density is spatially structured for a crystal—at least once the degeneracy due to translations and point-group operations has been lifted; we stress that $P_1\ne 1$ only provided that a specific determination of the crystal is taken, since otherwise $P_1=1$ also in the “delocalized” crystalline phase. In practice, in order to fix a crystal in space we should imagine to apply a suitable symmetry-breaking external potential, whose strength is sent to zero after statistical averages have been carried out (in line with Bogoliubov’s advice to interpret statistical averages of broken-symmetry phases as quasiaverages [28], which amounts to sending the strength of the symmetry-breaking potential to zero only after the thermodynamic limit has been taken). A way to accomplish this task is to constrain the position of just one particle. When periodic boundary conditions are applied, keeping one particle fixed will be enough to break the continuous symmetries of free space. As N grows, the effect of the external potential becomes weaker and weaker, since it does not scale with the size of the system.

3.1. One-Body Entropy

A reasonable form of one-body density for a three-dimensional Bravais crystal without defects is the Tarazona ansatz [29]:

$${\rho }^{\left(1\right)}\left(\mathbf{r}\right)={\left(\frac{\alpha }{\pi }\right)}^{3/2}\sum _{\mathbf{R}}{e}^{-\alpha {(\mathbf{r}-\mathbf{R})}^2}=\rho \sum _{\mathbf{G}}{e}^{-\frac{G^2}{4\alpha }}{e}^{i\mathbf{G}·\mathbf{r}},$$

(21)

where $\alpha >0$ is a temperature-dependent parameter, the R’s are direct-lattice vectors, and the G’s are reciprocal-lattice vectors (recall that $\mathbf{G}·\mathbf{R}=2\pi m$ with $m\in \mathbb{Z}$ and ${\int }_V{\mathrm{d}}^{3}rexp\{i(\mathbf{G}+{\mathbf{G}}^{\prime })·\mathbf{r}\}=V{\delta }_{{\mathbf{G}}^{\prime },-\mathbf{G}}$). Equation (21) is a rather generic form of crystal density, which recently we have also applied in a different context [30]. More generally, the one-body density appropriate to a perfect crystal must obey ${\rho }^{\left(1\right)}(\mathbf{r}+\mathbf{R})={\rho }^{\left(1\right)}\left(\mathbf{r}\right)$ for all $\mathbf{R}$, and is thus necessarily of the form

$${\rho }^{\left(1\right)}\left(\mathbf{r}\right)=\sum _{\mathbf{G}}{\tilde{u}}_{\mathbf{G}}{e}^{i\mathbf{G}·\mathbf{r}}\mathrm{with}{\tilde{u}}_{\mathbf{G}}^{*}={\tilde{u}}_{-\mathbf{G}}.$$

(22)

Since ${\int }_V{\mathrm{d}}^{3}r{\rho }^{\left(1\right)}\left(\mathbf{r}\right)=N$, it soon follows ${\tilde{u}}_0=\rho$. Calling $\mathcal{C}$ a primitive cell and $v_0$ its volume, ${\tilde{u}}_{\mathbf{G}}=v_0^{-1}{\int }_{\mathcal{C}}{\mathrm{d}}^{3}r{\rho }^{\left(1\right)}\left(\mathbf{r}\right)exp\{-i\mathbf{G}·\mathbf{r}\}\to 0$ as $G\to \infty$ (by the Riemann–Lebesgue lemma). In real space, a legitimate ${\rho }^{\left(1\right)}\left(\mathbf{r}\right)$ function is ${\sum }_{\mathbf{R}}\varphi (\mathbf{r}-\mathbf{R})$ with $\int {\mathrm{d}}^{3}r\varphi \left(\mathbf{r}\right)=1$ (integration bounds are left unspecified when the integral is over a macroscopic V). In the zero-temperature/infinite-density limit, particles sit at the lattice sites and the one-body density then becomes

$${\rho }^{\left(1\right)}\left(\mathbf{r}\right)=\sum _{\mathbf{R}}{\delta }^3(\mathbf{r}-\mathbf{R}).$$

(23)

Equation (23) is also recovered from Equation (21) in the $\alpha \to \infty$ limit.

For a crystalline solid, the one-body entropy, that is, the first term in the expansion of excess entropy, is (in units of $k_B$):

$$S_1\equiv -\frac{N}{V}\int {\mathrm{d}}^3{r}_1{P}_{1}ln{P}_1=-\int {\mathrm{d}}^3{r}_1{\rho }^{\left(1\right)}\left({\mathbf{r}}_1\right)ln\frac{{\rho }^{\left(1\right)}\left({\mathbf{r}}_1\right)}{\rho }.$$

(24)

One may wonder whether the integral in (24) is $\mathcal{O}\left(N\right)$ in the infinite-size limit. The answer is affirmative, and a simple argument goes as follows. Let ${\rho }^{\left(1\right)}\left(\mathbf{r}\right)$ be ${\sum }_{\mathbf{R}}\varphi (\mathbf{r}-\mathbf{R})$; if $\varphi \left(\mathbf{r}\right)$ is strongly localized near $\mathbf{r}=0$, then ${\rho }^{\left(1\right)}\left(\mathbf{r}\right)\simeq \varphi \left(\mathbf{r}\right)$ in the cell around $\mathbf{R}=0$ and $S_1\simeq -N{\int }_{\mathcal{C}}{\mathrm{d}}^{3}r\varphi \left(\mathbf{r}\right)ln\left\{\varphi \left(\mathbf{r}\right)/\rho \right\}=\mathcal{O}\left(N\right)$. Actually, we can provide a rigorous proof that $S_1$ is negative-semidefinite and its absolute value does not grow faster than N. Using $lnx\le x-1$ for $x>0$ and $xlnx\ge x-1$ for any $x\ge 0$, we obtain

$$0=\rho \int {\mathrm{d}}^3{r}_1\left(\frac{{\rho }_1}{\rho }-1\right)\le \int {\mathrm{d}}^3{r}_1{\rho }_{1}ln\frac{{\rho }_1}{\rho }\le \int {\mathrm{d}}^3{r}_1{\rho }_1\left(\frac{{\rho }_1}{\rho }-1\right)={\rho }^{-1}\int {\mathrm{d}}^3{r}_1{\rho }_1^2-N.$$

(25)

To estimate $\int {\mathrm{d}}^3{r}_1{\rho }_1^2$ we employ the one-body density in (21), which is sufficiently generic for our purposes:

$${\rho }^{-1}\int {\mathrm{d}}^{3}r{\left({\rho }^{\left(1\right)}\left(\mathbf{r}\right)\right)}^2=\rho \sum _{\mathbf{G},{\mathbf{G}}^{\prime }}{e}^{-\left(\frac{G^2}{4\alpha }+\frac{G^{\prime 2}}{4\alpha }\right)}\int {\mathrm{d}}^{3}r{e}^{i(\mathbf{G}+{\mathbf{G}}^{\prime })·\mathbf{r}}=N\sum _{\mathbf{G}}{e}^{-\frac{G^2}{2\alpha }}$$

(26)

(the above result is nothing but Parseval’s theorem as applied to (21)). The sum in the rhs of Equation (26) is the three-dimensional analog of a Jacobi theta function (see, e.g., [31]), whose value is $\mathcal{O}\left(1\right)$ for $\alpha >0$. Therefore, it follows from Equations (25) and (26) that the one-body entropy is at most $\mathcal{O}\left(N\right)$.

3.2. Two-Body Entropy

We now move to the problem of evaluating the two-body entropy $S_2$ for a crystal. For a homogeneous fluid, $S_2$ is an extensive quantity which, in $k_B$ units, is equal to

$$\mathrm{fluid}:S_2=-2\pi \rho N{\int }_0^{\infty }\mathrm{d}r{r}^2\left(g\left(r\right)lng\left(r\right)-g\left(r\right)+1\right).$$

(27)

For a crystal, we have from Equation (20) that

$$S_2=-\frac{1}{2}{\rho }^2\int {\mathrm{d}}^3{r}_1{\mathrm{d}}^3{r}_2{P}_1{P}_2\left(g_{12}ln{g}_{12}-g_{12}+1\right).$$

(28)

As $xlnx\ge x-1$ for $x>0$, $S_2$ is usually negative and zero exclusively for $g_{12}=1$. In terms of density functions, $S_2$ is written as

$$S_2=-\frac{1}{2}\int {\mathrm{d}}^3{r}_1{\mathrm{d}}^3{r}_2\left({\rho }^{\left(2\right)}({\mathbf{r}}_1,{\mathbf{r}}_2)ln\frac{{\rho }^{\left(2\right)}({\mathbf{r}}_1,{\mathbf{r}}_2)}{{\rho }^{\left(1\right)}\left({\mathbf{r}}_1\right){\rho }^{\left(1\right)}\left({\mathbf{r}}_2\right)}-{\rho }^{\left(2\right)}({\mathbf{r}}_1,{\mathbf{r}}_2)+{\rho }^{\left(1\right)}\left({\mathbf{r}}_1\right){\rho }^{\left(1\right)}\left({\mathbf{r}}_2\right)\right).$$

(29)

We show below that Equation (29) can be expressed as a radial integral, i.e., in a way similar to the two-body entropy for a fluid.

We can assign a radial structure to crystals by appealing to a couple of functions introduced in [32], namely

$${\rho }^2\tilde{g}\left(r\right)=\int \frac{{\mathrm{d}}^3{r}_1}{V}\int \frac{{\mathrm{d}}^2\Omega }{4\pi }{\rho }^{\left(2\right)}({\mathbf{r}}_1,{\mathbf{r}}_1+\mathbf{r})$$

(30)

and

$${\rho }^2{\tilde{g}}_0\left(r\right)=\int \frac{{\mathrm{d}}^3{r}_1}{V}\int \frac{{\mathrm{d}}^2\Omega }{4\pi }{\rho }^{\left(1\right)}\left({\mathbf{r}}_1\right){\rho }^{\left(1\right)}({\mathbf{r}}_1+\mathbf{r}),$$

(31)

where the inner integrals are over the direction of r. For a homogeneous fluid, $\tilde{g}\left(r\right)=g\left(r\right)$ and ${\tilde{g}}_0\left(r\right)=1$. The authors of reference [32] have sketched the profile of $\tilde{g}\left(r\right)$ and ${\tilde{g}}_0\left(r\right)$ for a crystal; both functions show narrow peaks at neighbor positions in the lattice, with an extra peak at zero distance for ${\tilde{g}}_0\left(r\right)$, and the oscillations persist till large distances. The following sum rules hold (cf. Equation (8) for $n=1$):

$$4\pi \int \mathrm{d}r{r}^2\rho \tilde{g}\left(r\right)=\frac{1}{\rho }4\pi \int \frac{{\mathrm{d}}^3{r}_1}{V}{\rho }^{\left(1\right)}\left({\mathbf{r}}_1\right)\underset{\frac{N-1}{4\pi }}{\underbrace{\int \frac{{\mathrm{d}}^3{r}_2}{4\pi }\frac{{\rho }^{\left(2\right)}({\mathbf{r}}_1,{\mathbf{r}}_2)}{{\rho }^{\left(1\right)}\left({\mathbf{r}}_1\right)}}}=N-1$$

(32)

and

$$4\pi \int \mathrm{d}r{r}^2\rho {\tilde{g}}_0\left(r\right)=\frac{1}{\rho }4\pi \int \frac{{\mathrm{d}}^3{r}_1}{V}{\rho }^{\left(1\right)}\left({\mathbf{r}}_1\right)\underset{\frac{N}{4\pi }}{\underbrace{\int \frac{{\mathrm{d}}^3{r}_2}{4\pi }{\rho }^{\left(1\right)}\left({\mathbf{r}}_2\right)}}=N.$$

(33)

When the latter two formulae are rewritten as

$$4\pi \rho \int \mathrm{d}r{r}^2(\tilde{g}\left(r\right)-1)=-1\mathrm{and}4\pi \rho \int \mathrm{d}r{r}^2({\tilde{g}}_0\left(r\right)-1)=0,$$

(34)

it becomes apparent that both $\tilde{g}\left(r\right)$ and ${\tilde{g}}_0\left(r\right)$ decay to 1 at infinity. Similarly, we define:

$${\rho }^2\tilde{h}\left(r\right)=\int \frac{{\mathrm{d}}^3{r}_1}{V}\int \frac{{\mathrm{d}}^2\Omega }{4\pi }{\rho }^{\left(2\right)}({\mathbf{r}}_1,{\mathbf{r}}_1+\mathbf{r})ln\frac{{\rho }^{\left(2\right)}({\mathbf{r}}_1,{\mathbf{r}}_1+\mathbf{r})}{{\rho }^{\left(1\right)}\left({\mathbf{r}}_1\right){\rho }^{\left(1\right)}({\mathbf{r}}_1+\mathbf{r})},$$

(35)

which obviously vanishes at infinity. While $\tilde{h}\left(r\right)=g\left(r\right)lng\left(r\right)$ for a homogeneous fluid, we expect that $\tilde{h}\left(r\right)\ne \tilde{g}\left(r\right)ln\tilde{g}\left(r\right)$ in the crystal. Putting Equations (30)–(35) together, we arrive at

$$\mathrm{crystal}:S_2=-2\pi \rho N{\int }_0^{\infty }\mathrm{d}r{r}^2\left(\tilde{h}\left(r\right)-\tilde{g}\left(r\right)+{\tilde{g}}_0\left(r\right)\right)=-2\pi \rho N{\int }_0^{\infty }\mathrm{d}r{r}^2\tilde{h}\left(r\right)-\frac{N}{2}.$$

(36)

Even though the integrand vanishes at infinity, $S_2=\mathcal{O}\left(N\right)$ only if the envelope of $\tilde{h}\left(r\right)$ decays faster than $r^{-3}$ ($r^{-2}$ in two dimensions). A slower decay may be sufficient if $S_2$ is computed through the first integral in (36). For a spherically-symmetric interaction potential, the excess energy (i.e., the canonical average of the total potential energy U) can also be written as a radial integral:

$$\begin{array}{ccc}\hfill 〈U〉& =& \frac{1}{2}\int {\mathrm{d}}^3{r}_1{\mathrm{d}}^3{r}_2{\rho }^{\left(2\right)}({\mathbf{r}}_1,{\mathbf{r}}_2)u\left(\right|{\mathbf{r}}_2-{\mathbf{r}}_1\left|\right)\hfill \\ & =& \frac{1}{2}{\int }_0^{\infty }\mathrm{d}r{r}^{2}u\left(r\right)\int {\mathrm{d}}^3{r}_1\int {\mathrm{d}}^2\Omega {\rho }^{\left(2\right)}({\mathbf{r}}_1,{\mathbf{r}}_1+\mathbf{r})=2\pi \rho N{\int }_0^{\infty }\mathrm{d}r{r}^{2}u\left(r\right)\tilde{g}\left(r\right).\hfill \end{array}$$

(37)

For the one-body density in (21), ${\tilde{g}}_0\left(r\right)$ can be obtained in closed form. First we have:

$$\int \frac{{\mathrm{d}}^2\Omega }{4\pi }{\rho }^{\left(1\right)}({\mathbf{r}}_1+\mathbf{r})=\rho \sum _{\mathbf{G}}{e}^{-\frac{G^2}{4\alpha }}{e}^{i\mathbf{G}·{\mathbf{r}}_1}\frac{sin\left(Gr\right)}{Gr}.$$

(38)

Then, multiplying by ${\rho }^{\left(1\right)}\left({\mathbf{r}}_1\right)=\rho {\sum }_{{\mathbf{G}}^{\prime }}{e}^{-\frac{G^{\prime 2}}{4\alpha }}{e}^{i{\mathbf{G}}^{\prime }·{\mathbf{r}}_1}$ and finally integrating over ${\mathbf{r}}_1$ we arrive at

$${\tilde{g}}_0\left(r\right)=1+\sum _{\mathbf{G}\ne 0}{e}^{-\frac{G^2}{2\alpha }}\frac{sin\left(Gr\right)}{Gr}.$$

(39)

We see that the large-distance decay of ${\tilde{g}}_0\left(r\right)$ is usually slow, and the same will occur for $\tilde{g}\left(r\right)$ since $\tilde{g}\left(r\right)\simeq {\tilde{g}}_0\left(r\right)$ for large r. In two dimensions, the one-body density and ${\tilde{g}}_0$ functions respectively read:

$${\rho }^{\left(1\right)}\left(\mathbf{r}\right)=\frac{\alpha }{\pi }\sum _{\mathbf{R}}{e}^{-\alpha {(\mathbf{r}-\mathbf{R})}^2}=\rho \sum _{\mathbf{G}}{e}^{-\frac{G^2}{4\alpha }}{e}^{i\mathbf{G}·\mathbf{r}}\mathrm{and}{\tilde{g}}_0\left(r\right)=1+\sum _{\mathbf{G}\ne 0}{e}^{-\frac{G^2}{2\alpha }}{J}_0\left(Gr\right),$$

(40)

where $J_0$ is a Bessel function of the first kind. Since the envelope of $J_0$ maxima decays as $r^{-1/2}$ at infinity, we see that the asymptotic vanishing of ${\tilde{g}}_0$ is slower in two dimensions than in three.

Equation (39) has a definite limit for $\alpha \to \infty$, corresponding to zero temperature. Indeed, using Poisson summation formula and the expression of Dirac’s delta in spherical coordinates, we obtain:

$$\begin{array}{ccc}\hfill \rho {\tilde{g}}_0\left(r\right)& =& \rho \left(1+\sum _{\mathbf{G}\ne 0}\frac{sin\left(Gr\right)}{Gr}\right)=\rho \int \frac{{\mathrm{d}}^2\Omega }{4\pi }\sum _{\mathbf{G}}{e}^{i\mathbf{G}·\mathbf{r}}={\delta }^3\left(\mathbf{r}\right)+\int \frac{{\mathrm{d}}^2\Omega }{4\pi }\sum _{\mathbf{R}\ne 0}{\delta }^3(\mathbf{r}-\mathbf{R})\hfill \\ & =& {\delta }^3\left(\mathbf{r}\right)+\sum _{\mathbf{R}\ne 0}\frac{1}{4\pi }{\int }_0^{2\pi }\mathrm{d}\varphi {\int }_0^{\pi }\mathrm{d}\theta sin\theta \frac{1}{r^{2}sin\theta }\delta (r-R)\delta (\theta -{\theta }_{\mathbf{R}})\delta (\varphi -{\varphi }_{\mathbf{R}})\hfill \\ & =& {\delta }^3\left(\mathbf{r}\right)+\sum _{\mathbf{R}\ne 0}\frac{1}{4\pi R^2}\delta (r-R).\hfill \end{array}$$

(41)

Hence, ${\tilde{g}}_0\left(r\right)$ reduces to a sum of delta functions centered at lattice distances (including the origin). The latter result is actually general. Inserting Equation (23) in (31), we obtain:

$$\begin{array}{ccc}\hfill \rho {\tilde{g}}_0\left(r\right)& =& \frac{1}{\rho }\int \frac{{\mathrm{d}}^3{r}_1}{V}\int \frac{{\mathrm{d}}^2\Omega }{4\pi }\sum _{\mathbf{R}}{\delta }^3({\mathbf{r}}_1-\mathbf{R})\sum _{{\mathbf{R}}^{\prime }}{\delta }^3({\mathbf{r}}_1+\mathbf{r}-{\mathbf{R}}^{\prime })\hfill \\ & =& \frac{1}{\rho }\int \frac{{\mathrm{d}}^3{r}_1}{V}\int \frac{{\mathrm{d}}^2\Omega }{4\pi }\left\{\sum _{\mathbf{R}}{\delta }^3({\mathbf{r}}_1-\mathbf{R}){\delta }^3({\mathbf{r}}_1+\mathbf{r}-\mathbf{R})+\sum _{\mathbf{R}\ne {\mathbf{R}}^{\prime }}{\delta }^3({\mathbf{r}}_1-\mathbf{R}){\delta }^3({\mathbf{r}}_1+\mathbf{r}-{\mathbf{R}}^{\prime })\right\}\hfill \\ & =& \frac{1}{\rho }\int \frac{{\mathrm{d}}^3{r}_1}{V}\int \frac{{\mathrm{d}}^2\Omega }{4\pi }\left\{\sum _{\mathbf{R}}{\delta }^3({\mathbf{r}}_1-\mathbf{R}){\delta }^3\left(\mathbf{r}\right)+\sum _{\mathbf{R}\ne {\mathbf{R}}^{\prime }}{\delta }^3({\mathbf{r}}_1-\mathbf{R}){\delta }^3({\mathbf{r}}_1+\mathbf{r}-{\mathbf{R}}^{\prime })\right\}\hfill \\ & =& \frac{1}{\rho }\sum _{\mathbf{R}}{\delta }^3\left(\mathbf{r}\right)\int \frac{{\mathrm{d}}^3{r}_1}{V}{\delta }^3({\mathbf{r}}_1-\mathbf{R})+\frac{1}{\rho }\sum _{\mathbf{R}\ne {\mathbf{R}}^{\prime }}\int \frac{{\mathrm{d}}^3{r}_1}{V}{\delta }^3({\mathbf{r}}_1-\mathbf{R})\frac{\delta (r-|{\mathbf{r}}_1-{\mathbf{R}}^{\prime }\left|\right)}{4\pi |{\mathbf{r}}_1-{\mathbf{R}}^{\prime }{|}^2}\hfill \\ & =& {\delta }^3\left(\mathbf{r}\right)+\sum _{\mathbf{R}\ne 0}\frac{1}{4\pi R^2}\delta (r-R),\hfill \end{array}$$

(42)

q.e.d. At zero temperature, $\rho \tilde{g}\left(r\right)$ is given by the same sum of delta-function terms as in (42), but for the first term, ${\delta }^3\left(\mathbf{r}\right)$, which is missing—see Equation (92) below.

We add a final comment on possible alternative formulations of $\tilde{g}\left(r\right)$ for a crystal. One choice is to replace (30) with

$$\mathrm{option}\mathrm{B}:\rho \tilde{g}\left(r\right)=\int \frac{{\mathrm{d}}^3{r}_1}{V}\int \frac{{\mathrm{d}}^2\Omega }{4\pi }\frac{{\rho }^{\left(2\right)}({\mathbf{r}}_1,{\mathbf{r}}_1+\mathbf{r})}{{\rho }^{\left(1\right)}\left({\mathbf{r}}_1\right)}.$$

(43)

Apparently, this is a good definition since (see Equation (8))

$$4\pi \int \mathrm{d}r{r}^2\rho \tilde{g}\left(r\right)=4\pi \underset{1}{\underbrace{\int \frac{{\mathrm{d}}^3{r}_1}{V}}}\underset{\frac{N-1}{4\pi }}{\underbrace{\int \frac{{\mathrm{d}}^3{r}_2}{4\pi }\frac{{\rho }^{\left(2\right)}({\mathbf{r}}_1,{\mathbf{r}}_2)}{{\rho }^{\left(1\right)}\left({\mathbf{r}}_1\right)}}}=N-1.$$

(44)

However, with this $\tilde{g}\left(r\right)$ we cannot write $S_2$ as a radial integral—hence, option B is discarded altogether. Another possibility is

$$\mathrm{option}\mathrm{C}:\tilde{g}\left(r\right)=\int \frac{{\mathrm{d}}^3{r}_1}{V}\int \frac{{\mathrm{d}}^2\Omega }{4\pi }\frac{{\rho }^{\left(2\right)}({\mathbf{r}}_1,{\mathbf{r}}_1+\mathbf{r})}{{\rho }^{\left(1\right)}\left({\mathbf{r}}_1\right){\rho }^{\left(1\right)}({\mathbf{r}}_1+\mathbf{r})},$$

(45)

but this option is useless too, since

$$4\pi \int \mathrm{d}r{r}^2\rho \tilde{g}\left(r\right)=\rho \int {\mathrm{d}}^3{r}_1\frac{1}{V}\int {\mathrm{d}}^3{r}_2{g}^{\left(2\right)}({\mathbf{r}}_1,{\mathbf{r}}_2)=?$$

(46)

(observe that the inner integral is different from the one appearing in Equation (8)).

3.3. Symmetries of the Two-Body Density

A general property of the two-body density for a crystal is the CE sum rule

$$\int {\mathrm{d}}^3{r}_2{\rho }^{\left(2\right)}({\mathbf{r}}_1,{\mathbf{r}}_2)=(N-1){\rho }^{\left(1\right)}\left({\mathbf{r}}_1\right).$$

(47)

Other constraints follow from the translational symmetry of local crystal properties. As for the one-body density, fulfilling ${\rho }^{\left(1\right)}({\mathbf{r}}_1+\mathbf{R})={\rho }^{\left(1\right)}\left({\mathbf{r}}_1\right)$ for every $\mathbf{R}$, we must have that

$${\rho }^{\left(2\right)}({\mathbf{r}}_1,{\mathbf{r}}_2)={\rho }^{\left(2\right)}({\mathbf{r}}_1+\mathbf{R},{\mathbf{r}}_2+\mathbf{R}),$$

(48)

in turn implying

$$g^{\left(2\right)}({\mathbf{r}}_1,{\mathbf{r}}_2)=g^{\left(2\right)}({\mathbf{r}}_1+\mathbf{R},{\mathbf{r}}_2+\mathbf{R}).$$

(49)

Now observe [33] that (i) any function of ${\mathbf{r}}_1$ and ${\mathbf{r}}_2$ can also be viewed as a function of $({\mathbf{r}}_1+{\mathbf{r}}_2)/2$ and ${\mathbf{r}}_2-{\mathbf{r}}_1$; (ii) under a $\mathbf{R}$-translation, only the former variable is affected, not the relative separation. Hence, the most general function consistent with (49) is:

$$g^{\left(2\right)}({\mathbf{r}}_1,{\mathbf{r}}_2)=\sum _{\mathbf{G}}{\tilde{v}}_{\mathbf{G}}({\mathbf{r}}_2-{\mathbf{r}}_1)e^{i\mathbf{G}·\frac{{\mathbf{r}}_1+{\mathbf{r}}_2}{2}},$$

(50)

where

$$g^{\left(2\right)}({\mathbf{r}}_1,{\mathbf{r}}_2)\in \mathbb{R}⟹{\tilde{v}}_{\mathbf{G}}^{*}({\mathbf{r}}_2-{\mathbf{r}}_1)={\tilde{v}}_{-\mathbf{G}}({\mathbf{r}}_2-{\mathbf{r}}_1)$$

(51)

and

$$g^{\left(2\right)}({\mathbf{r}}_1,{\mathbf{r}}_2)=g^{\left(2\right)}({\mathbf{r}}_2,{\mathbf{r}}_1)⟹{\tilde{v}}_{\mathbf{G}}({\mathbf{r}}_2-{\mathbf{r}}_1)={\tilde{v}}_{\mathbf{G}}({\mathbf{r}}_1-{\mathbf{r}}_2).$$

(52)

In order that ${lim}_{r\to \infty }{g}^{\left(2\right)}({\mathbf{r}}_1,{\mathbf{r}}_1+\mathbf{r})=1$ it is sufficient that

$$\underset{r\to \infty }{lim}{\tilde{v}}_0\left(\mathbf{r}\right)=1\mathrm{and}\underset{r\to \infty }{lim}{\tilde{v}}_{\mathbf{G}}\left(\mathbf{r}\right)=0\mathrm{for}\mathbf{G}\ne 0.$$

(53)

We may reasonably expect that the most relevant term in the expansion (50) is indeed the $\mathbf{G}=0$ one (also notice that ${\tilde{v}}_{\mathbf{G}}\to 0$ as $G\to \infty$ by the Riemann–Lebesgue lemma).

Equation (50) is still insufficient to establish the scaling of two-body entropy with the size of the crystal. Some general results can be obtained under the (strong) assumption that ${\tilde{v}}_{\mathbf{G}}\left(\mathbf{r}\right)=0$ for any $\mathbf{G}\ne 0$. If we change the notation from ${\tilde{v}}_0$ to $\mathcal{G}\left(\mathbf{r}\right)\equiv 1+\mathcal{H}\left(\mathbf{r}\right)$ (which, by Equations (51) and (52), is a real and even function), then a necessary condition for $\mathcal{H}$ is:

$$\int {\mathrm{d}}^3{r}_2{\rho }^{\left(1\right)}\left({\mathbf{r}}_2\right)\mathcal{H}({\mathbf{r}}_2-{\mathbf{r}}_1)=-1\mathrm{for}\mathrm{any}{\mathbf{r}}_1\mathrm{where}{\rho }^{\left(1\right)}\left({\mathbf{r}}_1\right)\ne 0.$$

(54)

The rationale behind Equation (54) is particularly transparent near $T=0$, where the peaks of the one-body density are extremely narrow. As argued below (see Equation (67) ff.), $\mathcal{H}$ as a function of ${\mathbf{r}}_2$ is roughly $-1$ in the primitive cell $\mathcal{C}$ centered in ${\mathbf{r}}_1\approx {\mathbf{R}}_1$, ${\mathbf{R}}_1$ denoting the only lattice site contained in $\mathcal{C}$ and roughly zero outside $\mathcal{C}$. Since the integral of ${\rho }^{\left(1\right)}$ over $\mathcal{C}$ equals 1, Equation (54) will immediately follow.

Now writing $\mathcal{H}\left(\mathbf{r}\right)$ as a Fourier integral,

$$\mathcal{H}\left(\mathbf{r}\right)=\int \frac{{\mathrm{d}}^{3}k}{{\left(2\pi \right)}^3}\tilde{\mathcal{H}}\left(\mathbf{k}\right)e^{i\mathbf{k}·\mathbf{r}},$$

(55)

and using (21) as one-body density, Equation (54) yields

$$\rho \sum _{\mathbf{G}}{e}^{-\frac{G^2}{4\alpha }}\tilde{\mathcal{H}}\left(\mathbf{G}\right)e^{-i\mathbf{G}·{\mathbf{r}}_1}=-1,$$

(56)

which can only hold for arbitrary ${\mathbf{r}}_1$ if

$$\tilde{\mathcal{H}}\left(\mathbf{G}\right)=-\frac{1}{\rho }{\delta }_{\mathbf{G},0}.$$

(57)

Next, from Equation (30) we obtain:

$${\rho }^2\tilde{g}\left(r\right)={\rho }^2{\tilde{g}}_0\left(r\right)+\int \frac{{\mathrm{d}}^3{r}_1}{V}{\rho }^{\left(1\right)}\left({\mathbf{r}}_1\right)\int \frac{{\mathrm{d}}^2\Omega }{4\pi }{\rho }^{\left(1\right)}({\mathbf{r}}_1+\mathbf{r})\mathcal{H}\left(\mathbf{r}\right).$$

(58)

For the one-body density in (21), the inner integral becomes:

$$\int \frac{{\mathrm{d}}^2\Omega }{4\pi }{\rho }^{\left(1\right)}({\mathbf{r}}_1+\mathbf{r})\mathcal{H}\left(\mathbf{r}\right)=\rho \sum _{\mathbf{G}}{e}^{-\frac{G^2}{4\alpha }}{I}_{\mathbf{G}}\left(r\right)e^{-i\mathbf{G}·{\mathbf{r}}_1}$$

(59)

with

$$I_{\mathbf{G}}\left(r\right)=\int \frac{{\mathrm{d}}^2\Omega }{4\pi }\mathcal{H}\left(\mathbf{r}\right)e^{-i\mathbf{G}·\mathbf{r}}=\int \frac{{\mathrm{d}}^{3}k}{{\left(2\pi \right)}^3}\tilde{\mathcal{H}}\left(\mathbf{k}\right)\frac{sin\left(|\mathbf{k}-\mathbf{G}|r\right)}{|\mathbf{k}-\mathbf{G}|r}.$$

(60)

It is evident that $I_{\mathbf{G}}\left(r\right)$ vanishes at infinity. Upon inserting (59) in (58), we finally obtain:

$$\tilde{g}\left(r\right)={\tilde{g}}_0\left(r\right)+\sum _{\mathbf{G}}{e}^{-\frac{G^2}{2\alpha }}{I}_{\mathbf{G}}\left(r\right).$$

(61)

As r increases, the second term gradually vanishes and the large-distance oscillations of $\tilde{g}\left(r\right)$ then exactly match those of ${\tilde{g}}_0\left(r\right)$. As a countercheck, let us compute the integral of $\rho \tilde{g}\left(r\right)-\rho {\tilde{g}}_0\left(r\right)$ over the macroscopic system volume (which, by Equations (32) and (33), should be $-1$):

$$\begin{array}{ccc}\hfill 4\pi \int \mathrm{d}r{r}^2\rho (\tilde{g}\left(r\right)-{\tilde{g}}_0\left(r\right))& =& \rho \sum _{\mathbf{G}}{e}^{-\frac{G^2}{2\alpha }}·4\pi \int \mathrm{d}r{r}^2{I}_{\mathbf{G}}\left(r\right)\hfill \\ & =& \rho \sum _{\mathbf{G}}{e}^{-\frac{G^2}{2\alpha }}\int \frac{{\mathrm{d}}^{3}k}{{\left(2\pi \right)}^3}\tilde{\mathcal{H}}\left(\mathbf{k}\right)·4\pi \int \mathrm{d}r{r}^2\frac{sin\left(|\mathbf{k}-\mathbf{G}|r\right)}{|\mathbf{k}-\mathbf{G}|r}\hfill \\ & =& \rho \sum _{\mathbf{G}}{e}^{-\frac{G^2}{2\alpha }}\int \frac{{\mathrm{d}}^{3}k}{{\left(2\pi \right)}^3}\tilde{\mathcal{H}}\left(\mathbf{k}\right)\underset{{\left(2\pi \right)}^3{\delta }^3(\mathbf{k}-\mathbf{G})}{\underbrace{\int {\mathrm{d}}^{3}r{e}^{i(\mathbf{k}-\mathbf{G})·\mathbf{r}}}}\hfill \\ & =& \rho \sum _{\mathbf{G}}{e}^{-\frac{G^2}{2\alpha }}\underset{-(1/\rho ){\delta }_{\mathbf{G},0}}{\underbrace{\tilde{\mathcal{H}}\left(\mathbf{G}\right)}}=-1.\hfill \end{array}$$

(62)

Under the assumption that

$${\rho }^{\left(2\right)}({\mathbf{r}}_1,{\mathbf{r}}_1+\mathbf{r})={\rho }^{\left(1\right)}\left({\mathbf{r}}_1\right){\rho }^{\left(1\right)}({\mathbf{r}}_1+\mathbf{r})\mathcal{G}\left(\mathbf{r}\right),$$

(63)

the entropy expansion for a crystal reads

$$\begin{array}{ccc}\hfill \frac{S}{k_B}=& & -\int {\mathrm{d}}^3{r}_1{\rho }^{\left(1\right)}\left({\mathbf{r}}_1\right)ln\frac{{\rho }^{\left(1\right)}\left({\mathbf{r}}_1\right)}{\rho }\hfill \\ & & -\frac{1}{2}\int {\mathrm{d}}^3{r}_1{\rho }^{\left(1\right)}\left({\mathbf{r}}_1\right)\int {\mathrm{d}}^{3}r{\rho }^{\left(1\right)}({\mathbf{r}}_1+\mathbf{r})\left[\mathcal{G}\left(\mathbf{r}\right)ln\mathcal{G}\left(\mathbf{r}\right)-\mathcal{G}\left(\mathbf{r}\right)+1\right]+\dots \hfill \end{array}$$

(64)

Providing that it vanishes sufficiently rapidly at infinity, the function

$$\mathcal{K}\left(\mathbf{r}\right)=\mathcal{G}\left(\mathbf{r}\right)ln\mathcal{G}\left(\mathbf{r}\right)-\mathcal{G}\left(\mathbf{r}\right)+1$$

(65)

can be written as a Fourier integral, and using (21) as one-body density, the two-body entropy becomes

$$S_2=-\frac{1}{2}{\rho }^2\sum _{\mathbf{G},{\mathbf{G}}^{\prime }}{e}^{-\frac{G^2+G^{\prime 2}}{4\alpha }}\underset{V{\delta }_{{\mathbf{G}}^{\prime },-\mathbf{G}}}{\underbrace{\int {\mathrm{d}}^3{r}_1{e}^{i(\mathbf{G}+{\mathbf{G}}^{\prime })·{\mathbf{r}}_1}}}\int {\mathrm{d}}^{3}r\mathcal{K}\left(\mathbf{r}\right)e^{i{\mathbf{G}}^{\prime }·\mathbf{r}}=-\frac{1}{2}N\rho \sum _{\mathbf{G}}{e}^{-\frac{G^2}{2\alpha }}\tilde{\mathcal{K}}\left(\mathbf{G}\right),$$

(66)

which is clearly $\mathcal{O}\left(N\right)$.

3.4. Two-Body Density at $T=0$

In the zero-temperature limit, particles will be sitting at lattice sites, and the two-body density then becomes (see Equation (23)):

$${\rho }^{\left(2\right)}({\mathbf{r}}_1,{\mathbf{r}}_2)=\underset{\mathbf{R},{\mathbf{R}}^{\prime }}{{\sum }^{\prime }}{\delta }^3({\mathbf{r}}_1-\mathbf{R}){\delta }^3({\mathbf{r}}_2-{\mathbf{R}}^{\prime })={\rho }^{\left(1\right)}\left({\mathbf{r}}_1\right){\rho }^{\left(1\right)}\left({\mathbf{r}}_2\right)\left(1-{\mathbf{1}}_{\mathcal{C}}({\mathbf{r}}_2-{\mathbf{r}}_1)\right),$$

(67)

which is of the form (63). In Equation (67), ${\mathbf{1}}_{\mathcal{C}}\left(\mathbf{r}\right)$ is the indicator function of a Wigner–Seitz cell $\mathcal{C}$ centered at the origin (i.e., ${\mathbf{1}}_{\mathcal{C}}\left(\mathbf{r}\right)=1$ if $\mathbf{r}\in \mathcal{C}$ and ${\mathbf{1}}_{\mathcal{C}}\left(\mathbf{r}\right)=0$ otherwise). While the factor ${\rho }^{\left(1\right)}\left({\mathbf{r}}_1\right){\rho }^{\left(1\right)}\left({\mathbf{r}}_2\right)$ forces particles to be located at lattice sites, the only role of the $\mathcal{G}$ in (67) is to prevent the possibility of double site occupancy. However, a $\mathcal{G}$ function with this property is not unique; the one provided in (67) has the advantage of exactly complying with condition (57) (see below). Equation (67) indicates that the pair-correlation structure of a low-temperature solid is very different from the structure of a dense fluid close to freezing.

For

$$\mathcal{H}\left(\mathbf{r}\right)=-{\mathbf{1}}_{\mathcal{C}}\left(\mathbf{r}\right)=\left\{\begin{array}{cc}\hfill -1,& \mathrm{for}\mathbf{r}\in \mathcal{C}\hfill \\ \hfill 0,& \mathrm{otherwise}\hfill \end{array}\right.$$

(68)

the Fourier transform reads:

$$\begin{array}{ccc}\hfill \tilde{\mathcal{H}}\left(\mathbf{k}\right)& =& \int {\mathrm{d}}^{3}r\mathcal{H}\left(\mathbf{r}\right)e^{-i\mathbf{k}·\mathbf{r}}=-{\int }_{\mathcal{C}}{\mathrm{d}}^{3}r{e}^{-i\mathbf{k}·\mathbf{r}}.\hfill \end{array}$$

(69)

Now observe that $f\left(\mathbf{r}\right)=1$ is trivially periodic, and can thus be expanded in plane waves as $1={\sum }_{\mathbf{G}}{\tilde{f}}_{\mathbf{G}}{e}^{i\mathbf{G}·\mathbf{r}}$, with ${\tilde{f}}_{\mathbf{G}}={\delta }_{\mathbf{G},0}$. On the other hand,

$${\tilde{f}}_{\mathbf{G}}=\frac{1}{v_0}{\int }_{\mathcal{C}}{\mathrm{d}}^{3}rf\left(\mathbf{r}\right)e^{-i\mathbf{G}·\mathbf{r}}=\rho {\int }_{\mathcal{C}}{\mathrm{d}}^{3}r{e}^{-i\mathbf{G}·\mathbf{r}}.$$

(70)

Comparing Equations (69) and (70), we conclude that

$$\tilde{\mathcal{H}}\left(\mathbf{G}\right)=-\frac{1}{\rho }{\delta }_{\mathbf{G},0}.$$

(71)

For $\mathcal{H}\left(\mathbf{r}\right)=-{\mathbf{1}}_{\mathcal{C}}\left(\mathbf{r}\right)$ the function $I_{\mathbf{G}}\left(r\right)$ at Equation (60) equals $-sin\left(Gr\right)/\left(Gr\right)$ for $r<r_m$ and 0 for $r>r_M$, where $r_m$ ($r_M$) is the radius of the largest (smallest) sphere inscribed in (circumscribed to) $\mathcal{C}$. It then follows from Equation (61) that $\tilde{g}\left(r\right)=0$ for $r<r_m$, while $\tilde{g}\left(r\right)={\tilde{g}}_0\left(r\right)$ for $r>r_M$ (for a triangular crystal with spacing a we have $r_m=a/2$ and $r_M=a/\sqrt{3}$, both comprised between the first, 0, and the second, a, lattice distance). $T=0$, where ${\tilde{g}}_0\left(r\right)$ consists of infinitely narrow peaks centered at lattice distances; this implies that $\tilde{g}\left(r\right)={\tilde{g}}_0\left(r\right)$ everywhere but at the origin, where $\tilde{g}\left(r\right)=0$ while ${\tilde{g}}_0\left(r\right)$ is non-zero.

3.5. Scaling of Two-Body Entropy with N

We henceforth discuss in fully general terms how the two-body entropy scales with N for a crystal, avoiding to make any simplifying hypothesis on the structure of $g^{\left(2\right)}({\mathbf{r}}_1,{\mathbf{r}}_2)$. Using an obvious short-hand notation, the two-body entropy reads

$$S_2=-\frac{1}{2}\int \mathrm{d}1\mathrm{d}2\left({\rho }_{12}ln\frac{{\rho }_{12}}{{\rho }_1{\rho }_2}-{\rho }_{12}+{\rho }_1{\rho }_2\right)=-\frac{1}{2}\int \mathrm{d}1\mathrm{d}2{\rho }_1{\rho }_2\left(g_{12}ln{g}_{12}-g_{12}+1\right).$$

(72)

As we already know, $S_2\le 0$. From the inequality $lnx\le x-1$, valid for all $x>0$, we derive $-xlnx\ge x-x^2$ for $x\ge 0$, and then obtain:

$$S_2=\frac{1}{2}\int \mathrm{d}1\mathrm{d}2{\rho }_1{\rho }_2\left(-g_{12}ln{g}_{12}+g_{12}-1\right)\ge -\frac{1}{2}\int \mathrm{d}1\mathrm{d}2{\rho }_1{\rho }_2{\left(g_{12}-1\right)}^2.$$

(73)

Clearly, estimating the size of the lower bound in Equation (73) is a much simpler problem than working with $S_2$ itself.

Taking $h_{12}\equiv g_{12}-1$, it is evident that $h_{12}$ shares all symmetries of $g_{12}$. Hence, we can write:

$$h_{12}=\sum _{\mathbf{G}}{\tilde{h}}_{\mathbf{G}}({\mathbf{r}}_1-{\mathbf{r}}_2)e^{i\mathbf{G}·\frac{{\mathbf{r}}_1+{\mathbf{r}}_2}{2}}\mathrm{with}{\tilde{h}}_{\mathbf{G}}^{*}\left(\mathbf{r}\right)={\tilde{h}}_{-\mathbf{G}}\left(\mathbf{r}\right)\mathrm{and}{\tilde{h}}_{\mathbf{G}}\left(\mathbf{r}\right)={\tilde{h}}_{\mathbf{G}}(-\mathbf{r}).$$

(74)

Observe that the ${\tilde{h}}_{\mathbf{G}}\left(\mathbf{r}\right)$ functions are nothing but Fourier coefficients, once the h function has been expressed in terms of $\mathbf{S}=({\mathbf{r}}_1+{\mathbf{r}}_2)/2$ and $\mathbf{r}={\mathbf{r}}_1-{\mathbf{r}}_2$:

$${\tilde{h}}_{\mathbf{G}}\left(\mathbf{r}\right)=\frac{1}{v_0}{\int }_{\mathcal{C}}{\mathrm{d}}^{3}Sh(\mathbf{S}+\mathbf{r}/2,\mathbf{S}-\mathbf{r}/2)e^{-i\mathbf{G}·\mathbf{S}}.$$

(75)

By the Riemann–Lebesgue lemma, ${\tilde{h}}_{\mathbf{G}}\left(\mathbf{r}\right)\to 0$ as $G\to \infty$ (for arbitrary $\mathbf{r}$). Moreover, ${\tilde{h}}_{\mathbf{G}}\left(\mathbf{r}\right)\to 0$ for $r\to \infty$ (for arbitrary $\mathbf{G}$) since ${\rho }_{12}\to {\rho }_1{\rho }_2$ for $|{\mathbf{r}}_1-{\mathbf{r}}_2|\to \infty$. Similarly, for $k_{12}\equiv h_{12}^2$ we have that

$$k_{12}=\sum _{\mathbf{G}}{\tilde{k}}_{\mathbf{G}}({\mathbf{r}}_1-{\mathbf{r}}_2)e^{i\mathbf{G}·\frac{{\mathbf{r}}_1+{\mathbf{r}}_2}{2}}\mathrm{with}{\tilde{k}}_{\mathbf{G}}\left(\mathbf{r}\right)=\sum _{{\mathbf{G}}^{\prime }}{\tilde{h}}_{\mathbf{G}-{\mathbf{G}}^{\prime }}\left(\mathbf{r}\right){\tilde{h}}_{{\mathbf{G}}^{\prime }}\left(\mathbf{r}\right).$$

(76)

Now observe that, for ${\rho }^{\left(1\right)}\left(\mathbf{r}\right)={\sum }_{\mathbf{G}}{\tilde{u}}_{\mathbf{G}}{e}^{i\mathbf{G}·\mathbf{r}}$,

$${\rho }^{\left(1\right)}\left(\mathbf{r}\right){\rho }^{\left(1\right)}\left({\mathbf{r}}^{\prime }\right)=\sum _{\mathbf{G}}\underset{{\tilde{v}}_{\mathbf{G}}^{\infty }(\mathbf{r}-{\mathbf{r}}^{\prime })}{\underbrace{\left(\sum _{{\mathbf{G}}^{\prime }}{\tilde{u}}_{\mathbf{G}-{\mathbf{G}}^{\prime }}{\tilde{u}}_{{\mathbf{G}}^{\prime }}{e}^{i(2{\mathbf{G}}^{\prime }-\mathbf{G})·\frac{\mathbf{r}-{\mathbf{r}}^{\prime }}{2}}\right)}}{e}^{i\mathbf{G}·\frac{\mathbf{r}+{\mathbf{r}}^{\prime }}{2}}.$$

(77)

Using the above equation, and changing the integration variables from ${\mathbf{r}}_1$ and ${\mathbf{r}}_2$ to $\mathbf{S}$ and $\mathbf{r}$, we obtain:

$$\begin{array}{ccc}\hfill -\frac{1}{2}\int \mathrm{d}1\mathrm{d}2{\rho }_1{\rho }_2{\left(g_{12}-1\right)}^2& =& -\frac{1}{2}\int {\mathrm{d}}^{3}S{\mathrm{d}}^{3}r\sum _{\mathbf{G}}{\tilde{v}}_{\mathbf{G}}^{\infty }\left(\mathbf{r}\right)e^{i\mathbf{G}·\mathbf{S}}\sum _{{\mathbf{G}}^{\prime }}{\tilde{k}}_{{\mathbf{G}}^{\prime }}\left(\mathbf{r}\right)e^{i{\mathbf{G}}^{\prime }·\mathbf{S}}\hfill \\ & =& -\frac{1}{2}\sum _{\mathbf{G},{\mathbf{G}}^{\prime }}\underset{V{\delta }_{{\mathbf{G}}^{\prime },-\mathbf{G}}}{\underbrace{\int {\mathrm{d}}^{3}S{e}^{i(\mathbf{G}+{\mathbf{G}}^{\prime })·\mathbf{S}}}}\underset{i_{\mathbf{G},{\mathbf{G}}^{\prime }}}{\underbrace{\int {\mathrm{d}}^{3}r{\tilde{v}}_{\mathbf{G}}^{\infty }\left(\mathbf{r}\right){\tilde{k}}_{{\mathbf{G}}^{\prime }}\left(\mathbf{r}\right)}}\hfill \\ & =& -\frac{1}{2}V\sum _{\mathbf{G}}{i}_{\mathbf{G},-\mathbf{G}},\hfill \end{array}$$

(78)

where

$$i_{\mathbf{G},-\mathbf{G}}=\sum _{{\mathbf{G}}^{\prime }}{\tilde{u}}_{\mathbf{G}-{\mathbf{G}}^{\prime }}{\tilde{u}}_{{\mathbf{G}}^{\prime }}\int {\mathrm{d}}^{3}r\left(\sum _{{\mathbf{G}}^{″}}{\tilde{h}}_{-\mathbf{G}-{\mathbf{G}}^{″}}\left(\mathbf{r}\right){\tilde{h}}_{{\mathbf{G}}^{″}}\left(\mathbf{r}\right)\right)e^{i(2{\mathbf{G}}^{\prime }-\mathbf{G})·\frac{\mathbf{r}}{2}}.$$

(79)

In the special case ${\tilde{h}}_{\mathbf{G}}=\mathcal{H}\left(\mathbf{r}\right){\delta }_{\mathbf{G},0}$, we have $h_{12}={\tilde{h}}_0=\mathcal{H}({\mathbf{r}}_1-{\mathbf{r}}_2)$ and ${\tilde{k}}_{\mathbf{G}}\left(\mathbf{r}\right)=\mathcal{H}{\left(\mathbf{r}\right)}^2{\delta }_{\mathbf{G},0}$. Then, from Equation (77) we derive

$$\sum _{\mathbf{G}}{i}_{\mathbf{G},-\mathbf{G}}=i_{0,0}=\int {\mathrm{d}}^{3}r{\tilde{v}}_0^{\infty }\left(\mathbf{r}\right){\tilde{k}}_0\left(\mathbf{r}\right)=\sum _{\mathbf{G}}{\left|{\tilde{u}}_{\mathbf{G}}\right|}^2\widetilde{{\mathcal{H}}^2}\left(\mathbf{G}\right).$$

(80)

An independent computation of the integral leads to the same result:

$$\begin{array}{ccc}\hfill -\frac{1}{2}\int \mathrm{d}1\mathrm{d}2{\rho }_1{\rho }_2{\left(g_{12}-1\right)}^2& =& -\frac{1}{2}\int {\mathrm{d}}^{3}r{\rho }^{\left(1\right)}\left(\mathbf{r}\right)\int {\mathrm{d}}^3{r}^{\prime }\underset{{\sum }_{\mathbf{G}}{\tilde{u}}_{\mathbf{G}}^{*}{e}^{-i\mathbf{G}·(\mathbf{r}+{\mathbf{r}}^{\prime })}}{\underbrace{{\rho }^{\left(1\right)}(\mathbf{r}+{\mathbf{r}}^{\prime })}}{\mathcal{H}}^2\left({\mathbf{r}}^{\prime }\right)\hfill \\ & =& -\frac{1}{2}\sum _{\mathbf{G}}{\tilde{u}}_{\mathbf{G}}^{*}\underset{V{\tilde{u}}_{\mathbf{G}}}{\underbrace{\int {\mathrm{d}}^{3}r{\rho }^{\left(1\right)}\left(\mathbf{r}\right)e^{-i\mathbf{G}·\mathbf{r}}}}\underset{\widetilde{{\mathcal{H}}^2}\left(\mathbf{G}\right)}{\underbrace{\int {\mathrm{d}}^3{r}^{\prime }{\mathcal{H}}^2\left({\mathbf{r}}^{\prime }\right)e^{-i\mathbf{G}·{\mathbf{r}}^{\prime }}}}\hfill \\ & =& -\frac{1}{2}V\sum _{\mathbf{G}}{\left|{\tilde{u}}_{\mathbf{G}}\right|}^2\widetilde{{\mathcal{H}}^2}\left(\mathbf{G}\right),\hfill \end{array}$$

(81)

which should be compared with Equation (66). For $\mathcal{H}\left(\mathbf{r}\right)=-{\mathbf{1}}_{\mathcal{C}}\left(\mathbf{r}\right)$ and ${\tilde{u}}_{\mathbf{G}}=\rho exp\{-G^2/\left(4\alpha \right)\}$, we readily obtain $S_2=-N/2$ from both Equations (66) and (78), meaning that in this case the two-body entropy coincides with its lower bound in Equation (73).

The quantity (80) is clearly $\mathcal{O}\left(1\right)$, since the summand is rapidly converging to zero; this implies that the two-body entropy of a crystal is, at least for ${\tilde{h}}_{\mathbf{G}}=\mathcal{H}\left(\mathbf{r}\right){\delta }_{\mathbf{G},0}$, bounded from below by a $\mathcal{O}\left(N\right)$ quantity. In the most general case, where Equations (78) and (79), rather, apply, we can only observe the following. As G grows in size, for any fixed ${\mathbf{G}}^{\prime }$ and ${\mathbf{G}}^{″}$ both ${\tilde{u}}_{\mathbf{G}-{\mathbf{G}}^{\prime }}{\tilde{u}}_{{\mathbf{G}}^{\prime }}$ and ${\tilde{h}}_{-\mathbf{G}-{\mathbf{G}}^{″}}\left(\mathbf{r}\right){\tilde{h}}_{{\mathbf{G}}^{″}}\left(\mathbf{r}\right)$ get smaller, suggesting that $i_{\mathbf{G},-\mathbf{G}}$ will decrease too. However, this is not enough to conclude that ${\sum }_{\mathbf{G}}{i}_{\mathbf{G},-\mathbf{G}}$ is $\mathcal{O}\left(1\right)$, and the only way to settle the problem is numerical.

3.6. Numerical Evaluation of the Structure Functions

The utility of (36) clearly relies on the possibility of determining the integrand in simulation with sufficient accuracy. First we see how the one-body entropy, Equation (24), is computed. We start dividing V into a large number $M=V/v_c$ of identical cubes of volume $v_c$, chosen to be small enough that a cube contains the center of at most one particle. Let $c_{\alpha }=0,1$ (with $\alpha =1,\dots ,M$) be the occupancy of the $\alpha$th cube in a given system configuration and $\langle c_{\alpha }\rangle$ its canonical average as computed in a long Monte Carlo simulation of the weakly constrained crystal (to fix the center of mass of the crystal in space it is sufficient to keep one particle fixed; then, periodic boundary conditions will contribute to keep crystalline axes also fixed in the course of simulation). Given this setup, the local density at ${\mathbf{r}}_1$ (a point inside the $\alpha$th cube) can be estimated as

$${\rho }^{\left(1\right)}\left({\mathbf{r}}_1\right)\approx \frac{〈c_{\alpha }〉}{v_c},$$

(82)

and the integral in (24) becomes

$$\int {\mathrm{d}}^3{r}_1{\rho }^{\left(1\right)}\left({\mathbf{r}}_1\right)ln\frac{{\rho }^{\left(1\right)}\left({\mathbf{r}}_1\right)}{\rho }\approx \sum _{\alpha =1}^M〈c_{\alpha }〉ln\frac{〈c_{\alpha }〉}{\rho v_c}$$

(83)

(notice that $\rho v_c=N/M\ll 1$; we need $v_c\to 0$ and an infinitely long simulation to make (82) an exact relation). Similarly, if ${\mathbf{r}}_2$ falls within the $\beta$th cube, then

$${\rho }^{\left(2\right)}({\mathbf{r}}_1,{\mathbf{r}}_2)\approx \frac{〈c_{\alpha }{c}_{\beta }〉}{v_c^2}$$

(84)

and from Equation (30) we derive

$$\begin{array}{ccc}\hfill \tilde{g}\left(r\right)& \approx & \frac{1}{{\rho }^{2}V}\sum _{\alpha =1}^M{v}_c\frac{1}{N_{\gamma }}\sum _{\left|\gamma \right|=r}\frac{〈c_{\alpha }{c}_{\alpha +\gamma }〉}{v_c^2}=\frac{1}{{\rho }^{2}V}〈\sum _{\alpha =1}^M{v}_c\frac{1}{N_{\gamma }}\sum _{\left|\gamma \right|=r}\frac{c_{\alpha }{c}_{\alpha +\gamma }}{v_c^2}〉\hfill \\ & =& \frac{1}{M{\rho }^2{v}_c^2}〈\sum _{\alpha =1}^M{\delta }_{c_{\alpha },1}\frac{1}{N_{\gamma }}\sum _{\left|\gamma \right|=r}{c}_{\alpha +\gamma }〉.\hfill \end{array}$$

(85)

In the above formula $N_{\gamma }\simeq 4\pi r^2\Delta r/v_c$ is the number of cubes whose center lies at a distance r from $\alpha$ (to within a certain tolerance $\Delta r\ll r$), and the inner sum is carried out over those cubes only. Since $\rho v_c{N}_{\gamma }=4\pi r^2\Delta r\rho$ and $M\rho v_c=N$, an equivalent formula for $\tilde{g}\left(r\right)$ is

$$\tilde{g}\left(r\right)\approx 〈\frac{1}{N}\sum _{i=1}^N\frac{{\mathcal{N}}_i(r\pm \Delta r/2)}{4\pi r^2\Delta r\rho }〉,$$

(86)

denoting ${\mathcal{N}}_i(r\pm \Delta r/2)$ the number of particles found at a distance between $r-\Delta r/2$ and $r+\Delta r/2$ from the ith particle in the given configuration. Equation (86) closely reflects the method of computing the radial distribution function in a CE simulation (see, e.g., Equation (11) in reference [34]).

The function $\tilde{g}\left(r\right)$ admits yet another expression, which further strengthens its resemblance to the $g\left(r\right)$ of a liquid (as reported e.g., in [35]). It follows from Equations (30) and (6) that

$${\rho }^2\tilde{g}\left(r\right)=\frac{1}{V}\int \frac{{\mathrm{d}}^2\Omega }{4\pi }〈\underset{ij}{{\sum }^{\prime }}\int {\mathrm{d}}^3{r}_1{\delta }^3({\mathbf{r}}_1-{\mathbf{R}}_i){\delta }^3({\mathbf{r}}_1+\mathbf{r}-{\mathbf{R}}_j)〉.$$

(87)

Note that, for any sufficiently smooth function $f\left(\mathbf{r}\right)$,

$$\begin{array}{ccc}\hfill \int {\mathrm{d}}^{3}rf\left(\mathbf{r}\right)\int {\mathrm{d}}^3{r}_1{\delta }^3({\mathbf{r}}_1-{\mathbf{R}}_i){\delta }^3({\mathbf{r}}_1+\mathbf{r}-{\mathbf{R}}_j)& =& \int {\mathrm{d}}^3{r}_1{\delta }^3({\mathbf{r}}_1-{\mathbf{R}}_i)\int {\mathrm{d}}^{3}rf\left(\mathbf{r}\right){\delta }^3({\mathbf{r}}_1+\mathbf{r}-{\mathbf{R}}_j)\hfill \\ & =& \int {\mathrm{d}}^3{r}_1{\delta }^3({\mathbf{r}}_1-{\mathbf{R}}_i)f({\mathbf{R}}_j-{\mathbf{r}}_1)=f({\mathbf{R}}_j-{\mathbf{R}}_i)\hfill \end{array}$$

(88)

and

$$\begin{array}{ccc}\hfill \int {\mathrm{d}}^{3}rf\left(\mathbf{r}\right)\int {\mathrm{d}}^3{r}_1{\delta }^3({\mathbf{r}}_1-{\mathbf{R}}_i){\delta }^3({\mathbf{R}}_i+\mathbf{r}-{\mathbf{R}}_j)& =& \int {\mathrm{d}}^3{r}_1{\delta }^3({\mathbf{r}}_1-{\mathbf{R}}_i)\int {\mathrm{d}}^{3}rf\left(\mathbf{r}\right){\delta }^3({\mathbf{R}}_i+\mathbf{r}-{\mathbf{R}}_j)\hfill \\ & =& f({\mathbf{R}}_j-{\mathbf{R}}_i)\underset{1}{\underbrace{\int {\mathrm{d}}^3{r}_1{\delta }^3({\mathbf{r}}_1-{\mathbf{R}}_i)}}=f({\mathbf{R}}_j-{\mathbf{R}}_i),\hfill \end{array}$$

(89)

we are allowed to replace ${\delta }^3({\mathbf{r}}_1-{\mathbf{R}}_i){\delta }^3({\mathbf{r}}_1+\mathbf{r}-{\mathbf{R}}_j)$ with ${\delta }^3({\mathbf{r}}_1-{\mathbf{R}}_i){\delta }^3({\mathbf{R}}_i+\mathbf{r}-{\mathbf{R}}_j)$ in Equation (87), and thus obtain

$${\rho }^2\tilde{g}\left(r\right)=\frac{1}{V}\int \frac{{\mathrm{d}}^2\Omega }{4\pi }〈\underset{ij}{{\sum }^{\prime }}{\delta }^3({\mathbf{R}}_i+\mathbf{r}-{\mathbf{R}}_j)\underset{1}{\underbrace{\int {\mathrm{d}}^3{r}_1{\delta }^3({\mathbf{r}}_1-{\mathbf{R}}_i)}}〉=\frac{1}{V}\int \frac{{\mathrm{d}}^2\Omega }{4\pi }〈\underset{ij}{{\sum }^{\prime }}{\delta }^3({\mathbf{R}}_i+\mathbf{r}-{\mathbf{R}}_j)〉,$$

(90)

which finally leads to

$$\rho \tilde{g}\left(r\right)=\int \frac{{\mathrm{d}}^2\Omega }{4\pi }〈\frac{1}{N}\sum _i\sum _{j\ne i}{\delta }^3({\mathbf{R}}_i+\mathbf{r}-{\mathbf{R}}_j)〉.$$

(91)

At zero temperature, we can neglect the average and simply write

$$\rho \tilde{g}\left(r\right)=\int \frac{{\mathrm{d}}^2\Omega }{4\pi }\frac{1}{N}\sum _i\sum _{j\ne i}{\delta }^3({\mathbf{R}}_i+\mathbf{r}-{\mathbf{R}}_j)=\int \frac{{\mathrm{d}}^2\Omega }{4\pi }\sum _{\mathbf{R}\ne 0}{\delta }^3(\mathbf{r}-\mathbf{R})=\sum _{\mathbf{R}\ne 0}\frac{1}{4\pi R^2}\delta (r-R),$$

(92)

where in the last step we have followed the same path leading to Equation (41).

We can similarly proceed for the functions at Equations (31) and (35), which can be computed by the following formulae:

$${\tilde{g}}_0\left(r\right)\approx \frac{1}{{\rho }^{2}V}\sum _{\alpha =1}^M{v}_c\frac{1}{N_{\gamma }}\sum _{\left|\gamma \right|=r}\frac{〈c_{\alpha }〉}{v_c}\frac{〈c_{\alpha +\gamma }〉}{v_c}=\frac{1}{M{\rho }^2{v}_c^2}\sum _{\alpha =1}^M〈c_{\alpha }〉\frac{1}{N_{\gamma }}\sum _{\left|\gamma \right|=r}〈c_{\alpha +\gamma }〉$$

(93)

and

$$\tilde{h}\left(r\right)\approx \frac{1}{M{\rho }^2{v}_c^2}\sum _{\alpha =1}^M\frac{1}{N_{\gamma }}\sum _{\left|\gamma \right|=r}〈c_{\alpha }{c}_{\alpha +\gamma }〉ln\frac{〈c_{\alpha }{c}_{\alpha +\gamma }〉}{〈c_{\alpha }〉〈c_{\alpha +\gamma }〉}.$$

(94)

While $\tilde{g}\left(r\right)$ is the statistical average of an estimator whose histogram can be updated in the course of the simulation (see Equation (86)), ${\tilde{g}}_0\left(r\right)$ can only be estimated at the end of simulation, once $\langle c_{\alpha }\rangle$ has been evaluated for every $\alpha$ with effort comparable to that made for the one-body entropy. Much more costly is the calculation of $\tilde{h}\left(r\right)$, which should also be performed at the end of simulation after evaluating $\langle c_{\alpha }{c}_{\beta }\rangle$ for every $\alpha$ and $\beta$.

Using translational lattice symmetry, the radial distribution functions and $\tilde{h}\left(r\right)$ of a crystal can also be written as:

$$\begin{array}{ccc}\hfill \rho \tilde{g}\left(r\right)& =& {\int }_{\mathcal{C}}{\mathrm{d}}^3{r}_1\int \frac{{\mathrm{d}}^2\Omega }{4\pi }{\rho }^{\left(2\right)}({\mathbf{r}}_1,{\mathbf{r}}_1+\mathbf{r});\rho {\tilde{g}}_0\left(r\right)={\int }_{\mathcal{C}}{\mathrm{d}}^3{r}_1\int \frac{{\mathrm{d}}^2\Omega }{4\pi }{\rho }^{\left(1\right)}\left({\mathbf{r}}_1\right){\rho }^{\left(1\right)}({\mathbf{r}}_1+\mathbf{r});\hfill \\ \hfill \rho \tilde{h}\left(r\right)& =& {\int }_{\mathcal{C}}{\mathrm{d}}^3{r}_1\int \frac{{\mathrm{d}}^2\Omega }{4\pi }{\rho }^{\left(2\right)}({\mathbf{r}}_1,{\mathbf{r}}_1+\mathbf{r})ln\frac{{\rho }^{\left(2\right)}({\mathbf{r}}_1,{\mathbf{r}}_1+\mathbf{r})}{{\rho }^{\left(1\right)}\left({\mathbf{r}}_1\right){\rho }^{\left(1\right)}({\mathbf{r}}_1+\mathbf{r})},\hfill \end{array}$$

(95)

leading to simplifying Equations (85), (93), and (94) into

$$\begin{array}{ccc}\hfill \tilde{g}\left(r\right)& =& 〈\sum _{\alpha =1}^{M/N}{\delta }_{c_{\alpha },1}\frac{{\sum }_{\left|\gamma \right|=r}{c}_{\alpha +\gamma }}{4\pi r^2\Delta r\rho }〉;{\tilde{g}}_0\left(r\right)=\sum _{\alpha =1}^{M/N}〈c_{\alpha }〉\frac{{\sum }_{\left|\gamma \right|=r}〈c_{\alpha +\gamma }〉}{4\pi r^2\Delta r\rho };\hfill \\ \hfill \tilde{h}\left(r\right)& =& \sum _{\alpha =1}^{M/N}\frac{1}{4\pi r^2\Delta r\rho }\sum _{\left|\gamma \right|=r}〈c_{\alpha }{c}_{\alpha +\gamma }〉ln\frac{〈c_{\alpha }{c}_{\alpha +\gamma }〉}{〈c_{\alpha }〉〈c_{\alpha +\gamma }〉}.\hfill \end{array}$$

(96)

In the above formulae, the $\alpha$ index only runs over the cubes contained in a Wigner–Seitz/Voronoi cell of the lattice, while the $\beta$ sum is still carried out over all cubes in the simulation box.

3.7. Numerical Tests

We first examine the shape of the structure functions $\tilde{g}\left(r\right)$ and ${\tilde{g}}_0\left(r\right)$ for hard spheres, choosing a r resolution of $\Delta r=0.05$ (in units of the particle diameter $\sigma$). We take a system of $N=4000$ particles arranged in a fcc lattice with packing fraction $\eta =0.600$ (recall that the melting value is approximately 0.545). Periodic conditions are applied at the system boundary. In order to constrain the crystal in space, we keep one particle fixed during the simulation. As for ${\tilde{g}}_0\left(r\right)$, we employ the Tarazona ansatz for $\alpha =95$ (see Equation (39)), a value providing the best fit to the one-body density drawn from simulation.

We use the standard Metropolis Monte Carlo (MC) algorithm, constantly adjusting the maximum shift of a particle during equilibration until the fraction of accepted moves becomes close to $50\%$ (then, the maximum shift is no longer changed). We produce 50,000 MC cycles in the equilibration run, whereas CE averages are computed over a total of further $2\times {10}^5$ cycles. Our results are plotted in Figure 1. While at short distances $\tilde{g}\left(r\right)$ and ${\tilde{g}}_0\left(r\right)$ are rather different, as r increases the oscillations of the two functions become closer and closer in amplitude.

To obtain the one-body density with sufficient accuracy, we use a grid of about 50 points along each space direction in the unit cell. However, this grid resolution is too high for allowing the computation of $\tilde{h}\left(r\right)$, as the memory requirements for processing the $〈c_{\alpha }{c}_{\beta }〉$ data are very huge. On the other hand, a coarser grid is incompatible with the $\Delta r$ chosen.

To get closer to achieving our goal, i.e., to ascertain the N dependence of the two-body entropy for a crystal, we consider a two-dimensional system—hard disks. For this system, the transformation from fluid to solid occurs in two stages, via an intermediate hexatic fluid phase [36] (the transition from isotropic to hexatic fluid is first-order, whereas the hexatic-solid transition is continuous and occurs at $\eta =0.700$). We consider a system of $N=1152$ hard disks, arranged in a triangular crystal with packing fraction $\eta =0.800$, and a mesh consisting of about 80 points along each direction in the unit cell. Even though translational correlations are only quasi-long-ranged in an infinite two-dimensional crystal, when one of the particles is kept artificially fixed this specificity is lost and the (finite) two-dimensional crystal is made fully similar to a three-dimensional crystal. Observe also that an infinite two-dimensional crystal shares at least the same breaking of rotational symmetry typical of an infinite three-dimensional crystal.

As before, we first look at the structure functions drawn from simulation, $\tilde{g}\left(r\right)$ and ${\tilde{g}}_0\left(r\right)$. Our results are plotted in Figure 2, together with the ${\tilde{g}}_0\left(r\right)$ function of Equation (40) for $\alpha =75$. For this $\alpha$ the matching between the two ${\tilde{g}}_0$ functions is nearly perfect, indicating that the peaks of the one-body density are (to a high level of accuracy) Gaussian in shape. For $\eta =0.800$ we find $S_1/N=-2.156$.

In Figure 3, we show our main result, $\tilde{h}\left(r\right)$, for $\eta =0.800$ and two different crystal sizes, $N=288$ and 1152. We point out that, in order to obtain these data, we had to run a separate simulation for each r, as the memory usage is rather extreme. To be sure, we have computed the ${\tilde{g}}_0$ values in an independent way, i.e., using the same program loop written for $\tilde{h}\left(r\right)$, eventually finding the same results as in Figure 2. Looking at Figure 3, we see that $\tilde{h}\left(r\right)$ shows a series of peaks at neighbor positions and in the valleys within, taking preferentially positive values (meaning that its oscillations are not centered around zero). However, the damping of large-distance oscillations is too gradual to allow us to assess the nature of the asymptotic decay of $\tilde{h}\left(r\right)$ and then compute $S_2$. We must attempt a few explanations for this behavior of $\tilde{h}\left(r\right)$: On one hand, the decay of $\tilde{h}\left(r\right)$ may really be slow (at least in two dimensions), but $S_2$ would nonetheless be extensive, which implies a large $S_2/N$ value. It may as well be that constraining the crystal in space by hinging the position of one particle has a strong effect on the speed of $\tilde{h}$ decay, which only a finite-size scaling of data can relieve. Indeed, when going from $N=288$ to $N=1152$ the values of $\tilde{h}$ are slightly shifted downwards.

In summary, we have not reached any clear demonstration of $S_2$ extensivity in a crystal. This task has proved to be very hard to settle numerically. Our hope is that, based on our preparatory work, other authors with more powerful computational resources at their disposal can push the numerical analysis forward and eventually come up with a definite solution of the problem.

4. Conclusions

In this paper, we inquired into the possibility of extending the zero-RMPE criterion, a popular one-phase criterion of freezing for simple fluids, to also cover the melting of a solid. After revisiting the derivation of the entropy MPCE in the canonical ensemble, we argued that the formula applies for a crystal too. We exploited lattice symmetries to constrain the structures of one- and two-body densities, so as to gain as much information as possible on the first few terms in the entropy expansion. While this was enough to prove that the crystal one-body entropy is an extensive quantity, the information obtained was not sufficient to hold the same for the two-body entropy, whose scaling with the size of the crystal remains elusive. We thus attempted to clarify the question numerically, but we faced an insurmountable obstacle in the computational and memory limitations. To alleviate the problem, we turned towards a two-dimensional case, namely, hard disks, but with poor results: the structure function that must be integrated over distances to obtain the two-body entropy is weakly convergent to zero. In the near future, we plan to check whether the situation is more favorable for a different two-dimensional interaction, either endowed with an attractive tail (e.g., the Lennard-Jones potential) or provided with a soft core (for example, a Gaussian repulsion).