Suppose we have a 2-dimensional simple random walk: we start at $(0,0)$, and at every step, we add a random unit vector in one of the four cardinal directions selected independently and uniformly.

It is well-known that this procedure will with probability $1$ hit every element of $\mathbb{Z}^2$ infinitely often. Thus, it makes sense to ask about the probability that such a walk will hit $(1,0)$ before $(2,0)$.

Running some Monte Carlo simulations, it looks like the walk first lands on $(1,0)$ something like $70\%$ of the time, but I don't have much confidence about the accuracy of these simulations since I cannot actually run them all to completion and have to either throw out the unfinished trials or make a guess as to how they will conclude. Some more precise simulations show that the walk first lands on $(1,0)$ with probability at least $0.607$ and on $(2,0)$ with probability at least $0.153$.

Is there a known exact value for this probability, and how can it be computed in general for any two points on the square lattice? I'm also interested in a general formula for the probability of encountering $n$ points in a given order.

  • 5
    $\begingroup$ What we are looking for is a discrete harmonic function on the square lattice, with boundary conditions $f(1,0)=1$, $f(2,0)=0$, and then evaluate $f(0,0)$. $\endgroup$ Commented Feb 6, 2021 at 21:59
  • 1
    $\begingroup$ @WillM.: I don't follow your symmetry argument. Why should the probability of reaching $(2,0)$ versus $(1,1)$ first be identical? This is certainly not the case two moves out. $\endgroup$ Commented Feb 6, 2021 at 23:32
  • 1
    $\begingroup$ @HagenvonEitzen In addition, it's safe to assume that $f(n,m) \approx 0.5$ when $n$ and $m$ are large. I've done some inefficient numerical calculations to compute $f(0,0)$ recursively, and it appears that the value is larger than $1 / \sqrt{2}$. I would guess that the value lies around $0.723$. $\endgroup$ Commented Feb 7, 2021 at 10:19
  • 7
    $\begingroup$ There is a good chance that it is $2-4/\pi$ (it agrees with Jeroen's calculation), which is the effective resistance between the origin and (2,0). I calculated what the equivalent 3-node circuit should be for the effective resistances $(1/2,1/2, 2-4/\pi)$, and the desired probability came out to be exactly that. You can read the effective resistance values here: mathpages.com/home/kmath668/kmath668.htm. $\endgroup$
    – E-A
    Commented May 24, 2021 at 3:51
  • 3
    $\begingroup$ For your more general question, you wanted to calculate the probability of hitting a sequence of points in a given chain. In the resistor analogy, to hit point $a$ before $b,c,d$, you’d want the voltage at $a$ to be 1 and the voltage at $b,c,d$ to be 0. The probability of hitting $a$ before $b,c,d$ starting at $x$ is then the voltage at $x$. If you put a point source at each of $a,b,c,d$ you could calculate the implied voltage differences using E-A’s link. By taking linear combinations of them, you can get the $1,0,0,0$ voltages at those points and the voltage at any point $x$. $\endgroup$
    – Eric
    Commented May 27, 2021 at 1:34

4 Answers 4


Update: I originally got the probability is $1-2/\pi$, but there seemed to be mistake taking the wrong order of limits. Taking the right order I got $2-4/\pi$.

Denote $f_n(\mathbf{x})$ the number of ways to end at $\mathbf{x} = (x,y)$ after $n$ steps. Denote $g_n(\mathbf{x})$ the number of ways to end at $\mathbf{x}$ after $n$ steps for the first time. Denote $h_n(\mathbf{x},\mathbf{y})$ the number of ways to end at $\mathbf{x}$ after $n$ steps for the first time having not passed through $\mathbf{y}$.


$$ h_n(\mathbf{x},\mathbf{y}) = g_n(\mathbf{x}) - \sum_{k=0}^{n-1} h_k(\mathbf{y},\mathbf{x}) g_{n-k}(\mathbf{x}-\mathbf{y}) \tag{1} $$ The desired probability is $\lim_{n \rightarrow \infty} \frac{1}{4^n} \sum_{k=0}^n h_k((1,0),(2,0))4^{n-k}$. Similarly

$$ g_n (\mathbf{x}) = f_n(\mathbf{x}) - \sum_{k=0}^{n-1} g_k(\mathbf{x}) f_{n-k}(\mathbf{0}) \tag{2} \\ $$

Thus in principle we can successively solve for $f_n$, $g_n$, $h_n$, sum, and take the limit.

To compute $f_n$:

For $n$ odd, $f_n((2l,0)) = 0, \forall l \in \mathbb{N}$. For $n = 2m$

$$ f_{2m} ((2l,0)) = \sum_{i=0}^{m} {2m \choose 2i} {2i \choose i+l} {2m-2i \choose m-i}. \\ $$ The terms are $0$ for $i < l$. ${2i \choose i+l}$ is the number of ways of stepping horizontally $2i$ times and ending at $l$. ${2m-2i \choose m-i}$ is similar, but for vertical steps. ${2m \choose 2i}$ is the number of ways of alternating between $2i$ horizontal steps and $2m-2i$ vertical steps.

Expanding the terms and simplifying

\begin{align} f_{2m} ((2l,0)) &= {2m \choose m+l} \sum_{i=0}^{m} {m-l \choose i-l} {m+l \choose m-i} \\ &= {2m \choose m+l}^2 \end{align}

by Vandermonde's Convolution. Similarly, for $n$ even, $f_n((1,0)) = 0$, and

\begin{align} f_{2m+1}((1,0)) &= \sum_{i=0}^m {2m+1 \choose 2i+1} {2i+1 \choose i+1} {2m-2i \choose m-i} \\ &= {2m+1 \choose m} \sum_{i=0}^m {m+1 \choose i+1} {m \choose i} \\ &= {2m+1 \choose m}^2 \end{align}

To compute $g_n$:

For even $n, g_n((1,0)) = 0$. For $n = 2m+1, \mathbf{x} = (1,0), (2)$ becomes

$$ \sum_{i=0}^m g_{2i+1}((1,0)) {2m - 2i \choose m-i }^2 = {2m+1 \choose m}^2 \tag{3} $$

Recursively computing a few cases $g_1((1,0)) = 1, g_3((1,0)) = 5, g_5((1,0)) = 44, g_7((1,0)) = 469$. This doesn't appear in OEIS and I don't recognize a pattern.

For odd $n$, $g_n((2,0)) = 0$. For $n=2m$, $(2)$ gives $$ \sum_{i=0}^m g_{2i}((2,0)) {2m - 2i \choose m-i }^2 = {2m \choose m+1}^2 \tag{4} $$

$(3)$ and $(4)$ are convolutions and can be solved using generating functions. Let

\begin{align} A(z) &= \sum_{k \geq 0} {2k \choose k}^2 z^{2k} \\ B(z) &= \sum_{k \geq 0} {2k+1 \choose k}^2 z^{2k+1} \\ C(z) &= \sum_{k \geq 0} {2k \choose k+1}^2 z^{2k} \\ G_{(1,0)}(z) &= \sum_{k \geq 0} g_{2k+1}((1,0))z^{2k+1} \\ G_{(2,0)}(z) &= \sum_{k \geq 0} g_{2k}((2,0))z^{2k} \\ H_{(\mathbf{x},\mathbf{y})}(z) &= \sum_{k \geq 0} h_k(\mathbf{x},\mathbf{y}) z^k \end{align}

The following have the same coefficients by $(3)$ therefore are equal

$$ G_{(1,0)}(z)A(z) = B(z) \tag{5} $$

Similarly by $(4)$ $$ G_{(2,0)}(z)A(z) = C(z) \tag{6} $$

Comment: It would be nice to find $h_k$ in its closed form, but it seems not easy. The complete elliptical integral of the first kind has series $$ K(t) = \frac{\pi}{2} \sum_{k=0}^\infty \left(\frac{(2k)!}{2^{2k}(k!)^2}\right)^2 t^{2k} $$ hence $A(z) = \frac{2}{\pi} K(4 z)$.

For $B(z)$, substituting ${2k+1 \choose k} = \frac{1}{2} {2(k+1) \choose k+1}$ gives $1 + 4zB(z) = A(z)$. By $(5)$

$$ 4zG_{(1,0)}(z) = -\frac{1}{A(z)} + 1 $$

By explicitly differentiating the RHS and setting $z=0$ I verified $g_1((1,0)) = 1, g_3((1,0)) = 5, g_5((1,0)) = 44$. So to find $g_k$ it seems we need the power series for $1/K(t)$.


We're almost done!

Multiplying $(1)$ by $z^n$ and summing over $n$ $$ H_{((1,0),(2,0))}(z) = G_{(1,0)}(z) - H_{((2,0),(1,0))}(z) G_{(1,0)}(z). $$ since $g_k(\mathbf{x}-\mathbf{y}) = g_k(\mathbf{y}-\mathbf{x})$ by symmetry (for even $n$ the equation says $0 = 0$). Similarly $$ H_{((2,0),(1,0))} = G_{(2,0)}(z) - H_{((1,0),(2,0))}(z) G_{(1,0)}(z) $$

Combining $$ H_{((1,0),(2,0))}(z) = G_{(1,0)}(z) - \left(G_{(2,0)}(z) - H_{((1,0),(2,0))}(z) G_{(1,0)}(z)\right) G_{(1,0)}(z) $$ So \begin{align} H_{(1,0),(2,0)}(z) &= G_{(1,0)}(z) \frac{1 - G_{(2,0)}(z)}{1 - G_{(1,0)}(z)^2} \\ &= B(z) \frac{A(z) - C(z)}{A(z)^2 - B(z)^2} \\ &= \frac{A(z) - C(z)}{2(A(z) - B(z))} - \frac{A(z) - C(z)}{2(A(z)+B(z))} \\ \end{align}

Update: Differentiating the RHS a few times gives $H_{((1,0),(2,0))} = z+5z^3 + \cdots$ so this seems correct. The probability is

$$ \lim_{n \rightarrow \infty} \frac{1}{4^n} \sum_{k=0}^n h_k((1,0),(2,0)) 4^{n-k} = H_{((1,0),(2,0))}(1/4) $$

Originally I took the wrong order of limits and got $1-\pi/2$. We have to be a bit careful since $A(1/4) = \frac{2}{\pi} K(1) = + \infty$. To compute this we need to use

\begin{align} H_{((1,0),(2,0))}(1/4) &= \lim_{z \rightarrow 1/4^{+}} H_{((1,0),(2,0))}(z) \\ &= \lim_{z \rightarrow 1/4^{+}} \left(\frac{A(z) - C(z)}{2(A(z) - B(z))} - \frac{A(z) - C(z)}{2(A(z)+B(z))}\right) \end{align} since $A(z), B(z), C(z)$ are convergent and continuous for $0\leq z<1/4$.

Inputing these into mathematica with Limit[(1 - 1/(4 x)) (2 EllipticK[16 x^2])/[Pi] + 1/(4 x) , x -> 1/4] and Limit[(2 EllipticK[16 x^2])/[Pi] - x^2 HypergeometricPFQ[{3/2, 3/2, 2, 2}, {1, 3, 3}, 16 x^2], x -> 1/4] gives $\lim_{z \rightarrow 1/4^+} (A(z) - B(z)) = 1, \lim_{z \rightarrow 1/4^+} (A(z) - C(z)) = 4 - 8/\pi$ respectively.

I don't know how Mathematica put the second in a closed form (the first I think I could do since its in terms of $K$), but since $\lim_{z \rightarrow 1/4^+} (A(z) + B(z)) = +\infty$ we get

$$ H_{((1,0),(2,0))}(1/4) = 2 - \frac{4}{\pi} $$

  • 1
    $\begingroup$ The first series evaluates to $2-4/\pi$; splitting into even-odd indices and using generating functions, the other two evaluate to $\frac{2}{\sqrt{3}}\mp \frac{2}{3} \left(2 \sqrt{3}-3\right)$, which means the last term that is the difference of their reciprocals is $1/2$. So, unless I've made an arithmetic error, your answer doesn't agree with the comments by a factor of $4$, unless the leading $1/2$ is meant to be a $2$, or there's some other error (or, plausibly, the original comment is off). Tl;dr : computations are mostly right, I think. $\endgroup$
    – Integrand
    Commented May 28, 2021 at 22:27
  • $\begingroup$ I did more calculation. A way to know for sure is to calculate $\lim_{z \rightarrow 1/4} \lim_{n \rightarrow \infty} A(z) - C(z)$, where they have been truncated to $n$ terms, instead of $\lim_{n \rightarrow \infty} \lim_{z \rightarrow 1/4} A(z) - C(z)$. This ensures convergence. $\endgroup$ Commented May 29, 2021 at 14:48
  • 1
    $\begingroup$ This looks really good but numerical simulations suggest the answer is closer to $.7$ than $.3$. I'm not sure where, but I strongly suspect a factor of $2$ was dropped somewhere in your work. $\endgroup$
    – QC_QAOA
    Commented May 29, 2021 at 17:48
  • $\begingroup$ Seconding all the comments about a missing factor of $2$ here: you can just directly add up the possible paths of length $\le5$ and see that the probability must be at least $\frac{189}{512}\approx 0.369$. $\endgroup$ Commented May 29, 2021 at 18:15
  • $\begingroup$ @tinklehoy: I agree with your explicit computations of $g_k((1,0))$, so I think everything up to that point is likely correct. Can you use the final expression to obtain coefficients of $h_k((1,0),(2,0))$? I have $1,0,5,0,42,0,429,0,4866,\ldots$ via direct computation. $\endgroup$ Commented May 30, 2021 at 18:40

Here's an answer coming from a potential-theoretic angle. I'll try to set it up for the more general problem of a random walk on $\mathbb{Z}^2$ initiated at the origin and its probability of arriving at one of the points $x_1,\ldots,x_n \in \mathbb{Z}^2$ before visiting any of the points $x_{n+1},\ldots,x_m \in \mathbb{Z}^2$.

  1. The sought-for probability is $p(0,0)$ where $p:\mathbb{Z}^2 \to [0,1]$ is a function that is harmonic on $\mathbb{Z}^2\setminus \{x_1,\ldots,x_m\}$ and evaluates to $p(\{x_1,\ldots,x_n\})=\{1\}$, $p(\{x_{n+1},\ldots,x_m\})=\{0\}$.

  2. The existence of a function $p$ solving the above described boundary value problem can be inferred from the existence of the random walk initiated at the origin.

  3. Concerning the uniqueness of the solution $p$ of that boundary value problem (and to get a hands-on expression for $p$), consider the 2D Green's function $G:\mathbb{Z}^2 \to \mathbb{R}$: $$G(x):=\left(\frac{1}{2\pi}\right)^2\int_{[-\pi,\pi]^2}dk\,\frac{1-\exp(ik\cdot x)}{4-2\cos(k_1)-2\cos(k_2)},$$ (which solves $(\Delta G)(.)= \chi_{\{(0,0)\}}(.)$, where we remind that the Laplacian $\Delta$ acts like $$(\Delta f)(x)=f(x+(1,0))+f(x+(0,1))+f(x-(1,0))+f(x-(0,1))-4f(x)\qquad).$$

  4. One can then check that the function $h(.):=p(.)-\sum_{k=1}^m (\Delta p)(x_k) G(.-x_k)$ is harmonic. One can see then that $\sum_{k=1}^m (\Delta p)(x_k) =0$ for if this sum was strictly positive or strictly negative, we would have that $h$ would asymptotically diverge to $+\infty$ resp. $-\infty$ and therefore $h$ would violate the maximum principle for harmonic functions. Next, $\sum_{k=1}^m (\Delta p)(x_k) =0$ implies that $h$ is bounded and therefore the Liouville theorem for harmonic functions on $\mathbb{Z}^2$ implies that $h$ is a constant function. So this leads to the conclusion that $p$ must necessarily be a linear combination of the functions $\{G(.-x_k)-G(.-x_1)\}_{2\leq k \leq m}\cup \{\chi_{\mathbb{Z}^2}\}$. The coefficients of this linear combination can be determined through the system of linear equations resulting from the evaluations $p(\{x_1,\ldots,x_n\})=\{1\}$, $p(\{x_{n+1},\ldots,x_m\})=\{0\}.$ and by our existence result of point 2 and the just-concluded uniqueness result, this system of linear equations has a solution.

  5. For the specific problem the OP asked about, it is clear then that $$p(.)=\frac{1}{2}\left(1+\frac{G(.-(1,0))-G(.-(2,0))}{G((0,0))-G((-1,0))}\right).$$ Maybe I will come back later to carry out the integrals in this expression.

  • $\begingroup$ (+1) Great work! Just for curious, how do you get $\sum_{k=1}^{m} \Delta p(x_k)=0$? I was able to establish this in OP's case using the symmetry in OP's setting, but I have no idea how I can show this for a general configuration of points $x_k$'s. $\endgroup$ Commented Jan 30 at 16:32
  • $\begingroup$ @SangchulLee A detailed explanation could be too long for a comment, but the knowledge that $p$ is bounded and that we're in 2 dimensions is key. The next tool behind this statement is a discrete analog of Stokes' theorem (used over a bounded but large domain) and the notion that the flux (surface) integral may reveal something about how fast $p(x)$ might be growing as $x$ tends toward infinity. (Hint: if the sum of the sources at the $x_k$ doesn't vanish, one should infer that $p$ must grow logarithmically) $\endgroup$
    – 5th decile
    Commented Jan 30 at 17:09
  • $\begingroup$ @SangchulLee Actually, I was a bit too fast and upon closer investigation I had trouble myself making that explanation explicit. I've modified the answer. $\endgroup$
    – 5th decile
    Commented Jan 30 at 20:52
  • $\begingroup$ Thanks, the last hint really worked for me! $\endgroup$ Commented Jan 31 at 5:35

Not an answer but actually trying this I got a probability of $0.73$. The code is Mathematica Code

An explanation is as follows:

$1)$ Start at $t=0$ in the complex plain

$2)$ Add one of the following randomly $(1,-1,i,-i)$ to $t$

$3)$ Check if $t=1$ or $t=2$ (if so, stop and record this result)

$4)$ Check if the magnitude of the resulting number is greater than $100$

$5)$ If not, got back to step $2)$. If it is, then choose randomly from $t=1$ or $t=2$, record the result, and stop

I do this $10,000$ times and record the result each time. The key assumption I made was that for $|t|>100$ the probability of getting $t=1$ or $t=2$ first was about $50\%$ each.

Update: I re-ran my code with $\lambda=10^6$ and got a final probability of $P=0.726992$. This agrees heavily with the conjectured answer of $2-\frac{4}{\pi}$ as



This is a community answer expanding @Vergilius's potential-theoretic approach. Feel free to improve it to your liking!

Potential Kernel. Let $\mathfrak{a} : \mathbb{Z}^2 \to \mathbb{R}$ denote the potential kernel given by

$$ \mathfrak{a}(x) = \frac{2}{4\pi^2} \int_{\mathbb{T}^2} \frac{1 - \cos(x\cdot k)}{2 - \cos(k_1) - \cos(k_2)} \, \mathrm{d}k, $$

where $\mathbb{T} = \mathbb{R}/2\pi\mathbb{Z}$. This satisfies the asymptotic formula

$$ \mathfrak{a}(x) = g\log\|x\| + \mathcal{O}(1) $$

as $\|x\| \to \infty$ with $g = \frac{2}{\pi}$ (see [Stö49]). Also, it is easy to verify that this solves the Poisson equation

$$ \Delta \mathfrak{a} = 4 \delta_0, $$

where $\Delta f(x) = \sum_{y \in \mathcal{N}(x)} (f(y) - f(x)) $ is the Laplacian and $\mathcal{N}(x)$ is the set of nearest neighbors of $x$ in $\mathbb{Z}^2$, and $\delta_0(x) = \mathbf{1}[x = 0]$ is the Kronecker delta.

Hitting Probabilities. Let $(S_n)$ denote the simple random walk (SRW) on $\mathbb{Z}^2$. Also, let $\mathbf{P}^x(\cdot)$ be the law of $(S_n)$ started at $x$. If

$$ \tau_W = \inf\{ n \geq 0 : S_n \in W \} $$

denotes the first hitting time of the set $W \subseteq \mathbb{Z}^2$, then the recurrence of $(S_n)$ tells that $\tau_W < \infty$ holds $\mathbf{P}^x$-a.s. for any non-empty subset $W$.

Now let us identify $\mathbb{R}^2$ with $\mathbb{C}$ for notational simplicity. Then with $W = \{1, 2\}$, we are interested in the function

$$ p(x) = \mathbf{P}^x( S_{\tau_W} = 1 ), $$

and in particular, its value $p(0)$ at the origin. To this end, we list key properties of $p$ to be used:

  • We have $\Delta p = 0$ on $\mathbb{Z}^2 \setminus W$. This is easily proved from the Markov property.

  • It satisfies $p(1) = 1$ and $p(2) = 0$.

  • By symmetry, we have $p(3-x) = 1 - p(x)$.

  • Using the above identity, we can easily verify that $\Delta p(1) + \Delta p(2) = 0$ holds.

Combining Together. Let $\phi : \mathbb{Z}^2 \to \mathbb{R}$ by

$$ \phi(x) = p(x) - c[\mathfrak{a}(x-1) - \mathfrak{a}(x-2)], $$

where $c = \frac{1}{4}\Delta p(1) = -\frac{1}{4}\Delta p(2)$. Then it is clear that $\Delta \phi = 0$ identically on $\mathbb{Z}^2$. Moreover, the potential kernel asymptotics guarantees that $\phi$ is bounded. So by the Liouville's theorem, $\phi$ is constant. Using the identity $p(3-x) + p(x) = 1$, we find that this constant is precisely $\frac{1}{2}$. Hence,

$$ p(x) = \frac{1}{2} + c [\mathfrak{a}(x-1) - \mathfrak{a}(x-2)]. $$

The value of constant $c$ can be determined by plugging $x = 1$ and $x = 2$, yielding

$$ p(x) = \frac{1}{2}\left( 1 - \frac{\mathfrak{a}(x-1) - \mathfrak{a}(x-2)}{\mathfrak{a}(1)} \right). $$

(Here, we utilized the symmetry $\mathfrak{a}(-x) = \mathfrak{a}(x)$ and $\mathfrak{a}(0) = 0$.) Then plugging $x = 0$ into the above identity, the desired probability is

$$ p(0) = \frac{\mathfrak{a}(2)}{2\mathfrak{a}(1)}. \tag{1} $$

So it suffices to compute $\mathfrak{a}(1)$ and $\mathfrak{a}(2)$. The first one is easily resolved using symmetry:

$$ \mathfrak{a}(1) = \frac{2}{4\pi^2} \int_{\mathbb{T}^2} \frac{1 - \cos(k_1)}{2 - \cos(k_1) - \cos(k_2)} \, \mathrm{d}k = \frac{1}{4\pi^2} \int_{\mathbb{T}^2} \mathrm{d}k = 1. $$

Moreover, invoking the formula $\frac{1}{2\pi} \int_{0}^{2\pi} \frac{\mathrm{d}\theta}{a-\cos\theta} = \frac{1}{\sqrt{a^2 - 1}}$, $a > 1$, we get

\begin{align*} \mathfrak{a}(2) &= \frac{2}{4\pi^2} \int_{\mathbb{T}^2} \frac{1 - \cos(2k_1)}{2 - \cos(k_1) - \cos(k_2)} \, \mathrm{d}k \\ &= \frac{1}{\pi} \int_{\mathbb{T}} \frac{1 - \cos(2k_1)}{\sqrt{(2 - \cos(k_1))^2 - 1}} \, \mathrm{d}k_1 \\ &= \frac{4}{\pi} \int_{0}^{\pi} \frac{\sin^2(k_1)}{\sqrt{(1 - \cos(k_1))(3 - \cos(k_1))}} \, \mathrm{d}k_1 \\ &= \frac{4}{\pi} \int_{-1}^{1} \frac{\sqrt{1-t^2}}{\sqrt{(1 - t)(3 - t)}} \, \mathrm{d}t \tag{$t=\cos k_1$} \\ &= \frac{4}{\pi} \int_{-1}^{1} \sqrt{\frac{1+t}{3-t}} \, \mathrm{d}t \\ &= \frac{4}{\pi} \int_{0}^{1} \frac{4\sqrt{u}}{(1+u)^2} \, \mathrm{d}u \tag{$u = \frac{1+t}{3-t}$} \\ &= \frac{4}{\pi} (\pi - 2). \end{align*}

Plugging these into $\text{(1)}$, it therefore follows that

$$ p(0) = 2 - \frac{4}{\pi}. $$

[Stö49] Alfred Stöhr. Über einige lineare partielle differenzengleichungen mit konstanten koeffizienten. Mathematische Nachrichten, 3(6):330–357, 1949.


You must log in to answer this question.

Not the answer you're looking for? Browse other questions tagged .