Ensemble averaged Madelung energies of finite volumes and surfaces

The cohesive energy of ionic crystals is dominated by the electrostatic energy between ionic point charges, known as the Madelung energy. The calculation of Madelung energies is a mathematically non-trivial problem because of the long-range nature of the Coulomb interaction. The Madelung constant of basic crystal structures was first successfully calculated by Ewald [1]. In the Ewald method, the Coulomb interaction is divided into a short range part for which the Madelung sum converges fast in real space and a long range part which can be summed in reciprocal space. The Ewald method is very accurate and widely used, but it is numerically rather involved and relies on periodic boundary conditions. Alternatively, the Madelung energy can be calculated through direct summation in real space, which is numerically simpler and can be used in finite system [2] and non-periodic structures. However, the lattice sums are only conditionally convergent, i.e. the result depends on the summation order. This reflects the physical fact that in a finite crystal, the potential at an inner site can be changed at will by choosing particular surface terminations [3]. In three dimensions, the Madelung sums diverge for the most natural, shell-like summation order [4]. Divergence can be avoided in two ways. In the first type of methods, the lattice is divided into neutral cells of vanishing dipole moment [5]. The sum over these cells converges absolutely because the quadrupole-quadrupole interaction decays as $1/r^{5}$. However, dipole moment free cells generally involve fractional ions at the corners and edges, and can be difficult to construct for low symmetry systems [6]. In the second type of method, the summation is done in the natural order of increasing distance but the system is neutralized at each step by adding a background or surrounding sphere [7, 8]. An example is the Wolf method [7] which keeps only the short-range part of the Ewald method but compensates each charge inside the summation sphere by an opposite charge at the cut-off radius.

Here we consider ensemble averaged quantities in finite subvolumes of a macroscopic system, especially the mean electrostatic potential at the sites of a given ionic species. This leads to an alternative definition of Madelung energies which converges unconditionally as a function of system size for any subvolume shape. We derive an exact expression of the Madelung constants in finite spheres. The leading term in the expansion over inverse size is found to be independent of the crystal structure. As a consequence, the Madelung contribution to the surface energy, averaged over surface orientations, is the same for all crystal structures, and provides a universal first order approximation of the surface energy of ionic systems. We apply the theory to the relative stability of CsCl and NaCl ionic structures as a function of system size. Under the assumption of equal nearest neighbor distance we find that the stable structure switches from NaCl to CsCl when the system size exceeds a few hundred ions, which explains the phase transition observed in cesium halide clusters.

We consider a collection of $N$ point charges at positions ${\bf r}_{i}^{\alpha}$, where $\alpha$ labels the different species with charge $q_{\alpha}$. The electrostatic energy is given by

$$U=\frac{1}{2}\sum_{j\beta\neq i\alpha}\frac{q_{\alpha}q_{\beta}}{|{\bf r}^{\alpha}_{i}-{\bf r}^{\beta}_{j}|}=\frac{1}{2}\sum_{\alpha}q_{\alpha}\sum_{i}\phi({\bf r}^{\alpha}_{i})$$

(1)

where

$$\phi({\bf r}^{\alpha}_{i})=\sum_{j\beta(\neq i\alpha)}\frac{q_{\beta}}{|{\bf r}^{\alpha}_{i}-{\bf r}^{\beta}_{j}|}\;.$$

(2)

is the electrostatic potential at site ${\bf r}^{\alpha}_{i}$. For crystals we take each ion in the unit cell as a different species $\alpha$ and $i$ is the cell index. In the infinite crystal, $\phi({\bf r}^{\alpha}_{i})$ is independent of $i$. The Madelung constant of the species $\alpha$ is commonly defined as

$$M_{\alpha}=-\phi({\bf r}^{\alpha}_{i}){d}/{q_{\alpha}}\;,$$

(3)

where $d$ is the nearest neighbor distance. In a practical real space summation, one computes the potential $\phi({\bf r}^{\alpha}_{i})$ at some site $i$ of a finite cluster (usually the central site), and lets the cluster size $N$ go to infinity. The problem with this approach is that the series in Eq. (2), converges only conditionally for $N\rightarrow\infty$. This means that the sum, i.e. the potential at site $i$, depends on the summation order and may diverge.

Instead of the potential at a given site $i$, we consider the average potential at the sites of the $\alpha$-ions,

$${\bar{\phi}}_{\alpha}=\frac{1}{N_{\alpha}}\sum_{ij\beta}^{i\neq j}\frac{q_{\beta}}{|{\bf r}^{\alpha}_{i}-{\bf r}^{\beta}_{j}|}$$

(4)

where $N_{\alpha}$ is the number of $\alpha$ ions in the cluster. The electrostatic energy in Eq. (1) can be rewritten as

$$U=\frac{1}{2}\sum_{\alpha}N_{\alpha}q_{\alpha}{\bar{\phi}}_{\alpha}\;.$$

(5)

Equations (1) and (5) are equivalent and hold for finite and infinite systems. For an infinite crystal, we obviously have ${\bar{\phi}}_{\alpha}=\phi({\bf r_{i}^{\alpha}})$. This suggests a redefinition of the Madelung constants as

$$M_{\alpha}=-{\bar{\phi}}_{\alpha}{d}/{q_{\alpha}}\;.$$

(6)

Of course, as with Eqs. (2,3), the true Madelung constants are obtained in the limit $N\rightarrow\infty$. The advantage of definition (6) is that the electrostatic energy can be expressed simply in terms of $M_{\alpha}$ as

$$U=-\frac{1}{2d}\sum_{\alpha}N_{\alpha}M_{\alpha}q_{\alpha}^{2}\;,$$

(7)

which holds for any system, finite or infinite.

We introduce the particle density of species $\alpha$

$$\rho_{\alpha}({\bf r})=\langle\sum_{i}\delta({\bf r}-{\bf r}_{i}^{\alpha})\rangle$$

(8)

where $\langle\dots\rangle$ denotes a statistical ensemble average. The pair distribution function is defined as

$$g_{\alpha\beta}({\bf r},{\bf r}^{\prime})=\frac{1}{\rho_{\alpha}({\bf r})\rho_{\beta}({\bf r}^{\prime})}\langle\sum_{ij}^{i\neq j}\delta({\bf r}-{\bf r}_{i}^{\alpha})\delta({\bf r}^{\prime}-{\bf r}_{j}^{\beta})\rangle\;.$$

(9)

With these functions, we can rewrite the ensemble average of the Eq. (4), as

$$\langle{\bar{\phi}}_{\alpha}\rangle=\frac{1}{N_{\alpha}}\sum_{\beta}q_{\beta}\int d{\bf r}\rho_{\alpha}({\bf r})\int d{\bf r}^{\prime}\rho_{\beta}({\bf r}^{\prime})\frac{g_{\alpha\beta}({\bf r},{\bf r}^{\prime})}{|{\bf r}-{\bf r}^{\prime}|}\;.$$

(10)

Now we consider a finite subvolume $V$ of an infinite ionic solid and average over all possible positions and orientations of $V$. Equivalently we can keep the subvolume fixed in space and perform statistical averaging in an ensemble of crystals of random position and orientation, which is realized in a powder sample. This ensemble is homogeneous and isotropic, i.e. it has the symmetry of a fluid, where the density $\rho_{\alpha}$ is a constant, and the pair distribution function $g_{\alpha\beta}(r)$ depends only on the distance $r=|{\bf r}_{i}^{\alpha}-{\bf r}_{j}^{\beta}|$ [9]. As a result, and using $\rho_{\alpha}=N_{\alpha}/V$, Eq. (10) simplifies to

$$\langle{\bar{\phi}}_{\alpha}\rangle=\frac{1}{V}\sum_{\beta}\rho_{\beta}q_{\beta}\int_{V}d{\bf r}\int_{V}d{\bf r}^{\prime}\frac{g_{\alpha\beta}(|{\bf r}-{\bf r}^{\prime}|)}{|{\bf r}-{\bf r}^{\prime}|}\;.$$

(11)

This expression holds for a finite subvolume $V$ of any homogeneous and isotropic ensemble, in particular for a fluid or powder sample. Further, as explained above, statistical averaging in such an ensemble is equivalent to considering a single crystal and averaging over the position and orientation of the subvolume $V$. Hence Eq. (11) also holds in this sense. In standard direct space lattice summation methods, a finite cluster size corresponds to a subvolume of the infinite crystal with fixed center and orientation as illustrated in Fig. 1 (a). In the present method, an average is done over all possible positions and orientations (Fig. 1 b).

In the special case of a spherical subvolume considered in the numerical applications below, only positional averaging is meaningful. However, the main results and conclusions are independent of the subvolume shape.

For stability, any ionic system must be globally charge neutral, i.e.

$$\sum_{\alpha}\rho_{\alpha}q_{\alpha}=0\;.$$

(12)

We do not assume any form of local charge neutrality. Averaging is done over all configurations in the ensemble, including those where the finite volume $V$ has a net charge. By virtue of Eq. (12), we may add to the integrand in Eq. (11) any function $f({\bf r},{\bf r}^{\prime})$ that is independent of $\alpha$, without changing the value of the expression. We shall add $-1/|{\bf r}-{\bf r}^{\prime}|$, i.e. we replace $g_{\alpha\beta}(r)$ by the pair correlation function $h_{\alpha\beta}(r)\equiv g_{\alpha\beta}(r)-1$ and obtain

$$\langle{\bar{\phi}}_{\alpha}\rangle=\sum_{\beta}\rho_{\beta}q_{\beta}H_{\alpha\beta}^{V}\;,$$

(13)

where

$$H_{\alpha\beta}^{V}=\frac{1}{V}\int_{V}d{\bf r}\int_{V}d{\bf r}^{\prime}{h_{\alpha\beta}(|{\bf r}-{\bf r}^{\prime}|)}{v(|{\bf r}-{\bf r}^{\prime}|)}$$

(14)

with $v(r)=1/r$. For $V\rightarrow\infty$, this simplifies to

$$H_{\alpha\beta}^{\infty}=\int_{0}^{\infty}h_{\alpha\beta}(r)v(r)4\pi r^{2}dr\;.$$

(15)

While we focus on the bare Coulomb potential $v(r)=1/r$, the present theory can in principle be applied to any potential form $v(r)$. For the special case $v=1$, $H_{\alpha\beta}^{\infty}$ are the Kirkwood-Buff integrals (KBI, usually denoted $G$ rather than $H$) of solution theory [10]. Previously we have extended KBI theory to finite volumes [11]. We showed that the bad convergence of the usual, truncated integrals can be avoided by an exact transformation from double volume integrals to one-dimensional radial integrals. Using this formalism [12, 13], we rewrite Eq. (14) as

$$H_{\alpha\beta}^{V}=\int_{0}^{\infty}h_{\alpha\beta}(r)v(r)w^{V}(r)dr$$

(16)

where

$$w^{V}(r)=\frac{1}{V}\int_{V}d{\bf r}_{1}\int_{V}d{\bf r}_{2}\delta(r-|{\bf r}_{1}-{\bf r}_{2}|)$$

(17)

is a weight function. Note that Eq. (16) is actually a finite integral, since $w^{V}(r)=0$ for $r>L_{\rm max}$ where $L_{\rm max}$ is the maximum distance in $V$. Equations (6,13,16,17) provide an exact expression for the ensemble averaged Madelung potential in finite volumes. The weight function is conveniently expressed as

$$w^{V}(r)=4\pi r^{2}y(r/L)$$

(18)

where $L=6V/A$ and $A$ is the surface area. The function $y(x)$ only depends on the shape of the volume [13]. For a sphere of diameter $L$, the exact expression is [11]

$$y_{s}(x)=(1-3x/2+x^{3}/2)\theta(1-x)$$

(19)

where $\theta(x)$ is the unit step function ($\theta$=0 for $x$$<$0 and $\theta$=1 for $x$$>$0). Analytic expressions of $y(r)$ are also known for cube and cuboid [13]. For any other shape, the function can be easily computed numerically [14]. In order to accelerate convergence to the thermodynamic limit, several extrapolations from the finite volume integrals have been proposed [11, 13, 15]. Here we focus on the second order expression of Ref. [13], with weight function $u_{2}(r)=4\pi r^{2}y_{2}(r/L)$, where

$$y_{2}(x)=(1-23x^{3}/8+3x^{4}/4+9x^{5}/8)\theta(1-x)$$

(20)

For further reference, we also define

$$y_{0}(x)=\theta(1-x)$$

(21)

which corresponds to simple truncation of radial integrals or lattice sums at $r=L$.

In the following, we consider a crystal with unit cell of volume $V_{c}$ containing $m$ point charges $q_{\alpha}$, $\alpha=1\dots m$, at positions ${\bf r}^{\alpha}$. We take each atom in the unit cell as a different species $\alpha$, such that $\rho_{\alpha}=1/V_{c}$ for all $\alpha$. The spherically averaged pair distribution function, appropriate for powder samples, is given by

$$g_{\alpha\beta}(r)=\frac{V_{c}}{4\pi r^{2}}\sum^{\prime}_{\bf T}\delta(r-|{\bf r}^{\alpha}-{\bf r}^{\beta}+{\bf T}|)\;,$$

(22)

where ${\bf T}$ are lattice vectors and the primed sum indicates that the term ${\bf T}=0$ is excluded for $\alpha=\beta$. Equivalently, Eq. (22) may be written as a sum over shells,

$$g_{\alpha\beta}(r)=\frac{V_{c}}{4\pi r^{2}}\sum_{k}n_{k}^{\alpha\beta}\delta(r-R_{k}^{\alpha\beta})$$

(23)

where $n_{k}^{\alpha\beta}$ is the number of $\beta$ ions on shell number $k$ of radius $R_{k}^{\alpha\beta}=|{\bf r}^{\alpha}-{\bf r}^{\beta}+{\bf T}|>0$ around an $\alpha$ ion. Inserting Eq. (23) into Eq. (16) we obtain

$$H_{\alpha\beta}^{V}={V_{c}}\sum_{k}n_{k}^{\alpha\beta}v(R_{k}^{\alpha\beta})y(R_{k}^{\alpha\beta}/L)-B^{V}$$

(24)

where

$$B^{V}=\int_{0}^{\infty}v(r)y(r/L)4\pi r^{2}dr\;.$$

(25)

The term $B^{V}$ comes from the constant $1$ which was subtracted when going from $g_{\alpha\beta}$ to $h_{\alpha\beta}$ in Eqs (11,13). By virtue of Eq. (12) the $B^{V}$ terms cancel upon summing over $\beta$ in Eq. (13). While $B^{V}$ has no direct physical meaning, it is important for the convergence of the individual terms $H_{\alpha\beta}$ in Eq. (24). For a sphere of diameter $L$, we obtain $B^{V}=4\pi L^{3}\int_{0}^{1}v(xL)y_{s}(x)x^{2}dx$ with $y_{s}$ given in Eq. (19). For the Coulomb interaction, $v(r)=1/r$, this yields $B^{V}=2\pi L^{2}/5=(12/5)V/L$. For $v=1$ (KBI) we have $B^{V}=\pi L^{3}/6=V$.

In the following we consider a sphere of diameter $L$ and calculate, with $y=y_{s}$ (Eq. 19), the exact finite volume integrals $H^{V}_{\alpha\beta}$ (Eq. 24) and the corresponding ensemble averaged Madelung constants (Eqs 6,13). We also compute extrapolations of $H_{\alpha\beta}^{V}$ to infinite volume, obtained with $y=y_{2}$ (Eq. 20) as well as the usual, unweigthed sums truncated at $r=L$, obtained with $y=y_{0}$ (Eq. 21). The convergence of the Madelung constant of the NaCl structure is shown in Fig. 2. The exact, infinite lattice value is $M^{\infty}=1.7475646$. The relative error $|M(L)/M^{\infty}-1|$ is plotted as a function of system size $L$. The exact finite volume result ($y_{s}$) converges very smoothly as $1/L$ while the $y_{2}$-extrapolation ($y_{2}$) converges as $1/L^{2}$. The truncated sum with $y=y_{0}$ corresponds to the most simple summation order over spherical shells. As it is well known, this summation order strongly diverges for the NaCl structure [4] and the result is not shown in Fig. 2.

We have examined the effect of charge neutralization using the shifted potential method by Wolf et al. [7], which consists in replacing the true Coulomb interaction $q_{i}q_{j}/r_{ij}$ by $q_{i}q_{j}(1/r_{ij}-1/L)$. For the usual, truncated sum ($y_{0}$), this is equivalent to adding a term $-Q/L$ to the potential, where $Q$ is the total charge in the sphere of radius $L$ [7, 8]. Charge neutralization makes the truncated sum converge, roughly as as $1/L$, as seen from the curve $y_{0}$-N in Fig. 2. Interestingly, for the exact finite volume integrals ($y_{s}$), charge neutralization has almost no effect. Indeed, the curve with charge neutralization ($y_{s}$-N) is almost identical to that without ($y_{s}$). This indicates that the positional averaging (see Fig. 1), which is implicit in $y_{s}$, largely cancels charge imbalances and makes explicit charge neutralization unnecessary. For the $y_{2}$ extrapolation, however, charge neutralization speeds up convergence considerably and the charge-neutralized sums ($y_{2}$-N) converge as $1/L^{3}$, in the same way as $y_{2}$-extrapolated, proper KBI integrals [13]. For comparison, the Wolf method [7] is also shown with an Ewald damping parameter $\alpha=0.5$. Because of the ${\rm erfc}(\alpha r)$-like damping term used in this scheme, the lattice sum converges exponentially fast and thus outperforms any method with power-law convergence, such as the present one. However, the Wolf method requires an empirical parameter ($\alpha$), whose optimum value depends on the problem at hand. As seen from Fig. 2, with the present scheme $y_{2}$-N, a practically reasonable error of $10^{-3}$ is achieved with a very moderate cut-off $L=10d$, which is comparable to $L=6d$ needed with the Wolf method. In conclusion of this section we have shown that ensemble averaged Madelung constants, Eqs (6,11), converge unconditionally even without charge compensation or damping factors. The $y_{2}$-extrapolation with charge neutralization converges as $1/L^{3}$, which is quite fast, albeit slower than the Wolf method.

In summary, we have developed a theory of ensemble averaged Madelung energies in finite volumes. The Madelung constants approach the thermodynamic limit in a universal way, independent of the structure, may it be crystalline or amorphous. A simple, general expression has been derived for the Madelung part of the termination averaged surface energy. The size dependence of the Madelung energies helps to understand the relative stability of different ionic structures as a function of system size, and provides a simple explanation for the phase transition from the rocksalt to the cesium chloride structure, which is observed in CsX clusters (X=Cl,Br,I).