The Madelung Constant in $N$ Dimensions

We introduce two convergent series expansions (direct and recursive) in terms of Bessel functions and representations of sums $r_N(m)$ of squares for $N$-dimensional Madelung constants, $M_N(s)$, where $s$ is the exponent of the Madelung series (usually chosen as $s=1/2$). The functional behavior including analytical continuation, and the convergence of the Bessel function expansion is discussed in detail. Recursive definitions are used to evaluate $r_N(m)$. Values for $M_N(s)$ for $s=\tfrac{1}{2}, \tfrac{3}{2}, 3$ and 6 for dimension up to $N=20$ and for $M_N(1/2)$ up to $N=100$ are presented. Zucker's original analysis on $N$-dimensional Madelung constants for even dimensions up to $N=8$ and their possible continuation into higher dimensions is briefly analyzed.


I. INTRODUCTION
The classical lattice energy E lat of an ionic crystal M + X − can be obtained from lattice summations of Coulomb interacting point charges and is usually presented by the Born-Lande form [1,2] where M lat is the Madelung constant for a specific lattice [3], N A is Avogadro's constant, and n is the Born exponent which corrects for the repulsion energy V = aR −n , a > 0 at nearest neighbor distance R 0 , Z is the ionic charge (+1 in the ideal case), e and 0 are the elementary charge and vacuum permittivity respectively. Values for Z 2 M and n have been tabulated for different crystals in the past [4]. For a simple cubic lattice with alternating charges in the crystal the Madelung constant (or function) M (s) ≡ M sc (s) is given by the 3D alternating lattice sum where the summation is over all integer values, the prime behind the sum indicates that i = j = k = 0 is omitted, s ∈ R, and s = 1 2 is chosen for a Coulomb-type of interaction. This sum is absolutely convergent for s > 3 2 , but only conditionally convergent for smaller s-values [5,6]. The problem with conditionally convergent series is that the Riemann Series Theorem states that one can converge to any desired value or even diverge by a suitable rearrangement of the terms in the series. This problem is well known for the Madelung constant (s = 1 2 ) and has been documented and analyzed in great detail by Borwein et al [6][7][8] and Crandall et al [9,10]. For example, one has to sum over expanding cubes and not spheres to arrive at the correct result of M ( 1 2 ) = −1.747 564 594 633 182 . . . [10]. It is currently not known if the Madelung constant can be expressed in terms of standard functions. The closest formula one can get is the one for s = 1 2 recently derived by Tyagi [11] following an approach by Crandall [10], which is correct to 10 significant figures if the sum is neglected (for more recent work and improvement of Tyagi's formula see Zucker [12]). Moreover, the sum converges relatively fast. Here r 3 (k) is the number of representations of k as a sum of three squares.
There are many expansions available leading to an accurate determination of the Madelung constant [10]. Perhaps the most prominent formulas are the ones by Benson-Mackenzie [13,14] M 1 2 = −12π i,j∈N sech 2 π 2 (2i − 1) 2 + (2j − 1) 2 (4) and by Hautot [15] (in modified form by Crandall [10]) The Madelung constant can easily be extended to a N dimensional series (N >0), and the prime after the sum denotes that the term corresponding to i 1 = i 2 = · · · = i N = 0 is omitted (in the shorter notation on the right 1 = (1, 1, . . . , 1) ). The sum is absolutely convergent for exponents s > N 2 . The Madelung series is as a special case of the more general Epstein zeta function [16].
Zucker has found analytical expressions in terms of standard functions for even dimensions up to N = 8 [17], Here η(s) is the Dirichlet eta function, β(s) the Dirichlet beta function, and ζ(s) the Riemann zeta function [17]. These standard functions are defined in the Appendix A together with their analytical continuations to the whole range of real (or complex) numbers, s ∈ R(C). By analogy with the three-dimensional case, an N -dimensional lattice can easily be constructed from its N linearly independent basis lattice vectors (or transformations of it). Higher dimensional lattices and their properties have been catalogued (up to certain dimensions) by Nebe and Sloane [18]. The simple cubic N -dimensional lattice can be drawn as an infinite graph with atoms (vertices) and edges connecting the nearest neighbor atoms (adjacent vertices). If we walk around the edges we alternate the charges (+/-sign or red/blue color of the vertices in the graph) in the ionic lattice corresponding to the alternating series for the Madelung constant. We can also derive the lattice from tiling the N -dimensional space with N -cubes by implying translational symmetry. Figure 1 shows the graphs for such N -cubes up to N = 5 together with the alternating color scheme. We notice that for dimensions N > 3 the graphs are not planar anymore. The number of nearest neighbor vertices for an N-dimensional cubic lattice is 2N and corresponds to the limit, For example, Crandall reports M 3 (50) = −5.999 999 999 999 989 341 . . . [10]. A general and relatively fast converging series expansion for the N -dimensional Madelung constant has been elusive for a very long time. For example, a recent suggestion was made by Mamode to use the Hadamard finite part of the integral representation of the underlying potential (e.g. a Coulomb potential) within the N -dimensional crystal [19], but computations are quite involved and results presented were only up to three dimensions. For the N -dimensional case one can explore expansions known for example for the Epstein zeta function [9,20,21] or similar techniques [22]. In this work, we introduce a general formula for the N -dimensional Madelung constant for a simple cubic crystal in terms of a fast convergent Bessel function expansion allowing for analytical continuation, which gives deep insight into the functional behavior of the N -dimensional Madelung constant. The derivation is given in the next section. The convergence of M N (s) with increasing dimension N is discussed in detail in the results section.

II. THEORY
In this section we derive two useful expansions for the N -dimensional Madelung constant. Consider M N +1 (s) and change the last summation index to k, and write Now separate the sum into the two cases k = 0 and k = 0 to get where By the gamma function integral in the form ( we have By using the modular transformation for the theta function [23], we get This can be rearranged further to give (20) where N 0 denotes the natural numbers including zero, and r odd N (m) is the number of representations of m as a sum of N odd squares. That is, r odd N (m) is the number of solutions of in integers. The integral in (20) can be evaluated in terms of Bessel functions by means of the formula to give On using this result back in (13) we obtain the recursion relation for the Madelung constant in terms of the dimension N , For fixed N , the term r odd N (8m + N ) can become very large for larger m and N values, but is more than compensated by the exponentially decreasing Bessel function, which we discuss in detail in the next section. The r odd N (m) values can be determined recursively which is described in the Appendix.
While the recursion relation (24) is useful if the Madelung constant of lower dimension is known, we seek for a second formula where the recursion relation has been resolved. Here, we proceed as above and separate the sum for M N +1 (s) into two cases according to whether i 1 = i 2 = · · · = i N = 0 or i 1 , i 2 , . . ., i N are not all zero. This gives where Applying the integral formula for the gamma function and then the modular transformation for the theta function we obtain where the last step follows by noting In terms of the modified Bessel function this becomes, by (22), On using this back in (26) we obtain For the case of N = 0 the sum of the right-hand side is zero (r 0 (m) = 0) and we have M 1 (s) = −2η(2s) in agreement with Zucker's formula (7). We can conveniently write the sum in the form, with Note that the coefficients c s (m) are independent of the dimension N . The sum in (32) converges fast because of the exponential asymptotic decay of the Bessel function. The more problematic part is the convergence with respect to the first sum (see eq.31) over m as we shall see.
As a special case we evaluate M N (1/2). Letting s → 1/2 in (30) gives a formula for the N +1 dimensional Madelung constant where r N (m) is the number of representations of m as a sum of N squares. The coefficient c 1/2 (m) becomes where the integral is obtained using the formula [24] K 0 (z) = and summing the resulting geometric series. For example, taking N = 2 gives On the other hand, using (24) and Zucker's equation (7) we get The coefficients c 1/2 (m) are listed in Table III Table III [10]. For larger exponents the series converges much faster, i.e. for M 3 (6) ( Table III) we need to sum only over 1 ≤ m ≤ 51 to reach convergence to 14 significant digits behind the decimal point. Note that we used backwards summation as small numbers add up. We also checked our values for the even dimensions up to N = 8 with the values obtained from the analytical function in (7) by Zucker [17], and they are in perfect agreement.
To discuss the convergence behavior of the series (31) we observe that the coefficients c 1/2 (m) are rapidly decreasing with increasing m. However, at the same time the r N (m) values increase also rapidly with increasing m (and increasing N ) shown in Figure 4. The asymptotic behavior of the Bessel functions is well known, i.e.they decrease exponentially with increasing m, K s (x) ∼ (π/2x) 1 2 e −x . On the other hand, the sum of squares representation increases polynomially for fixed N [25][26][27], e.g. we know from Ramanujan's work that r 2N (m) = O(m N ) (derived from eq.(14) in ref. [28]). This is also seen in the logarithmic behavior of log 10 r N (m) in Figure 4. This implies that the Madelung series expansion in terms of Bessel functions is converging, but very slowly for higher dimensions because of a very large dimensional prefactor. This can clearly seen from the m max values for M N (1/2) in Table III. For M N (s), s ≥ 1/2 we approximately have m max ≤ nint(1.16N 2 + 11.5N + 73), where nint represents the nearest integer function.
Perhaps more problematic is the appearance of large numbers due to the r N (m) values in the sum over m in eq.(31) where one reaches soon the limit with double precision arithmetic at large N values. This is clearly seen in Figure  5 for the case of dimension 16 and s = 1/2 which shows for the individual terms a strong oscillating behavior and poynomial increase up to rather large values around m = 14 followed by an exponential decay. For higher dimensions this maximum shifts to higher m values before the exponential decay sets in. However, if we add pairs of positive and negative terms in the oscillating series to obtain new coefficients b(2m) = a(2m) + a(2m − 1), we experience a far smoother and better convergence behavior.
By using the recursive formula (24) instead we obtain much fast convergence as we reach the exponential decay far earlier because of the argument 8m + N in the Bessel function, see Figure 6. Here we avoid such large values and the strong oscillating behavior as the sign change appears in the summation over k in .(31) rather than in (24). Hence, for accuracy reasons eq.(24) is preferred, and we used this equation instead for the values in Table III.
which is independent of the dimension N . This implies that all Madelung curves cross at these critical points. Moreover, from the Dirichlet eta function we know that η(2s) = 0 for s = −1, −2, . . .. This behavior is nicely seen in Figure 7. Comparing with Zucker's formulas we see that this is easily fulfilled for the specific dimensions given. Concerning the usual Madelung constant at s = 1/2 we see that they lie close to the crossing point at s = 0 which explains their rather slow decrease with increasing dimension N . Zucker was able to evaluate the Madelung series analytically for even dimensions up to N = 8 [17] based on previous work of Glasser [29,30]. He further conjectured that M 3 (s) may be expressed in terms of a yet unknown Dirichlet series (for a recent analysis of lattice sums arising from the Poisson equation see Ref. [31]). Of considerable help for future investigations will be the condition that M N (0) = −1 and M N (−n) = 0 for all n ∈ N. At these critical points we have the following properties where B n and E n are the Bernoulli and Euler numbers respectively [23]. For example, from Zucker's formulas (10) and (11)  functions can be obtained for higher even dimensions. For a detailed discussion we refer to Appendix B. In this sense, our expansions in terms of Bessel functions is perhaps the closest general form for a fast convergent series we can get for the N -dimensional Madelung constant.

IV. CONCLUSIONS
We presented fast convergent expressions for the Madelung constant in terms of Bessel function expansions which allow for an asymptotic exponential decay of the series. Even for higher dimensions the Madelung constants can be evaluated efficiently and accurately through the recursive expression or by using computer algebra to work with the generating functions. The number of representations of the N sum of squares can also be efficiently handled through recursive relations. The Madelung constants and their analytical continuations can be calculated easily by standard mathematical software packages to any precision. These numbers may be useful for future explorations of analytical formulas in higher dimensions. For s ≥ 1/2 a Fortran program with double precision accuracy is available from our CTCP website [32].

Appendix A: Special Functions
We give a brief overview over the special functions used in this work. More details can be found in the book by Andrews [23]. The Dirichlet (or L-) series (Riemann zeta, Dirichlet eta, and Dirichlet beta functions) are defined as Their analytic continuations to L-functions into the negative real part (or the whole complex plane) are given by [29] η(−s) = s(2 − 2 −s )π −s−1 sin( π 2 s)Γ(s)ζ(s + 1) .
Here, the gamma function is usually defined for real positive numbers as and when s = n ∈ N we have Γ(n) = (n−1)!. The gamma function on the whole real axis is then defined as the analytic continuation of this integral function to a meromorphic function by the simple recursion relation Γ(x) = Γ(x + 1)/x with 1/Γ(−n) = 0 for n ∈ N 0 [33]. The modified Bessel function of the second kind is defined as The higher-order Bessel functions can be successively reduced to lower order Bessel functions by and we use the symmetry K −ν (x) = K ν (x) for its analytical continuation.
In a similar fashion one obtains a recursive formula for the sum of odd squares, keeping in mind that r odd N (0) = 0 and we do not include this term in our summation. For completeness we mention that the sum of even squares is trivially related to the sum of squares by r even N (4m) = r N (m) and r even N (m) = 0 if m is not divisible by 4.

Appendix B: Why Zucker's analytical formulas do not continue into higher dimensions
Zucker's formulas (7)-(11) are equivalent to Jacobi's formulas for sums of 2, 4, 6 and 8 squares (e.g., see [34] pp. 177, 202, 238): Put q = e −u , multiply both sides by u s−1 and integrate, to obtain i1,i2,...,i8∈Z The integrals can be evaluated using eq.(A7) to give i1,i2,...,i8∈Z where the common factor Γ(s) has been cancelled from each side. In other words, we have obtained Thus we have obtained Zucker's formula (11) from the sum of squares formula (B4). The process is reversible, so (11) is equivalent to (B4). By similar calculations, each of Zucker's formulas (8)-(11) is equivalent to the respective formula in (B1)-(B4). By analogy with M 8 (s) in eq.(11), it is tempting to speculate that there might be expressions for M 10 (s) and M 12 (s) as finite sums of the forms for certain L-functions f i (s), g i (s), F i (s) and G i (s). However this is unlikely to be true for reasons that we shall now explain.
There are formulas for sums of 10, 12, 14, . . . squares that are similar to Jacobi's (B1)-(B4), but they involve other more complicated terms called cusp forms [35]. Glaisher found the formulas for 10, 12, 14, 16 and 18 squares, and a general formula for any even number of squares was obtained by Ramanujan. The formulas for sums of 10 and 12 squares are For a statement of the general formula, see [34] (p. 214). A proof of the general formula and references to other proofs can be found in ref. [36].
There is no simple formula for the coefficients in the expansions of E 10 (q) or E 12 (q), but they satisfy some remarkable properties. For example, if we write where the coefficients e 12 (n) are as above. It was known to Ramanujan that the Dirichlet series can be factored, and hence we obtain the formula The formula (B20) is the analogue of Zucker's formulas for the 12-dimensional lattice. Similar formulas can be given for sums of 2k squares for any positive integer k. The number of cusp forms is (k − 1)/4 . In particular, there are no cusp forms for 1 ≤ k ≤ 4 corresponding to Zucker's formulas for the lattice sums in 2, 4, 6 or 8 dimensions; there is one cusp form for 5 ≤ k ≤ 8 corresponding to the lattice sums in 10, 12, 14 or 16 dimensions; and there are two cusp forms for 9 ≤ k ≤ 12 corresponding to the lattice sums in 18, 20, 22 or 24 dimensions.
As a consequence of Ramanujan's conjectures and Deligne's proofs, we now know that the number of representations of N as a sum of an even number 2k squares is given by a dominant term that involves a sum of the (k −1)th powers of the divisors of N , plus a correction term (the coefficient in a cusp form) that is roughly the square root in magnitude of the dominant term. When the number of squares is 2, 4, 6 or 8 there is no cusp form, and the divisor sum formula is exact, and that is the reason the formulas of Zucker exist. When the number of squares is 10, 12, 14, . . ., there is an increasing number of cusp forms, and there is no easy formula for the coefficients in their power series expansions. That is the reason why Zucker's formulas stop at 8 dimensions, and why there are no similar formulas for dimensions 10, 12, 14, . . ..