An asymptotic expansion for the discrete harmonic potential

We give two algorithms that allow to get arbitrary precision asymptotics for the harmonic potential of a random walk.


INTRODUCTION
The discrete harmonic potential a of a random walk R, starting from 0, on a lattice Z can be most easily defined using where P n (z) is the probability of R to be at z at the nth step. It is easy to conclude from (1) that ∆a = δ {0} where ∆ is the discrete Laplacian for the walk R. Less obvious are the facts that a is positive, and that these two properties determine a up to the addition of a positive constant. "Positive" may be replaced by "with sub-linear growth" in dimension greater than 1. The discrete harmonic potential is an interesting quantity, strongly related to the discrete Green function g(x, y): on the one hand, a(x) = g(x, x) where g is the Green function corresponding to the set {0} [S76, P11.6, page 118]. On the other hand, g for arbitrary sets can be calculated if a is known -when the set is finite there is an explicit formula, [S76, T14.2, page 143].
Thus it becomes interesting to calculate or estimate a. Except in dimension 1 where a(x) = |x|, the problem is far from trivial. Take as an example the harmonic potential of the regular random walk on Z 2 which has a logarithmic nature. The estimate of a most common in the literature is where λ is some number, which can be expressed using Euler's γ constant, λ = 2 π γ + 1 π log 8. The earliest proof of (2) we are aware of is in Stöhr [S49]. A more accessible proof of a weaker result, giving an error estimate of o(1), can be found in Spitzer [S76, P12.3, page 124]. Spitzer's proof (which is not specific to the lattice Z 2 ), like Stöhr's, relies on the fact that the Fourier transform P n has an explicit formula, which can be summed to give a "pseudo-Fourier" representation of a as the explicit integral a(z) = 1 − e i ξ,z 1 − P 1 (ξ) dξ ; ( 3) indeed, this is how he proves that the sum in (1) converges in the first place. Fukai and Uchiyama used this technique to find the next order approximation -in the case of the simple random walk on Z 2 it is − Re z 4 6π|z| 6 -and to show that an asymptotic polynomial expansion exists. See [FU96] for the two dimensional case and [U98] for the case d ≥ 3.
The purpose of this note is to give two alternative approaches to the computation of a, which allow to get asymptotics of arbitrary precision. For example, here are the first few terms for the regular random walk on Z 2 : .
The first approach is a brute-force one which starts from (1): P n may be represented as a sum of multinom coefficients, and plugging in Stirling's expansion, and some elementary algebra gives the representation (4). Unfortunately, the algebra involved is too long. In effect, for all but the first coefficient (the − Re z 4 6π|z| 6 ), it is simply not practical to do the calculations by hand. We demonstrate this approach in two ways: we use it to prove that a polynomial approximation exists, reproducing the results of [FU96,U98] in the case the random walk is bounded (Fukai and Uchiyama show these results under the weaker assumption that each step of the random walk has K + d + ǫ moments, where K is the required precision in the expansion, d is the dimension and ǫ > 0). See theorem 1 on page 5; and we use it to get some explicit constants and high order expansions for other walks. See section 4.
The second approach is to approximate the discrete Laplacian ∆ using an appropriate differential operator. For example, for the regular random walk on Z 2 we have, from the Taylor expansion, In fact, ∆ = 1 2 (cosh ∂ ∂x + cosh ∂ ∂y ) − I, but we will not use this representation. Given that the expansion of a has an appropriate polynomial form, the expansion (5) allows to get many equations on the coefficients in (4). These equations are of course insufficient -the equation ∆a = δ {0} does not determine a uniquely, and the condition that a is positive or slowly increasing is hard to encode into the differential equations. However these differential equations allow to prove many of the apparent properties of (4), most notably that the coefficient of z l |z| −k for l < 1 2 k is always zero. This is theorem 2. Unfortunately, it is still not clear why all the coefficients are negative.
Returning to the case of the regular random walk on Z 2 , it can be seen that the equation ∆a = δ {0} combined with the π 2 -rotational symmetry gives that the values of a are uniquely defined by their values on one diagonal, say {m + im} ∞ m=0 , and vice versa, any arbitrary sequence on the diagonal can be extended to a symmetric function f with ∆f = δ {0} using a very simple recursion. On the other hand, it turns out that the values of a on the diagonal can be calculated explicitly using the integral (3) to get [S76, chapter 15]. Indeed we have just sketched an algorithm for calculating any value of a(z): this is the McCrea-Whipple algorithm [MW40] 1 . As is now evident, all the values are of the form n + 1 π q, n ∈ Z and q ∈ Q. Interestingly, when not on the diagonal, n and q increase exponentially, even though the result is only logarithmic in size 2 . As for the coefficients of (4), it turns out -perhaps unsurprisingly -that (6) can be used to complete the few coefficients that cannot be deduced from (5). This is theorem 3.
There is a third approach to the problem which we shall not discuss at length as it seems inferior to both former ones. Roughly it goes as follows (for the case of the regular random walk on Z 2 ). "Guess" that the solution is approximately 2 π log |z|. We discretize to Z 2 in the natural way: define Actually, it is not necessary to guess the exact value 2 π , other constants not too different also work. A calculation can show that This allows us to write a series of corrections of f as follows: f 1 = f and then and hence (7) gives that the f n 's converge 3 in the L 1 norm to an a which satisfies ∆a = δ {0} and |a| ≤ C log |z| + C so it is the harmonic potential up to an additive constant. Furthermore, the fact that ∆f − δ {0} = O(|z| −4 ) can be used to derive similar estimates for a. However, the best approximation we were able to get from that method is and getting this approximation was not significantly easier than the brute force approach. Moreover, this approach did not allow to get actual values for any of the constants -only that some constants exist, and their values were derived from (6). We have described the contents of this paper except for the last section. That section contains the results of some computer-aided simulations and a calculation of an explicit constant in the O(·) in (2).
We wish to thank Greg Lawler for some useful remarks.

Standard definitions.
A lattice is a discrete additive subgroup of R n . The dimension of a lattice Z (denoted by dim Z) is the dimension of the linear span of Z as a linear subspace of R n . A basis for a d-dimensional lattice is a set {e 1 , . . . , e d } such that Z = e 1 Z + · · · + e d Z. Z ′ is called a sublattice of Z if it is a subgroup of Z, and its index is the size of Z/Z ′ . The volume of (the cell of) a lattice Z, denoted by vol Z, is the volume of R d /Z (the discreteness of Z allows us to select a measurable set of representatives). Alternatively it can be defined as where B(r, 0) is a ball of radius r around 0 and |B(r, 0)| is its volume.
A random walk R on a lattice Z is a probabilistic process {R 0 , R 1 , . . . }, R i ∈ Z, such that R 0 = 0 and R i − R i−1 has a distribution independent of i, and such that for every z ∈ Z there exists some n such that P(z = R n ) > 0. 4 It is bounded if R i −R i−1 is bounded. We will typically confuse the process R with the distribution of any step (e.g. with the distribution of R 1 ) so that we can use notations such as ER comfortably. The dimension of R is the dimension of the lattice Z, and is denoted by dim R.
The drift of a random walk R is ER. The random walk is balanced if its drift is zero.
The discrete Laplacian ∆ of a random walk R is an operator on functions on the lattice Z defined by ||v|| will always refer to the L 2 norm. Similarly, for a matrix A, ||A|| refers to its norm as an operator from L 2 to L 2 .

THE DIRECT APPROACH
The aim of this section is to prove the following theorem.
Theorem 1. Let a be the harmonic potential of a d-dimensional balanced bounded random walk on a lattice Z. Let K > 0 be some integer. Then there exists a constant τ , a linear map A, a polynomial Q = Q K in d variables and an integer L = L K such that τ and A have explicit formulas -see the comments after the end of the proof (page 10 below). The expression Q(z)/||Az|| L is λ + O(||Az|| 1−d ) (this is classical, of course, but will also follow from the proof below). The number λ is of course also walk-specific.
In the case d = 1 the equation ∆a = δ {0} deteriorates into a simple recursion formula and the theorem is nothing but an exercise (and Q ≡ λ). Thus we concentrate on the case d > 1.
The proof of theorem 1 is replete with these rational functions Q obeying the condition above. Notice that this simply means that such that c 0 0 = 1 and otherwise c i j = 0 only if −L ≤ i < 0 and 0 ≤ j 1 + · · · + j k < −2i. It turns out (and this can be deduced from the proofs below with some care) that L K = 4K is always sufficient. However we will have no use for this fact.
Proof of lemma 1. Consider Stirling's series, Assume for the moment that ||w|| < C √ n log n where C will be fixed later. In this case we can insert monomials such as This proves the case ||w|| < C √ n log n. In the case ||w|| > C √ n log n both sides of (10) are O(n −K ) which follows easily from the fact that in this case, for C sufficiently large where Q is a polynomial. The coefficients of Q are fixed polynomial functions (depending on A) of the coefficients of P and of v. Q is independent of y and the O(·) is uniform in v and y.
The use of the parameter n in the formulation of the lemma is a bit artificial. Later, when we shall use the lemma, n will be an integer number (the number of steps), but this fact is not needed for the formulation or proof of the lemma.
where r = r 1 + · · · + r d , and as usual d r This means that as n → ∞ we have The exponent inside the O(·) kills of course all the other factors. The case r = 0 immediately gives the first part of the lemma. For the second part, put ξ = 0 in (12) and get that f (0) = 0 only when r is even, and that f (0) = n d/2 Q r (n) where Q r depends only on A. Writing P = r c r (z + v) r we get Q = r c r Q r which concludes the lemma: the fact that c r are polynomial in the coefficients of P and in v gives the same for the coefficients of Q; the fact that c r = O(||v|| deg P ) gives the same in the O(·) in (11). The result does not depend on y (except in the O(·) error) since f (n) does not depend on v -in effect the dependence on v appears only via the dependence of c r on v.
Lemma 3. Let R be a d-dimensional balanced bounded random walk on the lattice Z. Let K ∈ N. Then there exists a sublattice Y ⊂ Z × Z, a constant τ ′ , a quadratic function B and a Q such that with n L Q a polynomial in n and in the coordinates of z, and Q(n, z) = 1+o(1) uniformly in ||z|| ≤ C √ n log n.
Assume for simplicity that Z ⊂ R d . A simple comparison of the expectations of the left and right side of (13) shows that B(z) = 1 2 M −1 z, z , where M is the correlation matrix of R, M ij = E R, e i R, e j , where e i are unit vectors. Notice that the expression M −1 z, z is independent of the coordinate system chosen.
The lattice Y is clearly d+1 dimensional since otherwise we would have dim R < d. Thus Y ∩ Z is d-dimensional and vol Y ∩ Z is finite. This, with the requirement P n = 1 allows to calculate τ ′ and get We will not be using these equalities in the proof, though.
Proof. By applying a linear transformation we may assume R is a random walk on Z d . Assume R moves from z to z + v i with probability p i for i = 1, . . . , N + 1. Let W = (n, w 1 , . . . , w N +1 ) : p i n + w i ∈ Z ∀i and w i = 0 which is an (N + 1)-dimensional lattice. For an n ∈ N let k i be the number of times (between 0 and n − 1) that R moved from z to z + v i . We use lemma 1 with some K 2 which will be fixed later. This gives us where B 2 is a strictly positive quadratic form, C is some constant and n L Q 2 is a polynomial with Q 2 = 1 + o(1) uniformly for ||w|| ≤ C √ n log n. Now, for a point z ∈ R d we can write Since W (n, z) are intersections of a lattice with a translated subspace, we get that W (n, z) = ∅ for (n, z) outside some lattice Y , and for (n, z) ∈ Y we get that W (n, z) are translations of one fixed lattice X. Since we know Y is d + 1 dimensional (see the comment just after the statement of this lemma) we get that X is (N − d)dimensional. The calculation of the sum in (16) will be done of course with lemma 2 but we need some preparations first. By applying a linear transformation to the N -dimensional space of w's we may assume that W (n, z) = z 1 e 1 + · · · + z d e d + e d+1 Z + · · · + e N Z + y n where y n ∈ span{e d+1 , . . . , e N }. Combining (15), (16) and the linear transformation we get where K 3 = K 2 − N + d. This last error estimate follows easily from the following obvious estimates: From now on it would be easier to replace B 3 (w) with ||Aw|| 2 for some A ∈ GL(R N ) which is possible since B 3 is also strictly positive. We denote x = (0, . . . , 0, w d+1 , . . . , w N ) so that w = x + z and we can write ||Aw|| 2 = ||Ax|| 2 + ||Az|| 2 + 2 x, A * Az .
Since A * A is strictly positive we get that its lower right N − d minor is also strictly positive and thus invertible. Therefore there exists some . . , N . This of course implies x, A * Av = x, A * Az and then The importance of (18) is that its left part, ||Az|| 2 − ||Av|| 2 , is constant for w ∈ W (n, z), while exp − 1 n ||A(x + v)|| 2 can be summed by using lemma 2 on the last N − d coordinates. Notice that it depends on the z's only via v, which is linear in z. In other words, we can now define our quadratic form B: B(z) := ||Az|| 2 − ||Av(z)|| 2 and apply lemma 2. We get We may already remark that picking K 2 = K + N − d will give us a precision of O(n −K ) as required. Obviously O(e −cn ||v|| deg Q3 ) = O(n −K ) for any K as v is linear in z and ||z|| ≤ C √ n log n. To say something about Q, we need to consider Q 3 as a polynomial in x with coefficients which are polynomials in z and negative powers in n, i.e.
where n L c r are polynomials. Lemma 2 tells us that n L Q is a polynomial in n and that its coefficients are polynomial in the coefficients of Q 3 and in v (which is linear in z). Thus we get that n L Q is polynomial in both n and z. This concludes the lemma: the only claim we haven't shown is that Q = 1 + o(1) which can be deduced, for example, from the central limit theorem 5 .
Remark. For specific walks it is possible to prove lemma 3 using much simpler techniques. For example, for the simple random walk on Z 2 , we have where we make the notational convention that n α = 0 if α is not an integer. This formula allows to get lemma 3 from lemma 1 more-or-less directly, with no need to go through lemma 2. This works for other walks as well. However it does not seem to work for the very interesting case of the random 6-walk on the triangular lattice, a walk which can be considered no less natural than the random walk on Z 2 . We mean here a two dimensional random walk that goes at each step with probability 1 6 from z to one of the points z + ω i , i = 0, . . . , 5 where ω = 1 2 + i √ 3 2 = 6 √ 1.
Lemma 4. Let q ∈ N and 1 ≤ j ≤ q. For every s > 1 and for r → ∞ we have the superpolynomial estimate while for s = 1 we have Proof. Again we use the (one dimensional) Poisson summation formula, this time for the function Simple Fourier operations show that F s,j,q = 1 qr s−1 ξ∈Z e −2πijξ/q f 2π rξ q .
Since f ∈ C ∞ we get that f (ξ) = O(ξ −K ) for all K and hence F = r 1−s q f (0) + O(r −K ). This concludes the case s > 1. For the case s = 1 we notice that the O(·) in (19) is uniform in s ∈ (1, 2] since it depends only on the L 1 norm of the K'th derivative of f . Thus we need only to calculate (1) and ζ(s) = 1 s−1 + C + o(1) and a little use of l'Hôpital rule will give log r + C and conclude the lemma. In the other direction, there are known precise results connecting the "depth" of the zero of f at 0, i.e. the rate of convergence of f to zero, with the property that f is in some weighted H p space. See theorem 1.1 in [MS02]. We omit the details. This shows, for example, that f (ξ) = O(e −ξ 1/2+ǫ ) for any ǫ > 0.
Proof of theorem 1. Let B, τ ′ and Y be the quadratic form, constant and lattice given for R and K + 1 by lemma 3. Denoting a(z) := n:(n,z)∈Y we see thatã(z) − a(z) is a constant (indeed, the precaution of adding τ ′ n −d/2 is only necessary in the case d = 2 in order to make the sum converge). The fact that R is bounded shows that P(R n = z) is zero for n < c||z|| and of course n<c||z|| e −B(z)/n = O(e −c||z|| ). Therefore we can use (13) to get a(z) = n:(n,z)∈Y Lemma 4 allows to calculate the left sum, while the right one is obviously O(||z|| −K ). This concludes the proof of theorem 1.
We see now that the quadratic form of theorem 1 is the same as that of lemma 3, i.e. B(z) = 1 2 M −1 z, z where M is the correlation matrix of R. Thus A = M −1/2 -of course, this is not uniquely defined, but the quadratic form ||Az|| 2 is. As for τ , the factor vol Y ∩ Z of lemma 3 with the term 1 q of lemma 4 gives vol Z (this is easy to see) and we get τ = vol Z (2π) d/2 | det M | 1/2 · 2 d = 2 Γ( 1 2 d − 1) d = 2 (the factor 2 in dimension 2 comes from the fact that we formulated µ 2 = log ||z|| rather than log ||z|| 2 ).

THE DIFFERENTIAL APPROACH
We shall demonstrate the differential approach by proving a few results that would seem quite difficult using the direct approach only.
For simplicity of notation we shall always assume that the correlation matrix of R is constant. Thus we call a random walk spherical if it is bounded, balanced, and E R, e i R, e j = Cδ(i, j) (again, e i are unit vectors).
Theorem 2. Let R be a 2-dimensional spherical random walk on C. Then and α lk = 0 whenever l < 1 2 k. The sum here is over a finite collection of l ≥ 0 and k ≥ 0 with the condition −1 ≥ l − k ≥ −K.
Proof. The complex representation (20) is nothing but a restatement of theorem 1: writing x = 1 2 (z +z) and y = 1 2i (z −z) we get a representation in the form but of course if both i = 0 and j = 0 then we can cancel one of them with the |z| k factor. The fact that the end result is real means that α 0jk = α j0k and we get (20). Thus what we need to prove is the fact that α lk = 0 whenever l < k 2 . We change the representation again, writing and now we have to show that β kl = 0 whenever both k and l are negative (of course, β kl = 0 when k + l > 0). Now, the first thing to notice is that we may differentiate (21) formally. If R moves to v i with probability p i then where we used the one-dimensional Taylor expansion for the functions f (z + tv i ) for each i (the constant in the O(·) above is of course simply 1). Denote the nth differential operator in the sum by D n and denote the maximum inside the O(·) by M K . Clearly, Now, ∆a = 0 with (21) and (22) show that This is exactly what we mean when we say that (21) can be differentiated formally. Now to the actual calculation. Notice that D 1 = 0 (because R is balanced) and that D 2 is, up to multiplication with a constant, the continuous Laplacian, because R is spherical. Therefore for some constant γ = 0. Assume that we know β kl for all k + l > −m, and that β kl = 0 if both are negative. We write m n=1 D n β log zz + k+l>−m β kl z kzl = δ kl z kzl we get that δ kl = 0 if both k and l are negative. Moreover, (23) for K > m + 2 gives equations for β kl with k + l = −m, namely β kl = 1 2γkl δ k−1,l−1 . This shows inductively that β kl = 0 when they are both negative, and the theorem is proven.
The case l = 0 in theorem 2 has some potential uses for precise estimates of hitting probabilities, so it seems worth of special mention: Corollary. For every spherical random walk on R 2 , the high-order expansion of the harmonic potential a, considered as a function of a continuous variable, satisfies the superpolynomial estimate Remark. All the above discussion has analogues in dimension > 2, using the notions of spherical harmonics. Spherical harmonics are functions on S d−1 which are eigenvalues of the angular part of the continuous Laplacian (which is the angular momentum squared operator). It turns out that each is a restriction to S d−1 of a polynomial on R d and we shall use this representation. For a general introduction to the topic, see [SW71,sect. 4.2].
• In two dimensions, a factor P/|z| α (where P is a homogeneous polynomial) may be decomposed into a sum of the terms z kzl , k + l = deg P − α.
The high dimension analogue is a unique decomposition into a sum of terms Q/||z|| β where Q is a homogeneous (continuously) harmonic polynomial, deg Q − β = deg P − α. • The analogue of the fact that only z k andz l are harmonic is the fact that a term Q/||z|| β is harmonic if and only if β = 2 deg Q + d − 2. This follows from the following general formula: and from the harmonicity of Q (remember that ∆ C is the continuous Laplacian). (25) also gives the analogue of (24). • The analogue of "k and l are both negative" is "β < 2 deg Q + d − 2", or, equivalently, "∆ n C (Q/||z|| β ) is never zero, for any n". • The analogue of theorem 2 is the following: if R is a spherical random walk and We omit the proof, which is a similar induction.
• The corollary to theorem 2 is also true in higher dimensions. This follows from the fact that the terms Q/||z|| β for different β's are orthogonal in the sphere measure: Theorem 3. The coefficients of the high-order expansion of the harmonic potential a of the regular random walk on Z 2 can be calculated from the differential equations and (6).
Proof. There is almost nothing to prove here. We saw during the proof of the previous theorem that (23) allows to calculate all coefficients except those where k = 0 or l = 0 (i.e. the coefficients of the harmonic summands). There is only one of these (and its conjugate) in every level k + l = −n. Therefore the coefficient of z −n can be determined from the coefficients of z −n−1z , z −n−2z2 , . . . and from (6).
The proof of theorem 3 is also the second algorithm for getting (4).
Proof. Take as an example clause 2. We may assume that the walk is spherical. In two dimensions, the theorem is equivalent to showing that only terms of the form z kzl , 0 ≤ k ≤ 1 3 |l| (or 0 ≤ l ≤ 1 3 |k|) appear. On the one hand, k can only increase (by 1) in the step of taking D −1 2 (see the proof of theorem 2), and this operation also increases l by 1. On the other hand, the step of taking D 4 +D 5 +· · · decreases l by at least 4, and the theorem follows. In higher dimensions the same argument works using the decomposition into spherical harmonics discussed before theorem 3.
We note that clause 1 of theorem 4 is already mentioned in [FU96,U98]. Theorem 4 can be generalized to cases where higher symmetries exist. For example, for the 6-walk on the triangular lattice, D 3 ≡ D 5 ≡ 0 and D 4 = 9 4 ∆ 2 C which allow to get that the term of order |z| −k is P/|z| ν where ν ≤ 5 2 k and deg P ≤ 3 2 k. Theorem 5. The harmonic potential of a reversible random walk contains only terms of even order in even dimension and only terms of odd order in odd dimension.
Proof. Since D 2k+1 ≡ 0 for all k, we only need to handle the terms we cannot calculate, namely the (continuously) harmonic terms. However, a harmonic term of these orders must be P/||z|| α with deg P odd, and since a(z) = a(−z) such terms cannot appear in the expansion of a.

EXPLICIT CONSTANTS
This note is a spin-off from a different work [K] where it was necessary to have an explicit value for the constants in the O(·) in (2). How should one go about to calculate something like that? Generally, one would try to get some explicit constant in the O of an approximation of one order higher, and then use the McCrea-Whipple algorithm to check the first few values.
Which approach would be best to get an explicit constant in the O(·)? The differential approach doesn't seem to give any explicit constants whatsoever. Therefore it was necessary to use the direct approach. A program was written to handle all the algebra and the calculations of the explicit constants in the O(·). 7 After careful selection of parameters, the program got (4) with an error term which is smaller than ≤ 10 25 |z| −12 for any |z| > 1200. The program doing the algebra ran for half an hour on a PC. This may or may not be enough -it depends on the difference between the actual maximum and the asymptotic value. However, this difference is > 0.01 (see below) and 10 25 |z| −12 < 10 −5 |z| −2 for |z| > 1200 so this quantity is indeed negligible.
Interestingly, here the third approach (the one discussed briefly on page 3) shows a slight advantage over the other two. The third approach allows to get an explicit constant in the O(·) in (8) with much less effort -a few trivial programs to find optimal parameters were all that was necessary. Contrariwise to the simplicity of the programs, the estimates are not as good and it was necessary to run the McCrea-Whipple algorithm for |z| < 6000, which takes a few days.
The maximum of |z| 2 a(z) − 2 π log |z| − λ was found at z = 3 and so the constant implicit in the O(·) in (2) is