Let be the probability that a random walk on a d-D lattice returns to the origin. In 1921, Pólya proved that $p(1)=p(2)=1$
but $p(d)<1$ for $d>2$.
<http://mathworld.wolfram.com/PolyasRandomWalkConstants.html>


I wonder what we can say about the probability for $d \to \infty$
In other words, if there is a closed formula or approximation which presents the limit of the probability for very big $d$?
Does the limit exist?


Thank you in advance for any comments or approach to investigate the question.



Pólya's Random Walk Constants at infinity
The table in that Mathworld page suggests that $p(d) \rightarrow 0$ as $d \rightarrow \infty$.
That page also gives a formula for $p(d)$ in terms of a definite integral:
$$
p(d) = 1 - \left[ \int\_0^\infty I\_0(t/d)^d e^{-t} dt \right]^{-1},
$$
where $I\_0$ is a "modified Bessel function" with power series
$$
I\_0(x) = \sum\_{n=0}^\infty \frac{(x/2)^{2n}}{n!^2}
= 1 + \frac{x^2}{2^2} + \frac{x^4}{8^2} + \frac{x^6}{48^2} + \cdots .
$$
[For large $x$ it is known that $I\_0(x) \sim (2\pi x)^{-1/2} e^x$, so the integrand decays as a multiple of $x^{d/2}$ for $x \rightarrow \infty$, and the integral is finite **iff** $d>2$.] Substituting the power series for $I\_0(x)$ into the integral, and expanding termwise via $\int\_0^\infty t^m e^{-t} dt = m!$, yields
$$
1 + \frac1{2d} + \frac3{4d^2} + \frac3{2d^3} + \frac{15}{4d^4} + \frac{355}{32d^5} + \cdots
$$
(the coefficients in powers of $1/2d$ are [OEIS sequence A105227](https://oeis.org/A105227)). We then compute the asymptotic series
$$
p(d) \sim \frac1{2d} + \frac1{2d^2} + \frac7{8d^3} + \frac{35}{16d^4} + \frac{215}{32d^5} + \cdots,
$$
which seems consistent with the Mathworld table. The coefficients in powers of $1/2d$ are [OEIS sequence A043546](https://oeis.org/A043546); as often happens, finding a few terms makes it easier to hunt down some of the literature. [Added later: so does posting on Mathoverflow; a comment to this answer by Folkmar Bornemann gives a reference dating back 55+ years



> 
> Montroll, Elliot W: Random walks in multidimensional spaces, especially on periodic lattices, *J. Soc. Indust. Appl. Math.* **4** (1956), 241–260 (MR0088110)
> 
> 
> 


— thanks! Here's some **gp** code for this power series in $w = 1/2d$ and its coefficients, quite similar to Flajolet's Maple code reproduced in the [OEIS entry](https://oeis.org/A043546):



```
N = 20
I1 = sum(n=0,N,x^n/n!^2,O(x^(N+1)));
Iw = subst(I1,x,w^2*x)^(1/(2*w));
g = sum(n=0,N,(2*n)!*polcoeff(Iw,n,x)) + O(w^(N+1));
p = 1 - 1/g
vector(N,n,polcoeff(p,n))

```

This returns



```
[1, 2, 7, 35, 215, 1501, 11354, 88978, 675569, 4175664, 1725333, -687775083, -19848956619, -438027976068, -8715988203509, -161989586455204, -2784493824166078, -41530410660307610, -406672888265416456, 4420077014249902362]

```

and **gp** readily computes it for $N$ as large as 50, and with some more effort even for $N=100$.]


The form of the asymptotic series can be explained as follows: for each $k=1,2,3,\ldots$, the probability of return to the origin in $2k$ steps is $O(d^{-k})$ as $d \rightarrow\infty$; so the probability of return by the $2k$-step gives $p(d)$ to within $O(d^{-(k+1)})$, and this estimate is a polynomial in $1/2d$ with integer coefficients. For example, the $k=1$ probability is $1/2d$ exactly; for $k=2$, add $2!/(2d)^2 - O(d^{-3})$; "etc.".


