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:
N = 20I1 = 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/gThis 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.".

