Compressed lattice sums arising from the Poisson

Compressed lattice sums arising from the Poisson equation Dedicated to Professor Hari Sirvastava D H Bailey1 and J M Borwein2 1

Lawrence Berkeley National Laboratory, Berkeley, CA 94720,

E-mail: [email protected] Research supported in part by the Director, Office of Computational and Technology Research, Division of Mathematical, Information, and Computational Sciences of the U.S. Department of Energy, under contract number DE-AC02-05CH11231. 2

Laureate Professor and Director, Centre for Computer Assisted Research Mathematics and its Applications (CARMA), University of Newcastle, Callaghan, NSW 2308, Australia. E-mail: [email protected],[email protected] Research supported in part by the Australian Research Council.

Compressed lattice sums arising from the Poisson equation

2

Abstract: Recently, attention has been directed to the problem of solving the Poisson equation, either in engineering scenarios (computational) or in regard to crystal structure (theoretical). In [1] we studied a class of lattice sums that amount to solutions of Poisson’s equation. By virtue of striking connections with Jacobi ϑ-function values, we were able to develop new closed forms for certain solutions, and to extend such analysis to related lattice sums. We also alluded to results for the compressed sum √ 1 X cos(πmx) cos(πn d y) , (1) φ2 (x, y, d) := 2 π m,n∈O m2 + dn2 where d > 0, x, y are real number and O denotes the odd integers. In this paper we first survey the earlier work and then discuss the sum (1) more completely. 1. Introduction In [1], we analyzed various generalized lattice sums [6], which have been studied for many years in the mathematical physics community, for example in [6, 12, 13]. More recently interest was triggered by some intriguing research in image processing techniques [8]. These developments have underscored the need to better understand the underlying theory behind both lattice sums and the associated Poisson potential functions. In recent times, attention was directed to the problem of solving the Poisson equation, either in engineering scenarios (computationally, say for image enhancement) or in regard to crystal structure (theoretically). In [1] we thus studied a class of lattice sums that amount to Poisson solutions, namely the n-dimensional forms φn (r1 , . . . , rn ) =

1 π2

X m1 ,...,mn

eiπ(m1 r1 +···+mn rn ) . 2 2 m + · · · + m 1 n odd

By virtue of striking connections with Jacobi ϑ-function values, we were able to develop new closed forms for certain values of the coordinates rk , and extend such analysis to similar lattice sums. A primary result was that for rational x, y, the natural potential φ2 (x, y) is π1 log A where A is an algebraic number. Various extensions and explicit evaluations were given. We also touched on results for the compressed sum √ 1 X cos(πmx) cos(πn d y) φ2 (x, y, d) := 2 . (2) π m2 + dn2 2 m,n∈O

In this paper we discuss this (2) in more detail. More detailed motivation can be found in [1, 8]. Such work is made possible by number-theoretical analysis, symbolic computation and experimental mathematics, including extensive numerical computations using up to 20,000digit arithmetic.


3

1.1. About this paper In Section 2, we describe, consistent with [8], the underlying equations along with “natural” Madelung constants and relate them to the classical Madelung constants. In Section 3, we produce the solution φn which, especially with n = 2, provide the central objects of our study. In Section 4, we describe rapid methods of evaluating φn . In Section 5 we discuss closed forms for φ2 (x, y). In Section 6 we discuss closed forms for related sums in two dimensions. This leads to the introduction of theta function formalism and allows us to resolve much else. In Section 7 we are then able to evaluate compressed potential sums. Finally, we further discuss computational matters in Section 8. 2. Madelung type sums In a recent treatment of ‘natural’ Madelung constants [8] it is pointed out that the Poisson equation for an n-dimensional point-charge source, ∇2 Φn (r) = − δ(r),

(3)

gives rise to an electrostatic potential—we call it the bare-charge potential—of the form Cn Γ(n/2 − 1) 1 =: n−2 , if n 6= 2, (4) n/2 n−2 4π r r 1 Φ2 (r) = − log r =: C2 log r, (5) 2π where r := |r|. Since this Poisson solution generally behaves as r2−n , the previous work [8] defines a “natural” Madelung constant Nn as (here, m := |m|): Φn (r) =

0

X (−1)1·m Nn := Cn , n−2 m n m∈Z

if n 6= 2,

(6)

0

N2 := C2

X

(−1)1·m log m,

m ∈ Zn

where in cases such as this log sum one must infer an analytic continuation [8], as the literal sum is quite non convergent. This Nn coincides with the classical Madelung constant 0

X (−1)1·m Mn := m m ∈ Zn only for n = 3 dimensions, in which case N3 = obvious M, N relation.

(7) 1 M3 . 4π

But in all other dimensions there is no


4

A method for gleaning information about Nn is to contemplate the Poisson equation with a crystal charge source, modeled after NaCl (salt) in the sense of alternating lattice charges: X ∇2 φn (r) = − (−1)1·m δ(m − r). (8) m ∈ Zn

Accordingly—based on the Poisson equation (3)—solutions φn can be written in terms of the respective bare-charge functions Φn , as X φn (r) = (−1)1·m Φn (r − m). (9) m ∈ Zn

3. The crystal solutions φn In [8] it is argued that a solution to (8) is Qn 1 X k=1 cos(πmk rk ) , φn (r) = 2 π m ∈ On m2

(10)

where O denotes the odd integers (including negative odds). These φn do give the potential within the appropriate n-dimensional crystal, in that φn vanishes on the surface of the cube [−1/2, 1/2]n , as is required via symmetry within an NaCl-type crystal of any dimensions—thus we have solved the Poisson equation subject to Dirichlet boundary conditions. To render this representation more explicit and efficient, we could write equivalently φn (r) =

2n π2

X m1 ,...,mn > 0,

cos(πm1 r1 ) · · · cos(πmn rn ) . 2 2 m + · · · m 1 n odd

It is also useful that—due to the symmetry inherent in having odd summation indices—we can cavalierly replace the cosine product in (10) with a simple exponential: 1 X eiπm·r φn (r) = 2 . (11) π m ∈ On m2 Q Q This follows from the simple observation that exp(πimk rk ) = (cos(πmk rk ) + i sin(πmk rk )), so when the latter product is expanded out, the appearance of even a single sin term is annihilating, due to the bipolarity of every index mk . We observe that convergence of these conditionally convergent sums is by no means obvious but that results such as [6, Thm. 8.3 & Thm. 8.5] ensure that Qn 1 X k=1 cos(πmk rk ) φn (r, s) := 2 , (12) π m ∈ On m2s is convergent and analytic with abscissa σ0 for (n − 1)/4 ≤ σ0 ≤ (n − 1)/2 where Re(s) = σ. For the central case herein, summing over increasing spheres is analytic in two dimensions


5

for σ0 ≤ 23/73 < 1/2 and in three dimensions for σ0 ≤ 25/34 < 1, but in general the best estimate we have is σ0 ≤ n/2 − 1, so for n ≥ 5 to avoid ambiguity we work with the analytic continuation of (12) from the region of absolute convergence with σ > n/2. Indeed, all our transform methods are effectively doing just that. 4. Fast series for φn From previous work [8] we know a computational series Qn−1 cos(πRk rk+1 ) 1 X sinh(πR(1/2 − |r1 |) k=1 , φn (r) = 2π R cosh(πR/2) n−1

(13)

R∈O

suitable for any nonzero vector r ∈ [−1/2, 1/2]n . The previous work also gives an improvement in the case of n = 2 dimensions, namely the following form for which the logarithmic singularity at the origin has been siphoned off: 1 cosh(πx) + cos(πy) 2 X cosh(πmx) cos(πmy) φ2 (x, y) = log − . (14) 4π cosh(πx) − cos(πy) π m(1 + eπm ) + m∈O

These series, (13) and (14) are valid, respectively, for r1 , x ∈ [−1, 1]. For clarification, we give here the (n = 3)-dimensional case of the fast series: p π 2 2 2 X sinh 2 p + q (1 − 2|x|) cos(πpy) cos(πqz) p φ3 (x, y, z) = . p π p,q > 0, odd p2 + q 2 cosh π p2 + q 2

(15)

2

Though it may not be manifest in this asymmetrical-looking series, it turns out that for any dimension n the φn (r1 , · · · rn ) is invariant under permutations and sign-flips. For example, φ3 (x, y, z) = φ3 (−y, z, −x) and so on. It thus behooves the implementor to consider x—which appears only in the first sum of (15)—to be the largest in magnitude of the three coordinates, √for optimal convergence. A good numerical test case is the exact evaluation φ3 (1/6, 1/6, 1/6) = 4π3 , which we have confirmed to 500 digits. 5. Closed forms for the φ2 potential Provably we have the following results which were established by factorization of lattice sums after being empirically discovered by the methods described in the next few sections. Theorem 1 ([1]). We have √ 1 log(1 + 2/ 3), 8π √ 1 φ2 (1/4, 1/4) = log(1 + 2), 4π φ2 (1/3, 1/3) =

(16) (17)

Compressed lattice sums arising from the Poisson equation √ 1 log(3 + 2 3). 8π Using the integer relation method PSLQ [4] to hunt for results of the form, φ2 (1/3, 0)

=

?

exp (πφ2 (x, y)) = α,

6 (18)

(19)

for α algebraic we may obtain and further simplify many results such as: Conjecture 2 ([1]). We have discovered and subsequently proven 1 α + 1/α √ ? φ2 (1/4, 0) = log α where = 2+1 4π 2 q √ √ 1 ? φ2 (1/5, 1/5) = log 3 + 2 5 + 2 5 + 2 5 , 8π γ + 1/γ √ = 3+1 2 √ τ − 1/τ 1 ? log τ where = (2 3 − 3)1/4 φ2 (1/3, 1/6) = 4π 2 p √ ! 1+ 2− 2 1 ? √ log φ2 (1/8, 1/8) = 4 4π 2−1 q √ √ 1 µ + 1/µ ? φ2 (1/10, 1/10) = log µ where =2+ 5+ 5+2 5 4π 2 ?

φ2 (1/6, 1/6) =

1 log γ 4π

where

?

where the notation = indicates that we originally only had experimental (extreme-precision numerical) evidence of an equality. Such hunts are made entirely practicable by (14). Note that for general x and y we have φ2 (y, x) = φ2 (x, y) = −φ2 (x, 1 − y), so we can restrict searches to 1/2 > x ≥ y > 0. Remark 3 (Algebraicity). In light of our current evidence we conjectured that for x, y rational, log α (20) π for α algebraic. Theorem 5 proved this conjecture. We note that Theorem 5 proves that all values should be algebraic but does not, a priori, establish the precise values we have found. This will be addressed in §6.2. 3 ?

φ2 (x, y) =


7

6. Madelung and “jellium” crystals and Jacobi ϑ-functions We have studied φ2 (x, y) :=

1 X cos(πax) cos(πby) , π 2 a,b ∈ O a2 + b 2

(21)

as a ‘natural” potential for n = 2 dimensions in the Madelung problem. Here E denotes the even integers. There is another interesting series, namely 0

1 X cos(2πax) cos(2πby) ψ2 (x, y) := 2 4π a,b ∈ Z a2 + b 2 =

1 X cos(πax) cos(πby) , π 2 a,b ∈ E a2 + b 2

(22)

where E denotes the even integers. Now it is explained, and illustrated in [8] that this ψ2 function is the ‘natural” potential for a classical jellium crystal and relates to Wigner sums [6]. This involves a positive charge at every integer lattice point, in a bath—a jelly–of uniform negative charge density. As such, the ψ functions satisfy a Poisson equation but with different source term [8]. Note importantly that ψ2 satisfies Neumann boundary conditions on the faces of the Delord cube—in contrast with the Dirichlet boundary conditions satisfied in the Madelung case. Briefly, a fast series for ψn has been worked out [8] as: Q 1 1 2 1 2n−3 X cosh(πS(1 − 2|r1 |) n−1 k=1 cos(2πSk rk+1 ) ψn (r) = + r1 − |r1 | + , 12 2 2 π S sinh(πS) +(n−1) S∈Z

(23) valid on the Delord n-cube, i.e. for r ∈ [−1/2, 1/2]n . 6.1. Closed forms for ψ2 and φ2 Using the series (23) for high-precision numerics, it was discovered (see [8, Appendix]) that previous lattice-sum literature had concealed a longtime typographical error for certain 2dimensional sums, and that a valid closed form for ψ2 is actually 1 x2 1 Γ(1/4) + log √ − log ϑ1 π(ix + y), e−π . (24) ψ2 (x, y) = 2 4π 2π 8π Γ(3/4) In the Madelung case, it then became possible to cast φ2 likewise in closed form, namely 1 φ2 (x, y) = log |α(z)| 4π

where

ϑ22 (z, q)ϑ24 (z, q) α(z) := 2 ϑ1 (z, q)ϑ23 (z, q)

for q := e−π , z := π2 (y + ix). (See [8] and [1, Appendix II] for details.)

(25)

Compressed lattice sums arising from the Poisson equation Theorem 4 ([1]). For z := π2 (y + ix) λ(z) √ 1 − ϑ2 (z, q)ϑ4 (z, q) 1 1 2 φ2 (x, y) = log = log , 1/λ(z) 2π ϑ1 (z, q)ϑ3 (z, q) 4π 1 − √2 and

8

(26)

√ 1 log 2 µ(2z) 2λ(2z) − 1 , ψ2 (x, y) = − 4π

(27)

∞ ϑ24 (z, e−π ) Y (1 − 2 cos(2z)q 2n−1 + q 4n−2 )2 = λ(z) := 2 , ϑ3 (z, e−π ) n=1 (1 + 2 cos(2z)q 2n−1 + q 4n−2 )2

(28)

where

and 2 /2

µ(z) := e−πx

∞ Y (1 + 2 cos(2z)q 2n−1 + q 4n−2 )2 ϑ23 (z, e−π ) x2 /2 = q ϑ23 (0, e−π ) (1 + q 2n−1 )4 n=1

(29)

with q := e−π . We recall the general ϑ-transform giving for all z with Re t > 0 p 2 (30) ϑ3−k (πz, e−tπ ) = 1/t e−πz /t ϑ3+k (iπz/t, e−π/t ) p 2 for k = −1, 0, 1 (while ϑ1 (πz, e−tπ ) = −1/t e−πz /t ϑ1 (iπz/t, e−π/t )). In particular with t = 1 we derive that 2

ϑ3−k (π(ix + y), e−π ) = e−πz ϑ3+k (π(iy − x), e−π ),

(31)

which directly relates |µ(π(y + ix))| and |µ(π(x + iy))| in (29). Let κ(z) :=

ϑ22 (z, e−π ) . ϑ23 (z, e−π )

(32)

Then

√ κ(ix + y) + λ(ix + y) = 2 √ √ or 1 − 2 κ(ix + y) = 2 λ(ix + y) − 1, and (30) then shows √ λ(ix + y) = κ(−x + iy) = 2 − λ(−x + iy).

(33)

(34)

Likewise (26) is unchanged on replacing λ by κ. Hence, it is equivalent to (20) to prove that for all z = π2 (y + ix) with x, y rational λ(z) in (28) is algebraic. Theorem 5 (Algebraic values of λ and µ, [1]). For all z = π2 (y+ix) with x, y rational the values 1 of λ(z) and µ(2z) in (28) are algebraic. It follows that φ2 (x, y) = 4π log α with α algebraic. 1 Similarly ψ2 (x, y) = 4π log β for β algebraic. √ This type of analysis also work for any singular value: τ = −d and q = exp(2πiτ ) and so applies to sums with n2 + d m2 in the denominator, as we see in the next section.


9

6.2. Explicit equations for degree 2, 3 and 5 We illustrate the complexity of Ωm (u, v), the algebraic polynomial linking the input and output as we move from u = z to v = mz by first considering m = 2 and m = 3. Example 6 (λ(2z) and λ(3z)). As described in [1] ϑ44 (z, e−π ) − ϑ41 (z, e−π ) , ϑ43 (z, e−π ) + ϑ41 (z, e−π )

(35)

and so letting τ (z) = ϑ21 (z, e−π )/ϑ23 (z, e−π ) we obtain √ τ (z) = 2λ(z) − 1 2 2 λ (z) − τ 2 (z) 3/2 λ(2z) = 2 . 1 + τ 2 (z)

(36)

λ(2z)1/2 = 23/4

(37)

Iteration yields Ω2n . From this one may recursively compute λ(z/2n ) from λ(z) and watch the tower of radicals grow. Correspondingly, (λ (2 z) λ (z) − τ (2 z) τ (z))2 , λ(3z) = 2 (1 + τ (2 z) τ (z))2 λ (z)

(38)

where τ (z), τ (2z), λ(2z) are given by (36) and (37). After simplification we obtain !2 √ λ4 (z) − 6 λ2 (z) + 6 2 λ (z) − 3 √ . λ (3 z) = λ (z) 3 λ4 (z) − 6 2 λ3 (z) + 6 λ2 (z) − 1

(39)

√

1−λ/ 2 √ this allows us computationally In addition, since we have an algebraic relation α = 1−1/λ/ 2 to prove the evaluation of φ2 (x/2, y/2) and φ2 (x/3, y/3) once φ2 (x, y) is determined as it is for the cases of Theorem 1. 3

We may perform the same work inductively for d = 2n + 1 using (n + 1)z ± nz in the known addition formulas for ϑk [1] to obtain Ω5 and so on. The inductive step is 2 2 λ ((n + 1) z) λ (n z) − τ ((n + 1) z) τ (n z) λ((2n + 1) z) = . (40) λ (z) 1 + τ ((n + 1) z) τ (n z) Example 7 (Empirical computation of Ωd ). Given z and u = λ(z), v = λ(dz), we compute uj v k for 0 ≤ j ≤ J, 0 ≤ k ≤ K up to degree J, K and look for a relation to precision D. This is a potential candidate for Ωd . If Ωd (λ(w), λ(dw)) ≈ 0 at a precision significantly greater than D and for various choices of w, we can reliably determine Ωd this way [1]. 3 Example 8. We conclude this subsection by proving the following empirical evaluation 1 ? φ2 (1/5, 2/5) = log 51/4 . (41) 4π

Compressed lattice sums arising from the Poisson equation Let w =

π (1 10

10

+ 2i) so that

log |α(w)| log |α(2w)| while − φ2 (2/5, 1/5) = 4π 4π since φ2 (2/5, 4/5) = −φ2 (2/5, 1/5). We may—with help from a computer algebra system—solve for the stronger requirement that α(2w) = −α(w) using Example 6. We obtain r 2i − 1 α(2w) = 5 φ2 (2/5, 1/5) =

so that |α(2w)| = 5−1/4 and |α(w)| = 51/4 , as required. To convert this into a proof we can make an a priori estimate of the degree and length of α(w)—using (37), (39) and (40) while λ(5w) = 1—and then perform a high precision computation to show no other algebraic number could approximate the answer well enough. The underlying result we appeal to [5, Exercise 8, p. 356] is given next in Theorem 9. In this particular case, we use P (α) = α4 + 2 α2 + 5 with ` = 8, d = 4 and need to confirm that |P (α)| < 5D/4 L−3 8−D . A very generous estimate of L < 1012 and D < 103 shows it is enough to check |P (α(w)| ≤ 10−765 . This is very easy to confirm. Relaxing to L < 10100 , D < 104 requires verifying |P (α(w)| ≤ 10−7584 . This takes only a little longer. 3 Recall that the length of a polynomial is the sum of the absolute value of the coefficients. Theorem 9 (Determining a zero). Suppose P is an integral polynomial of degree D and length L. Suppose that α is algebraic of degree d and length `. Then either P (α) = 0 or max{1, |α|}D . Ld−1 `D Example 10 (λ(5z)). Making explicit the recipe for λ(5z) we eventually arrive at: |P (α)| ≥

√ √ 2 λ4 (z) − 4 2λ3 (z) + 14λ2 (z) − 10 2λ (z) + 5 λ (5z) = λ (z) (42) √ √ 2 5λ4 (z) − 10 2λ3 (z) + 14λ2 (z) − 4 2λ (z) + 1 √ √ √ √ 2 λ8 (z) + 4 2λ7 (z) − 32λ6 (z) + 36 2λ5 (z) − 34λ4 (z) + 4 2λ3 (z) + 8λ2 (z) − 4 2λ (z) + 1 × √ √ √ √ 2 . λ8 (z) − 4 2λ7 (z) + 8λ6 (z) + 4 2λ5 (z) − 34λ4 (z) + 36 2λ3 (z) − 32λ2 (z) + 4 2λ (z) + 1

This shows how generous our estimates were in Example 8.

3

In similar fashion we can now computationally confirm all of the exact evaluations in Conjecture 2 and such like. For example, we know that |α(π(1/3 + i/3))| = |1/α(π/2(2/3 + √ 2i/3))|. This solves to produce −α(π/2(1/3 + i/3))2 = 1 + 2/3 3 and establishes (16) of Theorem 1. Now Example 6 can be used to produce α(π/2(1/6 + i/6)) is as given in Conjecture 2, and so on. We now turn to our main focus:


11

7. “Compressed” potentials Yet another solution to the Poisson equation with crystal charge source, for d > 0, is √ 1 X cos(πmx) cos(πn d y) φ2 (x, y, d) := 2 . (43) 2 + dn2 π m m,n∈O2 √ This is the potential inside a crystal compressed by 1/ d on the y-axis, √ sense that √ in the the Delord-cube now becomes the cuboid {x, y} ∈ [−1/2, 1/2] × [−1/(2 d), 1/(2 d)]. Indeed, φ2 (x, y, d) vanishes on the faces of this 2-cuboid (rectangle) for d > 0 and integer. This is technically only a compression when d > 1. Along the same lines as involve the analysis of (13), we can posit a fast series √ √ 1 X sinh(πR d (1/2 − |x|) cos(πR d y) √ φ2 (x, y, d) = √ , (44) π d R ∈ O+ R cosh(πR d/2) √ √ valid on the Delord cuboid, i.e. x ∈ [−1/2, 1/2], y ∈ [−1/(2 d), 1/(2 d)]. Moreover, the log-accumulation technique of [8, Appendix] applied as with Theorem 4 yields: Theorem 11 (ϑ-representation for compressed potential). We have ϑ2 (z, q)ϑ4 (z, q) 1 , φ2 (x, y, d) = √ log ϑ1 (z, q)ϑ3 (z, q) 2 dπ √ √ where q := exp(−π d) and z := 12 π d(y + ix).

(45)

This theorem is consistent with d = 1, in the sense φ2 (x, y, 1) = φ2 (x, y). In addition, for integers a, b ≥ 1 we have: X b2 1 cos(πmx) cos(πny) 2 a φ2 ax, ay, 2 := 2 . (46) a π m∈aO,n∈bO m2 + n2 We note that there are various trivial evaluations such as φ2 (x, 1/4, 4) = 0 for all x. √ Moreover, for each positive rational d there is an analogue of Theorem 5 in which 1/ 2 is √ √ √ replaced by the d-th singular value kd [5, 6]. Thence, k2 = 2 − 1, k3 = ( 3 − 1)/ 8, and √ k4 = ( 2 − 1)2 . Precisely, we set √ √ ϑ24 (z, exp(−π d)) ϑ22 (z, exp(−π d)) √ , √ λd (z) := 2 κd (z) := 2 ϑ3 (z, exp(−π d)) ϑ3 (z, exp(−π d)) √ ϑ21 (z, exp(−π d)) √ τd (z) := 2 (47) ϑ3 (z, exp(−π d)) and have kd0 λd (z) + kd κd (z) = 1,

kd λd (z) − kd0 κd (z) = τd (z).

(48)


12

(see [5, Prop. 2.1]). We observe that (45) can be written for all x, y as 1 − kd0 λd (z) 1 , φ2 (x, y, d) = √ log 0 1 − k /λ (z) 4 dπ d

(49)

d

on appealing to (47) and (48). Hence we may now study λd exclusively. Fix an integer m > 1. The addition formulas for the ϑ’s as given in [5, §2.6] ϑ3 (z + w, q)ϑ3 (z − w, q)ϑ23 (0, q)

= ϑ23 (z, q)ϑ23 (w, q) + ϑ21 (z, q)ϑ21 (w, q) (50)

ϑ4 (z + w, q)ϑ4 (z − w, q)ϑ24 (0, q)

= ϑ24 (z, q)ϑ24 (w, q) − ϑ21 (z, q)ϑ21 (w, q) (51)

ϑ1 (z + w, q)ϑ1 (z − w, q)ϑ2 (0, q)ϑ3 (0, q) = ϑ1 (z, q)ϑ4 (z, q)ϑ2 (w, q)ϑ3 (w, q) + ϑ1 (w, q)ϑ4 (w, q)ϑ2 (z, q)ϑ3 (z, q)

(52)

when replacing z by (m − 1)z and w by z, and appealing to (48), allow one recursively to write λd (mz) algebraically in terms of λd (z). We first give the equation for λd (2z). Example 12 (λd (2z)). As for λ(2z), which is the case d = 1, we have analytically λd (z) − kd0 kd 2 2 λd (z) − τd2 (z) 0 −3 λd (2z) = kd . 1 + τd2 (z)

(53)

τd (z) =

(54)

This becomes λd (2z) =

kd0

1 + λ2d (z) − 2λd (z)/kd0 1 + λ2d (z) − 2λd (z) · kd0

2 ,

(55) 3

and λd (0) = kd0 . We turn to λd (3z). Example 13 (λd (3z)). As for λ(3z), which is the case d = 1, we have λd (z) − kd0 λd (2z) − kd0 , τd (2z) = , kd kd 2 λd (z)λd (2z) − τd (z)τd (2z) 1 λd (3z) = 0 2 . 1 + τd (z)τd (2z) k d λd (z)

τd (z) =

(56) (57)

In tandem with (55) this becomes λd (3z) = λd

4 λd − 3 kd0 − 6 λ3d kd0 + λ4d kd0 + 4 λd kd0 2 4 λ3d + kd0 − 6 λ2d kd0 − 3 λ4d kd0 + 4 λ3d kd0 2

2 2 ,

(58)

since λd (0) = kd0 . This can also be determined by using the methods of Example 7 to obtain the algebraic equation linking u = λd (z) and v = λd (3z) and then carefully factoring the result—guided by


13

results such as Example 10. It could also have been rigorously proven by the method used for corresponding formula with d = 1 was in [1]. 3 We next exhibit several sample compressed-sum evaluations of which (60) is the most striking: Theorem 14 (Some explicit ! √ 1 2 , ,2 φ2 3 3 ! √ 1 3 , ,3 φ2 3 3 1 1 φ2 , ,4 3 3 1 1 φ2 , ,9 4 4

compressed sums). √ 1 = √ log 9 − 6 2 8 2π 1 1 = √ log 3 8 3π √ √ √ √ 2 1 = log 3 2+ 3 2− 3 16π 3 √ √ √ √ 2 1 1 = log 3+5 3−4 2 3−1 1+ 2 48π 4 √ √ 2 3− 2 · .

(59) (60) (61)

(62)

√ √ Proof. Let us first sketch a proof of (60). This requires showing α = exp(8π 3ψ2 (1 3, 1/3, d)) = 1/3. Direct computer algebra using (49) shows that it suffices to prove that λ3 (z) = λ3 (z) = √ √ √ 1/ 2 + 1/ 6 when z3 := π(1/2 + i 3/6). Now (58) becomes 2 λ3 4 λ3 − 3 k30 − 6 λ3 2 k30 + λ3 4 k30 + 4 λ3 k30 2 λ3 (3z) = (63) 2 . 4 λ3 3 + k30 − 6 λ3 2 k30 − 3 λ3 4 k30 + 4 λ3 3 k30 2 √ However 3z is in the compressed lattice generated by 1, i/ 3 and so ) ( √ ! √ √ √ 3 + 1 3 3 3 − 1 √ = . λ3 (3z) ∈ 0, ∞, k30 = √ , λ3 iπ 2 2 2 2 2 2 From this, we may deduce that λ3 (z) is a zero of the denominator of (63), whose only real root √ √ is 1/ 2 + 1/ 6. √ To prove (59) we similarly consider (49) and (58) with d = 2 so that k2 = 2 − 1 and √ z2 := π(1 + i/ 2)/3 is in the compressed lattice. Indeed, λ2 (3z) = 0 and we now deduce that p p √ √ 2 + 2 2 + i 10 2 − 14 λ2 (z) = 2 and the result again follows. To prove (61) requires an application of this method with d = 4 and x = 1/3 and y = 1/6. We likewise may prove (62) (with d = 9). In this final case x = 1/4 and y = 1/12. In these cases we need to appeal to both (57) and (58) at least once.

Compressed lattice sums arising from the Poisson equation 14 √ √ In general the polynomial for α = exp(8π dψ2 (x, y d, d)) with x, y rational can be calculated as we describe in the next section. √ 7.1. Further results for φ2 (a/c, b/c d, d) √ As noted, we anticipated similarly algebraic evaluations for φ2 (a/c, b/c d, d) for nonnegative integer values of a, b, c, d. And, as we will show in more examples in the next section, this is indeed true. All such results can, in principle, be rigorously established by the techniques of Example 7 and Example 8 or indeed Theorem√14. But we forfend. We do however give the first few polynomials, Pd , for exp(8πφ2 (1/d, 1/ d, d)), that is, with a = b = 1, c = d: √ √ Example 15 (Some explicit compressed polynomials for exp(8π dφ2 (1/d, 1/ d, d))). We present here the polynomials Pd when 3 ≤ d ≤ 9: d 3 4 5 6 7 8

minimal polynomial 1 − 3α 1 − 644α + 3334α2 − 644α3 + α4 1 + 20α − 210α2 + 100α3 + 25α4 1 − 39368α + 309916α2 + 589192α3 + 639814α4 + 589192α5 + 309916α6 − 39368α7 + α8 1 − 14α + 21α2 − 7α3 1 − 972816α + 108957816α2 − 3146762800α3 + 87457048348α4 − 579962120464α5 −2624405729464α6 − 52677445809328α7 + 115992837288518α8 − 52677445809328α9 −2624405729464α10 − 579962120464α11 + 87457048348α12 − 3146762800α13 +108957816α14 − 972816α15 + α16 9 1 + 438α − 8997α2 + 8292α3 + 711α4 − 954α5 − 3α6

Note √ that reversing the order of the variables leads to a more intractable computation since φ2 (1/ d, 1/d, d)) is not of the form covered by by Theorem 11. 3 8. Computation techniques and results Numerous experimental evaluations of minimal polynomials satisfied by exp(8π φ2 (x, y)) (the case d = 1), and a detailed description of the underlying computational methodology employed to find these evaluations, are presented in [1]. A similar computational scheme was employed in √ √ this study to obtain minimal polynomials satisfied by exp(8π d φ2 (x, dy, d)). For instance, we recover among a myriad of other evaluations √ ! 1 1 3+2 3 1 φ2 . , ,9 = log 3 3 24π 9 This computational scheme is as follows:


15

(i) Given x, y and d, select a conjectured polynomial degree m and a precision level P . We typically set the numeric precision level P somewhat greater than m2 digits. √ (ii) Compute φ2 (x, y d, d) to P -digit precision using formula (45). Evaluate the four theta functions indicated using the very rapidly convergent formulas given in [5, pg. 52] or [14]. √ √ (iii) Generate the (m + 1)-long vector (1, α, α2 , · · · , αm ), where α = exp(8π d φ2 (x, y d, d)). Note: we have found that without the eight here, the degree of the resulting polynomial would be eight times as high (but the larger polynomials were in fourth or eighth powers). Given the very rapidly escalating computational cost of higher degrees, many of the results in [1] and herein would not be feasible without this factor. (iv) Apply the PSLQ algorithm (we employed the two-level multipair variant of PSLQ for d up to 19, and the three-level multipair variant for d ≥ 20 [2]) to find a nontrivial (m + 1)-long integer vector A = (a0 , a1 , a2 , · · · , am ) such that a0 + a1 α + a2 α2 + · · · + am αm = 0, if such a vector exists. PSLQ (or one of its variants) either finds a vector A, which then is the vector of coefficients of an integer polynomial satisfied by α (certified to the “epsilon” of the numerical precision used), or else exhausts precision without finding a relation, in which case the algorithm nonetheless provides a lower bound on the Euclidean norm of the coefficients of any possible degree-m integer polynomial A satisfied by α. (v) If no numerically significant relation is found, try again with a larger value of m and correspondingly higher precision. If a relation is found, try with somewhat lower m, until the minimal m is found that produces a numerically significant relation vector A. Here “numerically significant” means that the relation holds to at least 100 digits beyond the level of precision required to discover it. To obtain greater assurance that the polynomial produced by this process is in fact the minimal polynomial for α, use the polynomial factorization facilities in Mathematica or Maple to attempt to factor the resulting polynomial. √ √ Results that we have obtained for the special set of cases exp(8π d φ2 (1/d, 1/ d, d)), for integers d up to 40, are shown in Table 1. Our computations required up to 20,000digit precision, and, for large degrees and correspondingly high precision levels, were rather expensive (over 60 processor-hours in some cases). We employed the ARPREC arbitrary precision software [3]. √ √ The degree m(d) of the minimal polynomial for exp(8π d φ2 (1/d, 1/ d, d)), the number of zeroes z(d) among the minimal polynomial coefficients, the numeric precision level P , the run time in seconds T , and the approximate base-10 logarithm M of the absolute value of the central coefficient, are also shown in the table, together with the ratio (log10 M )/m(d). We were not able to obtain relations for all case d ≤ 40. Evidently the 12,000 or 20,000 digits we employed in these cases were still insufficient to recover the underlying relations. It is


16

interesting that in these cases, as opposed to the φ(x, y) constants we studied in [1], odd values and prime values of d generally yielded simpler relations than the even instances. In the earlier study [1], we included Jason Kimberley’s observation that the degree m(d) of the minimal polynomial for exp(8π φ2 (1/d, 1/d)) appeared to be given for odd primes by m(4k + 1) = (2k) · (2k) and m(4k + 3) = (2k + 2) · (2k + 1). If one sets m(2) = 1/2, for notational convenience, then it appears that for any prime factorization of an integer greater than 2: ! k k Y Y ? 2(e −1) m pei i = 4k−1 pi i m(pi ). (64) i=1

i=1

Unfortunately, we have not yet been able to find a similar formula for the minimal polynomial degree in the compressed cases under study in Table 1. We are actively investigating this at the present time. With regards to the run times listed in Table 1 (given to 0.01 second accuracy), it should be recognized that like all computer run times, particularly in a multicore or multiprocessor environment, they are only repeatable to two or three significant digits. They are listed here only to emphasize the extremely rapid increase in computational cost as the degree m and corresponding precision level P increase. Acknowledgements The authors thank many colleagues, for fruitful discussions about lattice sums and theta functions. [1] D.H. Bailey, J.M. Borwein, R.E. Crandall and I.J. Zucker, “Lattice sums arising from the Poisson equation.” Preprint 2102. Available at http://www.carma.newcastle.edu.au/jon/poisson.pdf. [2] D.H. Bailey and D.J. Broadhurst, “Parallel integer relation detection: Techniques and applications,” Mathematics of Computation, 70 236 (2000), 1719–1736. [3] D.H. Bailey, Y. Hida, X.S. Li and B. Thompson, “ARPREC: An arbitrary precision computation package,” Sept 2002, available at http://crd-legacy.lbl.gov/ dhbailey/dhbpapers/arprec.pdf. The ARPREC software itself is available at http://crd-legacy.lbl.gov/ dhbailey/mpdist. [4] J.M. Borwein and D.H. Bailey, Mathematics by Experiment: Plausible Reasoning in the 21st Century A.K. Peters Ltd, 2004, ISBN: 1-56881-136-5. Second expanded edition, 2008. [5] J.M. Borwein and P. B. Borwein, Pi and the AGM, John Wiley & Sons, New York, 1987. [6] J. Borwein, L. Glasser, R. McPhedran, J. Wan and J. Zucker, Lattice Sums: Then and Now. Cambridge University Press. To appear, 2013. [7] R. Crandall, “Unified algorithms for polylogarithms, L-functions, and zeta variants,” in Algorithmic reflections: Selected works, PSIpress, Portland, 2011. Also downloadable at http://www.psipress.com. [8] R. Crandall, “The Poisson equation and ‘natural’ Madelung constants,” preprint (2012). [9] P.J. Forrester and M.L. Glasser, “Some new lattice sums including an exact result for the electrostatic potential within the NaCl lattic,” J. Phys. A 15 (1982), 911–914. [10] M.L. Glasser,“The evaluation of lattice sums. III. Phase modulated sums,” J. Math. Phys. 15 188 (1974), 188–189.

Compressed lattice sums arising from the Poisson equation d m(d) 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40

1 4 4 8 3 16 6 16 15 16 12 48 8 32 32 24 27 64 24 40 33 64 20 128 27 48 84 64 45 128 40 72 72 72 36 72 48 72

P

T

200 0.00 200 0.01 200 0.01 200 0.02 200 0.01 1000 0.85 200 0.01 1000 0.76 1000 0.55 1000 0.85 1000 0.45 4000 100.31 200 0.03 2000 12.51 2000 8.93 2000 4.32 2000 5.27 5000 415.78 2000 4.52 3000 64.33 3000 10.92 5000 412.59 2000 3.77 20000 failed 2000 5.70 4000 131.29 6000 1375.38 5000 557.29 4000 38.89 20000 failed 4000 81.42 12000 failed 12000 1179.43 12000 failed 3000 29.11 12000 failed 4000 56.33 12000 failed

17

log10 M

log10 M m(d)

3.5230 2.3222 5.8061 1.1461 14.0644 3.9187 12.1262 7.5230 15.2686 9.5690 38.4539 2.1139 33.7985 21.9952 21.1366 16.8591 58.8250 15.6374 41.4566 14.4705 60.3300 20.0766

0.8807 0.5806 0.7258 0.3820 0.8790 0.6531 0.7579 0.5015 0.9543 0.7974 0.8011 0.2642 1.0562 0.6874 0.8807 0.6244 0.9191 0.6516 1.0364 0.4385 0.9427 1.0038

18.9234 58.4901 60.0921 56.9952 19.8425

0.7009 1.2185 0.7154 0.8906 0.4409

32.1363

0.8034

41.3569

0.5744

43.5933

1.2110

20.7849

0.4330

√ √ Table 1. PSLQ runs to recover minimal polynomials satisfied by exp(8π d φ2 (1/d, 1/ d, d)). Here m(d) is the degree, P is the precision level in digits, T is the run time in seconds, and log10 M is the size in digits of the central coefficient.


18

[11] M.L. Glasser and I.J. Zucker,“Lattice Sums,” pp. 67–139 in Perspectives in Theoretical Chemistry: Advances and Perspectives, Vol. 5 (Ed. H. Eyring). New York: Academic Press, 1980. [12] S. Kanemitsu, Y. Tanigawa, H. Tsukada and M. Yoshimoto, “Crystal symmetry viewed as zeta symmetry,” pp. 91–130 in T. Aoki, S. Kanemitsu, M. Nakahara and Y. Ohno, eds., Zeta Functions, Topology and Quantum Physics, Springer, New York, 2005 (reprinted 2010). [13] S. Kanemitsu and H. Tsukada, “Crystal symmetry viewed as zeta symmetry II,” pp. 275–292 in K. Alladi, J.R. Klauder and C.R. Rao, eds., The Legacy of Alladi Ramakrishnan in the Mathematical Sciences, Springer, New York, 2010. [14] E. Weisstein, “Jacobi theta functions,” MathWorld, available at http://mathworld.wolfram.com/JacobiThetaFunctions.html. [15] E.T. Whittaker and G.N. Watson, A Course of Modern Analysis, 4th edition, Cambridge University Press, 1946. [16] I.J. Zucker, “Exact values for some lattice sums in 2, 4, 6 and 8 dimensions,” J. Phys. A. 7 (1974), 1568–1575 . [17] I.J. Zucker, “Functional equations for poly-dimensional zeta functions and the evaluation of Madelung constants,” J. Phys. A., 9 (1976), 499–505. [18] I.J. Zucker, “Some infinite series of exponential and hyperbolic functions,” SIAM J. Math., 15, 2 (1984) , 406–413. [19] I.J. Zucker, “A systematic way of converting infinite series into infinite products,” J. Phys. A., 20 (1987), L13-L17. [20] I.J. Zucker, “Further relations amongst infinite series and products: II. The evaluation of three-dimensional lattice sums,” J. Phys. A., 23 (1990), 117–132. [21] I.J. Zucker and R.C. McPhedran. “Dirichlet L-series with real and complex characters and their application to solving double sums,” Proc.Roy.Soc.Lond. A 464 (2008), 1405–1422.