# What is the probability that a random walk on $\mathbb{Z}^2$ will hit $(1,0)$ before $(2,0)$?

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.

• 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)$. Commented Feb 6, 2021 at 21:59
• @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. Commented Feb 6, 2021 at 23:32
• @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$. Commented Feb 7, 2021 at 10:19
• 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.
– E-A
Commented May 24, 2021 at 3:51
• 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$.
– Eric
Commented May 27, 2021 at 1:34

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}$$.

Note

$$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)$$.

Probability:

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}$$

• 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. Commented May 28, 2021 at 22:27
• 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. Commented May 29, 2021 at 14:48
• 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. Commented May 29, 2021 at 17:48
• 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$. Commented May 29, 2021 at 18:15
• @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. 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.

• (+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. Commented Jan 30 at 16:32
• @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) Commented Jan 30 at 17:09
• @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. Commented Jan 30 at 20:52
• Thanks, the last hint really worked for me! Commented Jan 31 at 5:35

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

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

$$\left|2-\frac{4}{\pi}-\frac{726992}{10^6}\right|<10^{-3}$$

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.