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.

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