Periodic systems IV. Conditional convergence

Four years after the first three notes in this notebook, I am attempting to understand the electrostatics I was implicitly leaning on back then. This note is the explanation I wish I had at the time.

The problem

Let’s write down the Coulomb energy of a crystal and see where it goes wrong. Take a neutral unit cell and tile it over all space. The naive Coulomb energy per cell is a lattice sum of $1/r$ interactions between everything in the central cell and everything in every image cell. For point charges $q_i$ at positions $\mathbf{r}_i$ it reads

$$V = \frac{1}{2}\sum_{\mathbf{T}}{\vphantom{\sum}}' \sum_{ij} \frac{q_i q_j}{|\mathbf{r}_i - \mathbf{r}_j - \mathbf{T}|}$$

where $\mathbf{T}$ runs over lattice translations and the prime excludes the self-interaction at $\mathbf{T} = \mathbf{0}$. Notice that each class of terms (nuclear repulsion, electron-nuclear attraction, electron repulsion) diverges on its own. The number of images at distance $r$ grows like $r^2$ while the interaction only falls like $1/r$, and only the charge-neutral combination therefore has a chance.

Is the neutral combination convergent then? Conditionally, yes, and the qualifier matters. Recall from undergraduate analysis (the Riemann series theorem) that a conditionally convergent series can be rearranged to sum to any value you like. Here the rearrangement has a physical meaning. Summing the images in growing spheres, growing cubes, or slabs amounts to assembling a finite crystal of that shape and letting it grow, and the electrostatic energy of a finite body genuinely depends on its shape, because a dipole layer on the boundary produces a field in the interior no matter how large the body grows. The “ambiguity” of the lattice sum is real physics of surfaces that the infinite-lattice idealization has to decide about, one way or another.

The standard decision is to separate the shape-independent bulk energy from the shape-dependent surface term and keep only the former. For a spherical crystal in vacuum the surface term is $(2\pi/3\Omega)\,|\mathbf{P}|^2$ with $\mathbf{P}$ the cell dipole moment, and embedding the crystal in a perfect conductor (the “tinfoil” boundary condition) kills it entirely. Once that choice is made, the remaining sum is well defined.

Ewald Summation

It is a common misconception that Ewald summation fixes the divergence. The conditional convergence is settled by the boundary-condition choice above. What remains is a well-defined sum that converges too slowly to compute directly, and Ewald’s trick makes it fast. Recall the identity

$$\frac{1}{r} = \frac{\operatorname{erfc}(\eta r)}{r} + \frac{\operatorname{erf}(\eta r)}{r}$$

which you can think of as surrounding each point charge with a canceling Gaussian charge cloud and adding the cloud back separately. Notice what each piece buys us. The erfc piece dies off like a Gaussian, and its real-space lattice sum thus converges within a few shells of neighbors. The erf piece is smooth, and smooth functions have rapidly decaying Fourier components, such that its reciprocal-space sum picks up factors of $e^{-G^2/4\eta^2}/G^2$ and converges in a few $\mathbf{G}$ vectors. The splitting parameter $\eta$ just balances the two costs. This is the same divide-and-evaluate logic that range-separated density fitting applies to the DF metric in note I.

The Madelung constant

What happens to a charge’s interaction with its own images? Evaluating the Ewald potential $v_\mathrm{E}(\mathbf{r})$ at the charge’s own position, with the bare $1/r$ singularity removed, gives

$$v_\mathrm{M} = \lim_{r\to 0}\,\big[v_\mathrm{E}(\mathbf{r}) - 1/r\big]$$

the interaction of a unit charge with all of its periodic images and the neutralizing background. For ionic crystals the same construction packaged per ion pair is the classic Madelung constant. Rock salt has $\alpha \approx 1.7476$, i.e. each ion pair in the infinite crystal is bound 1.7476 times more strongly than an isolated Na⁺Cl⁻ pair at the same distance, and the partial sums of the naive shell-by-shell version oscillate without settling.

The same quantity sits behind the “Ewald correction” from note II. Recall that the exchange energy carries an integrable $1/G^2$ divergence at $\mathbf{G} = \mathbf{0}$, the reciprocal-space fingerprint of the same long-range Coulomb tail, and shifting the occupied orbitals by $v_\mathrm{M}$ removes the slowest-decaying finite-size error. The pieces of machinery from the first three notes, the range separation, the Madelung shift, and the finite-size corrections, are all consequences of the same conditionally convergent sum.