Programming Tips

Secrets of the Madelung Constant

by Joe Buhler and Stan Wagon

A constant that describes the energy inherent in salt has connections with subtleties of analytic functions and number theory.

The Madelung constant is truly an interdisciplinary number. It is defined using ideas at the crossroads of physics and chemistry, but it is really a purely mathematical object, one that calls upon some quite subtle mathematical ideas for its evaluation. While this column will not present any exciting programming tips, it provides many illustrations of one of the fundamental tenets of sophisticated programming: to dramatically speed up a program, change the algorithm, not the implementation.

Consider a rock salt crystal with ions of sodium (positive) and chloride (negative) arranged alternately in the familiar square lattice structure. The Madelung constant measures the amount of electrostatic energy that holds a single ion in place. For simplicity, we imagine that the charge at each ion is ±1 and the distance between neighboring ions is 1 (in physical units, the Madelung constant has to be multiplied by q/R2 where q is the charge at a single ion and R the distance between adjacent lattice points). We imagine the lattice as being infinite in all directions with a sodium ion at the origin. Then, under our assumptions, the electrical force exerted on that ion by one at distance r is ±1/r2. The energy required to move the first ion all the way to infinity is obtained by integrating the force from r to infinity, the result being ±1/r.

So the Madelung constant is the sum
                   
                                   (-1)x+y+z
         ∑      -----------------------
     x,y,z not all 0          √x2 + y2 + z2

It is useful to be more general and define the constant in any dimension, and also with a variable complex exponent in the denominator instead of 1/2.

                                                  (-1)x+y+z
       Md(s)   =    ∑     -----------------------
                    x,y,z not all 0         (x2 + y2 + z2)s

Now, the problem, at least in the physically interesting cases of d = 2 or 3 and s = 1/2, is that this series is conditionally convergent. The sum can be rearranged to converge to any desired target! (For some reminders and illustrations regarding rearrangements of conditionally convergent series, see [Packel and Wagon 1994].)

The conditional convergence can be viewed as a reflection of the fact that the slow decrease of the Coulomb potential 1/r causes the shape of an actual finite lattice and its surface to be important. See solid-state physics texts for more details: [Kittel 1986] contains a discussion of experimental measurements; [Ashcroft and Mermin 1976] and [Slater 1967] give physical interpretations of the constant and of the conditional convergence of the sum defining it.

Warmup: Two Dimensions

It is worthwhile to begin with some 2-dimensional Madelung computations. Let's start with  M2(1/2), the alternating sum of 1/√x2 + y2, and take a naive brute force approach using expanding squares; symmetry allows us to consider only one side and multiply by 4. The following code examines the square of side-length 201, and displays that result and the nine preceding partial sums by using FoldList and Take.

M2SquareShell[k_] :=
  4 Sum[(-1)^(k+y)/N[Sqrt[k^2 + y^2]], {y, -k, k - 1}]
 
terms = Map[M2SquareShell, Range[100]];

Take[FoldList[Plus, 0, terms], -10]

Convergence to something near -1.6 seems plausible. Now consider expanding circles. In this case there is an important simplification: easy reasoning shows that the lattice points on each circle have the same charge, and so a circle of radius √n, where n denotes x2 + y2, contributes (-1)n r(d, n)/√n to the sum, where r(d, n) is the number of ways of writing n as a sum of d squares. Wagon's previous column in this journal [Wagon 1996] contains a method and some Mathematica code for computing r(d, n). In fact, this function will be included in version 3.0 as SumOfSquaresR in the NumberTheoryFunctions package. The piece that we need here is small, so we give the necessary code in full without comment.

Attributes[r] = Listable;
r[_, _?Negative] := 0
r[_, 0] := 1
r[0, _] := 0
r[2, n_Integer?Positive] := r[2, n] = 
   4 Apply[Plus, Im[I ^ Divisors[n]]]
r[d_, n_Integer?Positive] := r[d, n] = 
      r[d-1, n] + 2 Apply[Plus, r[d-1, n - Range[Sqrt[N[n]]]^2]]

We now define the circular shell's contribution and look at the partial sums to radius 100. Because of the circular setup, we have to examine each integer between 1 and 1002 to see how many ways it is a sum of two squares.

M2Circle[k_] := (-1)^k r[2, k] / Sqrt[N[k]]
terms = Map[M2Circle, Range[10000]];
Take[FoldList[Plus, 0, terms], -10]

This computation is inconclusive. Of course, the pervasive problem in this whole area is the slow convergence of the series involved; the physical interpretation of this slow-convergence is the long range effect that a given ion can have on the energy. We can go a little farther and examine the results graphically.

nMax = 20000;
terms = Map[M2Circle, Range[nMax]];
partialSums = FoldList[Plus, 0, terms];
ListPlot[partialSums[[Range[1, nMax, 20]]], 
  PlotJoined -> True,
  PlotStyle -> AbsoluteThickness[0.5],
  Ticks -> {Table[{i, 20 i}, {i, 200, 1000, 200}], Automatic}];

FIGURE 1:

This view also supports the conjecture of convergence to roughly -1.6. This brings us to our first curiosity: It seems as if expanding squares and expanding disks might yield the same sum. Before getting to the bottom of this, let's consider the related problem of M2(1). The brute force approach seems to converge.

M2Circle[k_] := (-1)^k r[2, k] / N[k]
terms = Map[M2Circle, Range[1000]];
Take[FoldList[Plus, 0, terms], -10]

A seasoned investigator always divides by π.

% / N[Pi]

Curious: the resulting number seems close to -log 2. The reason this is weird is that the average value of r(2, n) is π; by this we mean that the limit of (1/n) ∑i≤n r(2, i) is π. The proof is not difficult and depends only on the fact that the area of a circle is πr2 (see [Honsberger, essay 8]). Here's a quick computation in support.

Sum[r[2, n] / 2000., {n, 1, 2000}]

Assuming this result, the apparent convergence to -π log 2 is both bizarre and attractive. For if we simply replace r(2, n) by π, its average value, then the series reduces to ∑(-1)n π/n, which is an alternating harmonic series with sum -π log 2, our tentative estimate of the actual circular sum for M2(1). Of course, such a replacement is by no means legal in general. Yet it has been proved that it works here: the circular partial sums for M2(1) do converge to -π log 2 ([Borwein and Borwein 1987, chap 9], [Borwein, Borwein and Pinner forthcoming]).

Now we refocus on M2(1/2) and try to make sense out of the -1.6 we saw earlier. In this case, replacing r(2, n) by π is not at all fruitful. But in fact there is something of a closed form. Analytic methods [BB] yield that, for re(s) > 0, M2(s) = –4 a(s) b(s) where 

                   ∞     (-1)n+1                              ∞       (-1)n     
      a(s) =  ∑   ---------      and     b(s) =  ∑    -----------  
                  n=1        ns                                               n=0      (2n + 1)s

Specializing to s = 1/2, it seems quite miraculous that an alternating sum of r(2, n)/√n should be the product of two very simple alternating sums of square roots. The functions a(s) and b(s) can be expressed in terms of the (generalized) Riemann zeta function. More precisely, it is well known that a(s) = (1 - 21-s) z(s) (although SymbolicSum gives a slightly different form), and b(s) can be found by using SymbolicSum as follows.

Needs["Algebra`SymbolicSum`"]  (* unnecessary in version 3.0 *)
Sum[(-1)^n / (2n + 1)^s, {n, 1, Infinity}]

We remind the reader that the Riemann zeta function Zeta[s] is defined to be ∑n=1∞  n-s and that the generalized (or Hurwitz) zeta function Zeta[s, a] is ∑n=0∞  (n + a)-s, for 0 < a ≤ 1. Although these series only converge for s > 1, they have analytic continuations to (almost all of) the complex plane. This means, for instance, that there is a function z(s) defined for all complex numbers (except s = 1) such that z(s) is equal to the sum ∑n=1∞ n-s for any complex number with real part bigger than 1. This notion of analytic continuation plays an important role in analytic number theory [Apostol 1976] and will play a role in the Madelung problem, as described below.

The summations a(s) and b(s) actually converge at s = 1 and it is well known that a(1) = log 2, and b(1) = π/4, which agrees with the earlier -π log 2 result. The ab formula can now be used for a rapid numerical computation.

alpha[s_] := (1 - 2 ^ (1-s)) Zeta[s]
beta[s_]  := (Zeta[s, 1/4] - Zeta[s, 3/4])/4^s

-4 alpha[1/2] beta[1/2]

% //N

The fact that this is within 1% of the golden ratio is merely a coincidence: the law of small numbers (i.e., there aren't enough of them) strikes again.

We emphasize that the results of this section are all known, and the computations do not shed any light on the proofs. However, note that changing the algorithm caused an enormous speedup, illustrating our fundamental precept: to dramatically speed up a program, change the algorithm!

Three Dimensions: The Constant for Salt

For the three-dimensional case we again begin with the expanding cubes method and a brute force sum. We first sum the points on the interior of a face and multiply by 6, then sum along the interior of an edge and multiply by 12, and then add in the 8 vertices.

M3Cubes[k_] :=
  6 Sum[(-1)^(k + y + z) / Sqrt[N[k^2 + y^2 + z^2]], 
              {y, -k + 1, k - 1}, {z, -k + 1, k - 1}] +
 12 Sum[(-1)^(2 k + y) / Sqrt[N[2 k^2 + y^2]], {y, -k + 1, k - 1}] +
  8 (-1)^(3 k) / Sqrt[N[3 k^2]]

Now we build up the sum for a cube of side-length 101 by computing the contribution of each shell and looking at the last 10 partial sums.

terms = Map[M3Cubes, Range[50]];

Take[partialSums = FoldList[Plus, 0, terms], -10]

ListPlot[partialSums, PlotJoined -> True];

FIGURE 2:

Convergence is, as usual, slow, but there is some indication of oscillatory convergence to about -1.75. In fact, this slow convergence is what makes the Madelung constant so interesting: it is another way of saying that an accurate computation of the electrostatic energy holding an ion in place requires, because of the small denominators (contrast with the more familiar inverse square scenario), that a very large number of ions be taken into account.

Now let us switch to expanding spheres, which perhaps seems more natural. In this case the computation is simpler, because the parity on a given spherical shell does not change. If we use n to denote x2 + y2 + z2, then the contribution for the shell of radius √n is just (-1)n r(3, n)/√n.

M3Spheres[k_] := (-1)^k r[3, k]/Sqrt[N[k]]

terms = Map[M3Spheres, Range[100]];

ListPlot[FoldList[Plus, 0, terms], PlotJoined -> True];

FIGURE 3:

This doesn't look good at all. The partial sums are oscillating wildly, and in fact it can be proved that they diverge [Borwein, Borwein and Pinner forthcoming, Buhler and Crandall 1990].

If the sum is so temperamental, how should it really be defined? As discovered by P. P. Ewald around 1920, there are several ways to specify the constant, and they are all equivalent to saying that the value can be obtained by analytic continuation (see [Slater 1967]). The sum M3(s) can be proved to converge for s > 3/2, and it can also be proved that there is an analytic continuation of this function to the entire complex plane. The Madelung constant is the value of this function at s = 1/2. This idea can be derived by solving Poisson's equation for the potential energy, or by replacing point charges by Gaussian distributions and letting the width of the distribution go to zero. The general idea, namely that the values of nonconvergent sums should be obtained by analytic continuation from regions where the sums do converge, is used elsewhere in physics.

Advanced Formulas

It turns out that there are many beautiful formulas, some fairly new, that can help us evaluate the Madelung constant M3(1/2) with great accuracy. Here are four such.Throughout, we use R to represent √i2 + j2. For reference, here is the constant to 50 digits, obtained by the last of the four formulas that follow: M3(1/2) = -1.74756459463318219063621203554439740348516143662474... .

Benson, 1956 [Borwein and Borwein 1987]:    M3(1/2)  =  -12 π ∑ sech2(πR/2), where the sum is over	
                                        															                   positive odd integers i, j.

Hautot, 1975 [Crandall 1994]:        M3(1/2)  =  -π/2 + 3 ∑ (-1)i  csch(π R)/R, where the sum is over
                                                        positive and negative integers i and j, excluding (0, 0).
 
Crandall–Buhler, 1987 [Crandall and Buhler 1987]   
                                     M3(1/2)  =  -π + √2 + ∑k=1∞ (-1)k r(3, k) (1 - tanh(π √k))/√k

Crandall–Buhler, 1987 [Crandall and Buhler 1987]     
                                      M3(1/2)  =  -2π + 27/4(∑k=1,∞  exp(-π/√2 (k - 1/2)2))2 + 
                                       2 ∑k=1∞  (-1)k r(3, k)/(√k (1 + exp(4 π √k)) )
                                                                      
These are all easy to implement, though there are tricks to speed things up. For Benson's formula, symmetry allows us to assume i < j, provided we multiply by 2 and add back the main diagonal. The definition of R using a delayed assignment makes it work inside the various sums. The 19th partial sum (i.e., allowing i and j to range up to 19), gives 26 correct digits. These formulas have similar general shapes, which guarantee that the truncation error isn't too much larger than the first omitted term; thus if the sum has stopped changing in the 27th digit, it is probably accurate to 27 digits; more rigorous estimates on the error will be given later in this section.

R := Sqrt[i^2 + j^2]

Benson[n_] := -12 Pi *
   (2 Sum[Sech[Pi R/2]^2, {i, 1, n, 2}, {j, 1, i-1, 2}] + 
      Sum[Sech[Pi i/Sqrt[2]]^2, {i, 1, n, 2}])

N[Benson[19], 27]

Note that, by agglomerating the lattice points at the same distance from the origin, we can rewrite Benson's formula as

       -12π + ∑k=2,even∞  sech(π √k/ 2)2 rodd(2, k) / 4

where rodd(2, n) is the number of representations of k as a sum of 2 odd squares. It is easy to compute rodd: it equals 0 unless n ∫ 2 (mod 4), in which case it equals r(2, n). But for small computations this is slower, because, to get the same region as computed by Benson[n], one must use n2 in the rodd formula. Perhaps for large computations this approach would be worthwhile. These comments apply to the Hautot formula as well.

The Hautot formula can be split into pieces according as i and j are both even or both odd (if they have opposite parity the contribution is 0 because the (-1)i term leads to full cancellation). Using the same ideas as for Benson's formula, with the additional complication of the points on the axes, leads to the following implementation. As with Benson, n = 19 gives 26 correct digits; but Benson is a little faster.

HautotOdd[n_] := 
   (2 Sum[Csch[Pi R]/R, {i, 1, n, 2}, {j, 1, i-1, 2}] + 
      Sum[Csch[Pi (Sqrt[2] i)]/(Sqrt[2] i), {i, 1, n, 2}])
      
HautotEven[n_] := 
   (2 Sum[Csch[Pi R]/R, {i, 2, n, 2}, {j, 2, i-1, 2}] + 
      Sum[Csch[Pi (Sqrt[2] i)]/(Sqrt[2] i), {i, 2, n, 2}])
      
HautotAxes[n_] := Sum[Csch[Pi i]/i, {i, 2, n, 2}]

Hautot[n_] := -Pi/2 + 12 * 
  (HautotAxes[n] + HautotEven[n] - HautotOdd[n])

N[Hautot[19], 27]

The Crandall-Buhler tanh formula is cute because simply ignoring the sum reduces it to √2 - π, which is in fact correct to 1%! Here 86 summands are required for 26 digits, and it is a little slower than Benson and Hautot.

CBTanh[n_] := -Pi + Sqrt[2] + 
  Sum[(-1)^k r[3, k]/Sqrt[k] (1 - Tanh[Pi Sqrt[k]]), {k, n}]
  
N[CBTanh[86], 27]

To implement the Crandall-Buhler exponential formula, we divide it into a main part and a correction part. This is the fastest of the formulas: 3 terms of the main sum and one correction term yield 6 digits. And though it is essentially tied with Benson's formula in terms of time to get 26 digits, for larger computations (100 or more digits) it is by far the fastest of the formulas (because of the extra 4 in the exponent in the correction part).

CBMain[n_] := -2 Pi +  
  2^(7/4) Sum[ Exp[-Pi / Sqrt[2] (k - 1/2)^2], {k, -n, n}] ^ 2

CBCorrection[m_] := 
  2 Sum[ (-1)^k r[3, k] / Sqrt[k] / (1 + Exp[4 Pi Sqrt[k]]),
    {k, 1, m}]

CB[n_, m_] := CBMain[n] + CBCorrection[m]

CB[3, 1] //N

N[CB[5, 21], 27]

All of the sums in the four formulas above can be written as one-dimensional sums whose kth terms are bounded in absolute value, for large k, by a kb exp(-c kd), where a, b, c, d are constants. The most important constant is d; it is 1/2 for all of the formulas (except for the CBMain formula, for which it is 1). The constant c is equal to π/2 for the Benson formula, π for the Hautot formula, 2π for the first Crandall-Buhler formula, and 4π for the correction term in the second Crandall-Buhler formula. Although the precise rate of convergence of the sum is determined by many factors (such as the constants a and b, whether or not the sum is alternating, the coarseness of the estimates leading to the a kb exp(-c kd) bound), in our summations the constant c has the biggest effect. It gives us a pretty good idea of the relative rates of convergence, which is borne out by our empirical investigations.

The error in stopping the sum at the nth term is equal to the sum of the kth terms for k > n. These sums can usually be bounded by using coarse estimates, such as taking absolute values and bounding r[3, k] by k3, and then using the idea of the integral test to bound the resulting sum. The details are similar in spirit in all cases; we illustrate by applying this idea to the two terms in the second Crandall-Buhler formula.

Suppose that we want 100 digits. The CBMain term is bounded by an integral in a fairly obvious way (an extra factor of about 4 is needed at the front because the sum is squared, but this is negligible for the case at hand).

TailBound1[n_] := 2^(7/4) Integrate[
  Exp[-Pi / Sqrt[2] (k - 1/2)^2], {k, n, Infinity}]
  
N[TailBound1[11], 110]

So 11 terms of the main formula are adequate for 100-digit precision. For the correction sum we take the absolute value, bound r[3, k] by k3, and ignore the "1 +" in the denominator. It is convenient to substitute u = √k manually.

TailBound2[n_] :=
  4 Integrate[u^6 Exp[-4 Pi u], {u, Sqrt[n], Infinity}]
  
N[Log[10, TailBound2[400]]]

So we can, with confidence, get 100 Madelung digits as follows.

Madelung = N[CB[11, 400], 100]

The power of these formulas obtained by analytic methods becomes evident upon comparison with the expanding cubes method. While it can be proved that that method does give the proper result, a sum over a million lattice points (the computation on page p) returned a result accurate to only two digits. In fact, as we just saw, for 6-digit accuracy the Crandall-Buhler formula will do the job instantaneously via CB[3, 1].

As described in [Slater 1967], the physical model giving rise to the Madelung sum is an idealization (e.g., various quantum effects are ignored and the physical dimensions of point charges are taken to be negligible). The value of the Madelung constant is probably physically relevant only to a few decimal places at most. Nonetheless, the mathematical expressions are interesting in and of themselves and are relevant to further understanding of the underlying physics.

The Delord Interpretation

Physicist Jean Delord (Reed College) suggested that the problem with the spherical summation approach for M3(1/2) was that the partial sums do not correspond to an electrically neutral lattice. Thus he conjectured that a convergent summation method could be obtained by adding back to the surface of the sphere an amount of charge that neutralizes the total charge of the sphere. Since the total spherical charge, assuming radius R, is the sum of (-1)x+y+z over all points inside or on the sphere of radius R, this requires adding ∑k≤R2  (-1)k r(3, k). Letting n denote R2 as usual, this leads to an energy-correction term of (1/√n) ∑k≤n  (-1)k r(3, k). Let us investigate this idea numerically. First we compute some spherical partial sums, as on page q; we let n run up to 2000, which means we are examining spheres up to radius 44.7.

nMax = 2000;
terms = Map[M3Spheres, Range[nMax]];
partialSums = Rest[FoldList[Plus, 0, terms]];

Now we define the missing charge by first accumulating the electrical corrections into a list, and then dividing by the distances (the square roots) to get the energy corrections.

chargeCorrex = 
  FoldList[Plus, -r[3, 1], Table[(-1)^n r[3, n], {n, 2, nMax}]];

energyCorrex = chargeCorrex /  Sqrt[Range[1., nMax]];

We subtract the corrections and note that, indeed, numbers in the general vicinity of the Madelung constant do show up.

Take[partialSums - energyCorrex, -10]

Let's look at the entire sequence of corrected spherical partial sums.

ListPlot[partialSums - energyCorrex, PlotJoined -> True,
  Epilog -> Line[{{1, Madelung}, {nMax, Madelung}}],
  PlotRange -> {-1.85, -1.4}, AxesOrigin -> {1, -1.85},
  Ticks -> {Range[500, nMax, 500], Range[-1.8, -1.4, 0.1]},
  PlotStyle->{AbsoluteThickness[0.5]}];

FIGURE 4:

As usual, the slow convergence leads to a somewhat inconclusive picture, but it is clear that the electrical neutralization has decreased the wild oscillations that we saw when trying to sum by expanding spheres. In fact, Delord's conjecture has been proved by Buhler and Crandall [Buhler and Crandall 1987]. As an exercise, the reader might enjoy investigating the effect of charge neutralization in the 2-dimensional case when summing by expanding circles. It seems to be an interesting open problem to determine what correction terms need to be introduced for d > 3 so that Md(1/2) can be summed by expanding balls in d dimensions. Finally, we mention that [Borwein, Borwein and Pinner forthcoming] has some results, and some open problems, on the summability of Md(s) when the summation is done over expanding regions of various shapes other than the cubes and spheres considered here.

REFERENCES

APOSTOL, TOM.  Introduction to Analytic Number Theory. Springer, New York, 1976.
ASHCROFT, NEIL W.  and MERMIN, N. DAVID. Solid State Physics. Saunders, Philadelphia, 1976.
BORWEIN, JONATHAN  and BORWEIN, PETER.Pi and the AGM — A Study in Analytic Number Theory and  Computational Complexity. Wiley, New York, 1987.
BORWEIN, DAVID, BORWEIN, JONATHAN and PINNER, CHRISTOPHER. Convergence of Madelung-like lattice sums, forthcoming.
BUHLER, J. P.  and CRANDALL, RICHARD E.  On the convergence problem for lattice sums, Journal of Physics, Series A: Mathematical and General  23 (1990) 2523-2528.
CRANDALL, RICHARD E. Projects in Scientific Computation. Springer/TELOS, New York, 1994.
CRANDALL, RICHARD E., and BUHLER, J. P.  Elementary function expansions for Madelung constants, Journal  of Physics, Series A: Mathematical and General 20 (1987) 5497-5510.
HONSBERGER, R.Ingenuity in Mathematics. New Mathematical Library, vol. 23, Random House, New York, 1970.
KITTEL, CHARLES.Introduction to Solid State Physics, sixth ed. Wiley, New York, 1986.
PACKEL, ED, and WAGON, STAN. Rearrangement patterns for the alternating harmonic series, Mathematica in Education 3:2 (Spring, 1994) 5-10.
SLATER, J.C.  Insulators, Semiconductors, and Metals, vol. 3 of Quantum Theory of Molecules and Solids, McGraw-Hill, New York, 1967.
WAGON, STAN. The magic of imaginary factoring. Mathematica in Education and Research 5:1 (Winter, 1996) 43-47.

Author information:

Stan Wagon uses Mathematica extensively in his research, teaching, and exposition and is especially appreciative of the new ways in which Mathematica allows us to look at mathematical objects. He is the author of several books, including Mathematica in Action and, with Ed Packel, Animating Calculus. He also teaches a high-altitude Mathematica course every summer in the mountains of Colorado.

Stan Wagon
Department of Mathematics
Macalester College
St. Paul, MN 55105
wagon@macalstr.edu

Joe Buhler is interested in many aspects of computational number theory and algebra, and has used Mathematica for many years as part of that research. He teaches at Reed College, and when not doing mathematics has been known to juggle and play Go.

Joe Buhler
Department of Mathematics
Reed College
Portland, OR 97202
jpb@reed.edu

ELECTRONIC SUBSCRIPTIONS

Included in the distribution for each electronic subscription is the file madelung.ma containing the text and code of the material described in this column.