Return to Answer

As I mentioned in my comment - what you are suggesting implies that the probability of being at (3,4) is the same as being at (5,0) for all large n. That seems unlikely, and would guess that $C_n=5$ for n large.

The answer to your modified question is yes! $\tilde C_n$ tends to infinity as n goes to infinity. (Phew! It took me a couple of revisions to prove this, but hopefully the calculations below are now correct).

In fact, $\tilde C_n\ge c\sqrt{n}$ for some positive constants c. I think that you can also show that $\tilde C_n\le C\sqrt{n}$ for some other constant C but I'm not completely sure yet, although it should follow from a closer examination of my expression below for $p_n(x)$.

You can derive an asymptotic expansion for $p_n(x)$ in 1/n. Evaluating this to second order is enough to answer your question. After n steps the distribution of the particle will be approximately normal with variance 2n/5 in both dimensions, so we expect to get $p_n(x)=\frac{5}{4\pi n}e^{-\frac{5}{4n}\vert x\vert^2}$ to leading order.

The idea is to note that you are repeatedly applying a linear operator, $$ p_{n+1}=Lp_n,\ Lp(x) \equiv (p(x)+p(x-e_1)+p(x+e_1)+p(x-e_2)+p(x+e_2))/5 $$ where $e_1=(1,0)$, $e_2=(0,1)$. In finite dimensional spaces, you would solve this by decomposing $p_0$ into a sum of eigenvectors and for large n, the dominant term of $L^np_0$ will be that corresponding to the largest eigenvalue. In this case, the infinite dimensional operator L has a continuous spectrum, and is diagonalized by a Fourier transform. $$ p_0(x)=1_{\lbrace x=0\rbrace}=\int_{-[\frac12,\frac12]^2}e^{2\pi ix\cdot u}\\,du. $$ Noting that $e^{2\pi ix\cdot u}$ (as a function of x) is an eigenvector of L, $$ Le^{2\pi ix\cdot u}=\left(\frac15+\frac25\cos(2\pi u_1)+\frac25\cos(2\pi u_2)\right)e^{2\pi ix\cdot u} $$ gives the following for $p_n$, $$ p_n(x)=L^np_0(x)=\int_{[-\frac12,\frac12]^2}\left(\frac15+\frac25\cos(2\pi u_1)+\frac25\cos(2\pi u_2)\right)^ne^{2\pi ix\cdot u}\\,du. $$ The term inside the parentheses is less than 1 in absolute value everywhere away from the origin, so looks like a Dirac delta when raised to a high power n. Using a Taylor expansion to second order, $$ \left(\frac15+\frac25\cos(2\pi u_1)+\frac25\cos(2\pi u_2)\right)^n =e^{-\frac45\pi^2n\vert u\vert^2}\left(1+\frac{8\pi^4n}{75}(7\vert u\vert^4-20u_1^2u_2^2)+O(n\vert u\vert^6)\right). $$ This expansion is valid over any domain on which $n\vert u\vert^6$ is bounded. Say, $\vert u\vert\le n^{-1/6}$. Outside of this domain, the integrand above is bounded by $e^{-cn(n^{-1/6})^2}=e^{-cn^{2/3}}$ for a constant c, which is much smaller than O(1/n^3) and can be neglected. Then, $$ p_n(x)=\int_{\mathbb{R}^2}\left(1+\frac{8\pi^4n}{75}(7\vert u\vert^4-20u_1^2u_2^2)+O(n\vert u\vert^6)\right)e^{-\frac45\pi^2n\vert u\vert^2+2\pi ix\cdot u}\\,du. $$ Here I not only substituted in the second order approximation to the integrand, but also extended the range of integration out to infinity. This is fine, because it can be shown that the value of this integral over $\vert u\vert\ge n^{-1/6}$ has size of the order of no more than $e^{-cn^{2/3}}$, so vanishes much faster than $O(1/n^3)$. Substituting in $v=\sqrt{\frac{8n}{5}}\pi u$ also shows that the $O(nu^6)$ term in the integrand vanishes at rate $1/n^3$, giving the following. $$ p_n(x)=\frac{5}{8\pi^2n}\int_{\mathbb{R}^2}\left(1+\frac{1}{24n}(7\vert v\vert^4-20v_1^2v_2^2)\right)e^{-\frac12\vert v\vert^2+i\sqrt{\frac{5}{2n}}x\cdot v}\\,dv+O(n^{-3}). $$ This integral can be computed, $$ p_n(x)=\frac{5}{4\pi n}e^{-\frac{5}{4n}\vert x\vert^2}\left(1+\frac{1}{24n}\left(36-\frac{90}{n}\vert x\vert^2+\frac{175}{4n^2}\vert x\vert^4-\frac{125}{n^2}x_1^2x_2^2\right)\right)+O(n^{-3}). $$ This is a bit messy, but the exact coefficients are not too important. What matters is the general form of the expression. The leading order term also agrees with the guess above based on it being approximately normal. Also, for any fixed $\vert x\vert \lt\vert y\vert$, the leading order term in $p_n(x)-p_n(y)$ will dominate for large n, giving $p_n(x)\gt p_n(y)$. So, $\tilde C_n\to\infty$.

Consider $\vert x\vert\le c\sqrt{n}$ for some $c\le1$. Then, $$ p_n(x)=\frac{5}{4\pi n}e^{-\frac{5}{4n}\vert x\vert^2}\left(1+\frac{3}{2n}\right)+O(c^2n^{-2}). $$ If $\vert x\vert\lt\vert y\vert\le c\sqrt{n}$ then $\vert y\vert^2-\vert x\vert^2\ge 1$ (as it is a nonzero integer) $$ \begin{align} p_n(x)-p_n(y)&=\frac{5}{4\pi n}\left(1+\frac{3}{2n}\right)e^{-\frac{5}{4n}\vert x\vert^2}\left(1-e^{-\frac{5}{4n}(\vert y\vert^2-\vert x\vert^2)}\right)+O(c^2n^{-2})\\\\ &\ge\frac{5}{4\pi n}e^{-\frac{5}{4n}\vert x\vert^2}(1-e^{-\frac{5}{4n}})+O(c^2n^{-2})\\\\ &=\frac{25}{16\pi n^2}e^{-\frac{5}{4n}\vert x\vert^2}\left(1+O(c^2)\right). \end{align} $$ As long as c is chosen small enough that the $O(c^2)$ term is always greater than -1, this expression will be positive. So $p_n(x)\gt p_n(y)$ for all $\vert x\vert\lt\vert y\vert\le c\sqrt{n}$, giving $\tilde C_n\ge c\sqrt{n}$.

As I mentioned in my comment - what you are suggesting implies that the probability of being at (3,4) is the same as being at (5,0) for all large $n$. That seems unlikely, and would guess that $C_n=5$ for $n$ large.

The answer to your modified question is yes! $\tilde C_n$ tends to infinity as $n$ goes to infinity. (Phew! It took me a couple of revisions to prove this, but hopefully the calculations below are now correct).

In fact, $\tilde C_n\ge c\sqrt{n}$ for some positive constants c. I think that you can also show that $\tilde C_n\le C\sqrt{n}$ for some other constant $C$ but I'm not completely sure yet, although it should follow from a closer examination of my expression below for $p_n(x)$.

You can derive an asymptotic expansion for $p_n(x)$ in 1/n. Evaluating this to second order is enough to answer your question. After $n$ steps the distribution of the particle will be approximately normal with variance 2n/5 in both dimensions, so we expect to get $p_n(x)=\frac{5}{4\pi n}e^{-\frac{5}{4n}\vert x\vert^2}$ to leading order.

The idea is to note that you are repeatedly applying a linear operator, $$ p_{n+1}=Lp_n,\ Lp(x) \equiv (p(x)+p(x-e_1)+p(x+e_1)+p(x-e_2)+p(x+e_2))/5 $$ where $e_1=(1,0)$, $e_2=(0,1)$. In finite dimensional spaces, you would solve this by decomposing $p_0$ into a sum of eigenvectors and for large n, the dominant term of $L^np_0$ will be that corresponding to the largest eigenvalue. In this case, the infinite dimensional operator $L$ has a continuous spectrum, and is diagonalized by a Fourier transform. $$ p_0(x)=1_{\lbrace x=0\rbrace}=\int_{-[\frac12,\frac12]^2}e^{2\pi ix\cdot u}\,du. $$ Noting that $e^{2\pi ix\cdot u}$ (as a function of $x$) is an eigenvector of $L$, $$ Le^{2\pi ix\cdot u}=\left(\frac15+\frac25\cos(2\pi u_1)+\frac25\cos(2\pi u_2)\right)e^{2\pi ix\cdot u} $$ gives the following for $p_n$, $$ p_n(x)=L^np_0(x)=\int_{[-\frac12,\frac12]^2}\left(\frac15+\frac25\cos(2\pi u_1)+\frac25\cos(2\pi u_2)\right)^ne^{2\pi ix\cdot u}\,du. $$ The term inside the parentheses is less than 1 in absolute value everywhere away from the origin, so looks like a Dirac delta when raised to a high power n. Using a Taylor expansion to second order, $$ \left(\frac15+\frac25\cos(2\pi u_1)+\frac25\cos(2\pi u_2)\right)^n =e^{-\frac45\pi^2n\vert u\vert^2}\left(1+\frac{8\pi^4n}{75}(7\vert u\vert^4-20u_1^2u_2^2)+O(n\vert u\vert^6)\right). $$ This expansion is valid over any domain on which $n\vert u\vert^6$ is bounded. Say, $\vert u\vert\le n^{-1/6}$. Outside of this domain, the integrand above is bounded by $e^{-cn(n^{-1/6})^2}=e^{-cn^{2/3}}$ for a constant c, which is much smaller than O(1/n^3) and can be neglected. Then, $$ p_n(x)=\int_{\mathbb{R}^2}\left(1+\frac{8\pi^4n}{75}(7\vert u\vert^4-20u_1^2u_2^2)+O(n\vert u\vert^6)\right)e^{-\frac45\pi^2n\vert u\vert^2+2\pi ix\cdot u}\,du. $$ Here I not only substituted in the second order approximation to the integrand, but also extended the range of integration out to infinity. This is fine, because it can be shown that the value of this integral over $\vert u\vert\ge n^{-1/6}$ has size of the order of no more than $e^{-cn^{2/3}}$, so vanishes much faster than $O(1/n^3)$. Substituting in $v=\sqrt{\frac{8n}{5}}\pi u$ also shows that the $O(nu^6)$ term in the integrand vanishes at rate $1/n^3$, giving the following. $$ p_n(x)=\frac{5}{8\pi^2n}\int_{\mathbb{R}^2}\left(1+\frac{1}{24n}(7\vert v\vert^4-20v_1^2v_2^2)\right)e^{-\frac12\vert v\vert^2+i\sqrt{\frac{5}{2n}}x\cdot v}\,dv+O(n^{-3}). $$ This integral can be computed, $$ p_n(x)=\frac{5}{4\pi n}e^{-\frac{5}{4n}\vert x\vert^2}\left(1+\frac{1}{24n}\left(36-\frac{90}{n}\vert x\vert^2+\frac{175}{4n^2}\vert x\vert^4-\frac{125}{n^2}x_1^2x_2^2\right)\right)+O(n^{-3}). $$ This is a bit messy, but the exact coefficients are not too important. What matters is the general form of the expression. The leading order term also agrees with the guess above based on it being approximately normal. Also, for any fixed $\vert x\vert \lt\vert y\vert$, the leading order term in $p_n(x)-p_n(y)$ will dominate for large n, giving $p_n(x)\gt p_n(y)$. So, $\tilde C_n\to\infty$.

Consider $\vert x\vert\le c\sqrt{n}$ for some $c\le1$. Then, $$ p_n(x)=\frac{5}{4\pi n}e^{-\frac{5}{4n}\vert x\vert^2}\left(1+\frac{3}{2n}\right)+O(c^2n^{-2}). $$ If $\vert x\vert\lt\vert y\vert\le c\sqrt{n}$ then $\vert y\vert^2-\vert x\vert^2\ge 1$ (as it is a nonzero integer) $$ \begin{align} p_n(x)-p_n(y)&=\frac{5}{4\pi n}\left(1+\frac{3}{2n}\right)e^{-\frac{5}{4n}\vert x\vert^2}\left(1-e^{-\frac{5}{4n}(\vert y\vert^2-\vert x\vert^2)}\right)+O(c^2n^{-2})\\ &\ge\frac{5}{4\pi n}e^{-\frac{5}{4n}\vert x\vert^2}(1-e^{-\frac{5}{4n}})+O(c^2n^{-2})\\ &=\frac{25}{16\pi n^2}e^{-\frac{5}{4n}\vert x\vert^2}\left(1+O(c^2)\right). \end{align} $$ As long as $c$ is chosen small enough that the $O(c^2)$ term is always greater than -1, this expression will be positive. So $p_n(x)\gt p_n(y)$ for all $\vert x\vert\lt\vert y\vert\le c\sqrt{n}$, giving $\tilde C_n\ge c\sqrt{n}$.

The answer to your modified question is yes, $\tilde C_n$ tends to infinity as n goes to infinity. Edit: I'm not so sure now, there was a mistake in my original calculation. What you can say is that there is a positive constant $a$ such that the result holds as long as you restrict to $\vert x\vert^2\lt\vert y\vert^2-a$. Whether or not we can take $a=1$ (which would imply the full result, as $\vert x\vert^2-\vert y\vert^2\in\mathbb{Z}$) requires a more careful consideration of the second order terms in my calculation. In any case, the rate will be of the order of $\sqrt{n}$. That is, $\tilde C_n\ge b\sqrt{n}$ for some constant $b$.

This can be proven by evaluating $p_n(x)$ to leading order in 1/n. (assuming I haven't made any errors below).

The idea is to note that you are repeatedly applying a linear operator, $$ p_{n+1}=Lp_n,\ Lp(x) \equiv (p(x)+p(x-e_1)+p(x+e_1)+p(x-e_2)+p(x+e_2))/5 $$ where $e_1=(1,0)$, $e_2=(0,1)$. In finite dimensional spaces, you would solve this by decomposing $p_0$ into a sum of eigenvectors and for large n, the dominant term of $L^np_0$ will be that corresponding to the largest eigenvalue. In this case, we can diagonalize by applying a Fourier transform $$ p_0(x)=1_{\lbrace x=0\rbrace}=\int_{-[\frac12,\frac12]^2}e^{2\pi ix\cdot u}\\,du. $$ Noting that $e^{2\pi ix\cdot u}$ (as a function of x) is an eigenvector of L, $$ Le^{2\pi ix\cdot u}=\left(\frac15+\frac25\cos(2\pi u_1)+\frac25\cos(2\pi u_2)\right)e^{2\pi ix\cdot u} $$ gives the following for $p_n$, $$ p_n(x)=L^np_0(x)=\int_{[-\frac12,\frac12]^2}\left(\frac15+\frac25\cos(2\pi u_1)+\frac25\cos(2\pi u_2)\right)^ne^{2\pi ix\cdot u}\\,du. $$ The term inside the parentheses is less than 1 in absolute value everywhere away from the origin, so looks like a Dirac delta when raised to a high power n.

The integral can be computed to leading order. After n steps, the standard deviation of the particle's distance from the origin grows of the order of $\sqrt{n}$, so most of its distribution is spread over an area of size the order of n. This means we expect $p_n(x)$ to go to zero at rate 1/n. Let's discard all terms which vanish faster than this.

The range of integration can be replaced by $[-\epsilon,\epsilon]$ for any $0<\epsilon<1/2$, creating an error only of order $e^{-cn}$ (some positive c). On such a range, $$ \frac15+\frac25\cos(2\pi u_1)+\frac25\cos(2\pi u_2)=\exp\left(-\frac{4\pi^2}{5}u^2+O(u^4)\right) $$ Substituting into the integral and changing variables, $$ \begin{align} p_n(x)&=\int_{[-\epsilon,\epsilon]^2}\exp\left(-\frac{4n\pi^2}{5}\vert u\vert^2+O(n\vert u\vert^4)+2\pi ix\cdot u\right)\\,du+O(e^{-cn})\\\\ &=\frac{5}{8n\pi^2}\int_{[-\epsilon'\sqrt{n},\epsilon'\sqrt{n}]^2}\exp\left(-\frac12\vert v\vert^2+O(\vert v\vert^4/n)+i\sqrt{\frac{5}{2n}}x\cdot v\right)\\,dv+O(e^{-cn}) \end{align} $$ Where $\epsilon'=\epsilon\sqrt{8\pi^2/5}$. Rearranging the v^4/n term in the exponential, it is possible to show that it only contributes of the order of 1/n to the integral. Integrals like $\int_{\epsilon\sqrt{n}}^\infty\exp(-x^2/2)\\,dx$ vanish at rate $\exp(-\epsilon^2n/2)/n$, giving $$ \begin{align} p_n(x)&=\frac{5}{8n\pi^2}\int_{\mathbb{R}^2}\exp\left(-\frac12\vert v\vert^2+i\sqrt{\frac{5}{2n}}x\cdot v\right)\\,dv+O(1/n^2)\\\\ &=\frac{5}{4n\pi}\exp\left(-\frac{5\vert x\vert^2}{4n}\right)+O(1/n^2) \end{align} $$ So, if $\vert x\vert\lt\vert y\vert$ then $p_n(x)\gt p_n(y)$ for large n. You can check that if $\vert x\vert\lt\vert y\vert-a$ for some constant a, $$ \begin{align} p_n(x)-p_n(y)&=\frac{5}{4n\pi}e^{-\frac{5\vert x\vert^2}{4n}}(1-e^{\frac{5(\vert x\vert^2-\vert y\vert^2)}{4n}})+O(1/n^2)\\\\ &\ge\frac{5}{4n\pi}e^{-\frac{5\vert x\vert^2}{4n}}(1-e^{-\frac{5a}{4n}})+O(1/n^2)\\\\ &\ge\frac{25}{4n^2\pi}\left(ae^{-\frac{5\vert x\vert^2}{4n}}+O(1)\right) \end{align} $$ To be certain that this is positive, we would need to evaluate the constant term in the final expression and compare to $a$.

The answer to your modified question is yes! $\tilde C_n$ tends to infinity as n goes to infinity. (Phew! It took me a couple of revisions to prove this, but hopefully the calculations below are now correct). 

In fact, $\tilde C_n\ge c\sqrt{n}$ for some positive constants c. I think that you can also show that $\tilde C_n\le C\sqrt{n}$ for some other constant C but I'm not completely sure yet, although it should follow from a closer examination of my expression below for $p_n(x)$.

You can derive an asymptotic expansion for $p_n(x)$ in 1/n. Evaluating this to second order is enough to answer your question. After n steps the distribution of the particle will be approximately normal with variance 2n/5 in both dimensions, so we expect to get $p_n(x)=\frac{5}{4\pi n}e^{-\frac{5}{4n}\vert x\vert^2}$ to leading order.

The idea is to note that you are repeatedly applying a linear operator, $$ p_{n+1}=Lp_n,\ Lp(x) \equiv (p(x)+p(x-e_1)+p(x+e_1)+p(x-e_2)+p(x+e_2))/5 $$ where $e_1=(1,0)$, $e_2=(0,1)$. In finite dimensional spaces, you would solve this by decomposing $p_0$ into a sum of eigenvectors and for large n, the dominant term of $L^np_0$ will be that corresponding to the largest eigenvalue. In this case, the infinite dimensional operator L has a continuous spectrum, and is diagonalized by a Fourier transform. $$ p_0(x)=1_{\lbrace x=0\rbrace}=\int_{-[\frac12,\frac12]^2}e^{2\pi ix\cdot u}\\,du. $$ Noting that $e^{2\pi ix\cdot u}$ (as a function of x) is an eigenvector of L, $$ Le^{2\pi ix\cdot u}=\left(\frac15+\frac25\cos(2\pi u_1)+\frac25\cos(2\pi u_2)\right)e^{2\pi ix\cdot u} $$ gives the following for $p_n$, $$ p_n(x)=L^np_0(x)=\int_{[-\frac12,\frac12]^2}\left(\frac15+\frac25\cos(2\pi u_1)+\frac25\cos(2\pi u_2)\right)^ne^{2\pi ix\cdot u}\\,du. $$ The term inside the parentheses is less than 1 in absolute value everywhere away from the origin, so looks like a Dirac delta when raised to a high power n. Using a Taylor expansion to second order, $$ \left(\frac15+\frac25\cos(2\pi u_1)+\frac25\cos(2\pi u_2)\right)^n =e^{-\frac45\pi^2n\vert u\vert^2}\left(1+\frac{8\pi^4n}{75}(7\vert u\vert^4-20u_1^2u_2^2)+O(n\vert u\vert^6)\right). $$ This expansion is valid over any domain on which $n\vert u\vert^6$ is bounded. Say, $\vert u\vert\le n^{-1/6}$. Outside of this domain, the integrand above is bounded by $e^{-cn(n^{-1/6})^2}=e^{-cn^{2/3}}$ for a constant c, which is much smaller than O(1/n^3) and can be neglected. Then, $$ p_n(x)=\int_{\mathbb{R}^2}\left(1+\frac{8\pi^4n}{75}(7\vert u\vert^4-20u_1^2u_2^2)+O(n\vert u\vert^6)\right)e^{-\frac45\pi^2n\vert u\vert^2+2\pi ix\cdot u}\\,du. $$ Here I not only substituted in the second order approximation to the integrand, but also extended the range of integration out to infinity. This is fine, because it can be shown that the value of this integral over $\vert u\vert\ge n^{-1/6}$ has size of the order of no more than $e^{-cn^{2/3}}$, so vanishes much faster than $O(1/n^3)$. Substituting in $v=\sqrt{\frac{8n}{5}}\pi u$ also shows that the $O(nu^6)$ term in the integrand vanishes at rate $1/n^3$, giving the following. $$ p_n(x)=\frac{5}{8\pi^2n}\int_{\mathbb{R}^2}\left(1+\frac{1}{24n}(7\vert v\vert^4-20v_1^2v_2^2)\right)e^{-\frac12\vert v\vert^2+i\sqrt{\frac{5}{2n}}x\cdot v}\\,dv+O(n^{-3}). $$ This integral can be computed, $$ p_n(x)=\frac{5}{4\pi n}e^{-\frac{5}{4n}\vert x\vert^2}\left(1+\frac{1}{24n}\left(36-\frac{90}{n}\vert x\vert^2+\frac{175}{4n^2}\vert x\vert^4-\frac{125}{n^2}x_1^2x_2^2\right)\right)+O(n^{-3}). $$ This is a bit messy, but the exact coefficients are not too important. What matters is the general form of the expression. The leading order term also agrees with the guess above based on it being approximately normal. Also, for any fixed $\vert x\vert \lt\vert y\vert$, the leading order term in $p_n(x)-p_n(y)$ will dominate for large n, giving $p_n(x)\gt p_n(y)$. So, $\tilde C_n\to\infty$.

Consider $\vert x\vert\le c\sqrt{n}$ for some $c\le1$. Then, $$ p_n(x)=\frac{5}{4\pi n}e^{-\frac{5}{4n}\vert x\vert^2}\left(1+\frac{3}{2n}\right)+O(c^2n^{-2}). $$ If $\vert x\vert\lt\vert y\vert\le c\sqrt{n}$ then $\vert y\vert^2-\vert x\vert^2\ge 1$ (as it is a nonzero integer) $$ \begin{align} p_n(x)-p_n(y)&=\frac{5}{4\pi n}\left(1+\frac{3}{2n}\right)e^{-\frac{5}{4n}\vert x\vert^2}\left(1-e^{-\frac{5}{4n}(\vert y\vert^2-\vert x\vert^2)}\right)+O(c^2n^{-2})\\\\ &\ge\frac{5}{4\pi n}e^{-\frac{5}{4n}\vert x\vert^2}(1-e^{-\frac{5}{4n}})+O(c^2n^{-2})\\\\ &=\frac{25}{16\pi n^2}e^{-\frac{5}{4n}\vert x\vert^2}\left(1+O(c^2)\right). \end{align} $$ As long as c is chosen small enough that the $O(c^2)$ term is always greater than -1, this expression will be positive. So $p_n(x)\gt p_n(y)$ for all $\vert x\vert\lt\vert y\vert\le c\sqrt{n}$, giving $\tilde C_n\ge c\sqrt{n}$.

As I mentioned in my comment - what you are suggesting implies that the probability of being at (3,4) is the same as being at (5,0) for all large n. That seems unlikely, and would guess that $C_n=5$ for n large.

The answer to your modified question is yes, $\tilde C_n$ tends to infinity as n goes to infinity. Also, the rate is at least $\sqrt{\frac45n\log n}$. That is, $$ \liminf_{n\to\infty}\frac{\tilde C_n}{\sqrt{\frac45n\log n}}\ge1. $$ This can be proven by evaluating $p_n(x)$ to leading order in 1/n. (assuming I haven't made any errors below).

The idea is to note that you are repeatedly applying a linear operator, $$ p_{n+1}=Lp_n,\ Lp(x) \equiv (p(x)+p(x-e_1)+p(x+e_1)+p(x-e_2)+p(x+e_2))/5 $$ where $e_1=(1,0)$, $e_2=(0,1)$. In finite dimensional spaces, you would solve this by decomposing $p_0$ into a sum of eigenvectors and for large n, the dominant term of $L^np_0$ will be that corresponding to the largest eigenvalue. In this case, we can diagonalize by applying a Fourier transform $$ p_0(x)=1_{\lbrace x=0\rbrace}=\int_{-[\frac12,\frac12]^2}e^{2\pi ix\cdot u}\\,du. $$ Noting that $e^{2\pi ix\cdot u}$ (as a function of x) is an eigenvector of L, $$ Le^{2\pi ix\cdot u}=\left(\frac15+\frac25\cos(2\pi u_1)+\frac25\cos(2\pi u_2)\right)e^{2\pi ix\cdot u} $$ gives the following for $p_n$, $$ p_n(x)=L^np_0(x)=\int_{[-\frac12,\frac12]^2}\left(\frac15+\frac25\cos(2\pi u_1)+\frac25\cos(2\pi u_2)\right)^ne^{2\pi ix\cdot u}\\,du. $$ The term inside the parentheses is less than 1 in absolute value everywhere away from the origin, so looks like a Dirac delta when raised to a high power n.

The integral can be computed to leading order. After n steps, the standard deviation of the particle's distance from the origin grows of the order of $\sqrt{n}$, so most of its distribution is spread over an area of size the order of n. This means we expect $p_n(x)$ to go to zero at rate 1/n. Let's discard all terms which vanish faster than this.

The range of integration can be replaced by $[-\epsilon,\epsilon]$ for any $0<\epsilon<1/2$, creating an error only of order $e^{-cn}$ (some positive c). On such a range, $$ \frac15+\frac25\cos(2\pi u_1)+\frac25\cos(2\pi u_2)=\exp\left(-\frac{4\pi^2}{5}u^2+O(u^4)\right) $$ Substituting into the integral and changing variables, $$ \begin{align} p_n(x)&=\int_{[-\epsilon,\epsilon]^2}\exp\left(-\frac{4n\pi^2}{5}\vert u\vert^2+O(n\vert u\vert^4)+2\pi ix\cdot u\right)\\,du+O(e^{-cn})\\\\ &=\frac{5}{8n\pi^2}\int_{[-\epsilon'\sqrt{n},\epsilon'\sqrt{n}]^2}\exp\left(-\frac12\vert v\vert^2+O(\vert v\vert^4/n)+i\sqrt{\frac{5}{2n}}x\cdot v\right)\\,dv+O(e^{-cn}) \end{align} $$ Where $\epsilon'=\epsilon\sqrt{8\pi^2/5}$. Rearranging the v^4/n term in the exponential, it is possible to show that it only contributes of the order of 1/n to the integral. Integrals like $\int_{\epsilon\sqrt{n}}^\infty\exp(-x^2/2)\\,dx$ vanish at rate $\exp(-\epsilon^2n/2)/n$, giving $$ \begin{align} p_n(x)&=\frac{5}{8n\pi^2}\int_{\mathbb{R}^2}\exp\left(-\frac12\vert v\vert^2+i\sqrt{\frac{5}{2n}}x\cdot v\right)\\,dv+O(1/n^2)\\\\ &=\frac{5}{4n\pi}\exp\left(-\frac{5\vert x\vert^2}{4n}\right)+O(1/n^2) \end{align} $$ So, if $\vert x\vert\lt\vert y\vert$ then $p_n(x)\gt p_n(y)$ for large n. You can check from this expression that $\vert x\vert\lt c\sqrt{\frac45n\log n}$ (any $c\lt1$) is enough for the first term to dominate in the inequality.

As I mentioned in my comment - what you are suggesting implies that the probability of being at (3,4) is the same as being at (5,0) for all large n. That seems unlikely, and would guess that $C_n=5$ for n large.

The answer to your modified question is yes, $\tilde C_n$ tends to infinity as n goes to infinity. Edit: I'm not so sure now, there was a mistake in my original calculation. What you can say is that there is a positive constant $a$ such that the result holds as long as you restrict to $\vert x\vert^2\lt\vert y\vert^2-a$. Whether or not we can take $a=1$ (which would imply the full result, as $\vert x\vert^2-\vert y\vert^2\in\mathbb{Z}$) requires a more careful consideration of the second order terms in my calculation. In any case, the rate will be of the order of $\sqrt{n}$. That is, $\tilde C_n\ge b\sqrt{n}$ for some constant $b$.

This can be proven by evaluating $p_n(x)$ to leading order in 1/n. (assuming I haven't made any errors below).

The idea is to note that you are repeatedly applying a linear operator, $$ p_{n+1}=Lp_n,\ Lp(x) \equiv (p(x)+p(x-e_1)+p(x+e_1)+p(x-e_2)+p(x+e_2))/5 $$ where $e_1=(1,0)$, $e_2=(0,1)$. In finite dimensional spaces, you would solve this by decomposing $p_0$ into a sum of eigenvectors and for large n, the dominant term of $L^np_0$ will be that corresponding to the largest eigenvalue. In this case, we can diagonalize by applying a Fourier transform $$ p_0(x)=1_{\lbrace x=0\rbrace}=\int_{-[\frac12,\frac12]^2}e^{2\pi ix\cdot u}\\,du. $$ Noting that $e^{2\pi ix\cdot u}$ (as a function of x) is an eigenvector of L, $$ Le^{2\pi ix\cdot u}=\left(\frac15+\frac25\cos(2\pi u_1)+\frac25\cos(2\pi u_2)\right)e^{2\pi ix\cdot u} $$ gives the following for $p_n$, $$ p_n(x)=L^np_0(x)=\int_{[-\frac12,\frac12]^2}\left(\frac15+\frac25\cos(2\pi u_1)+\frac25\cos(2\pi u_2)\right)^ne^{2\pi ix\cdot u}\\,du. $$ The term inside the parentheses is less than 1 in absolute value everywhere away from the origin, so looks like a Dirac delta when raised to a high power n.

The integral can be computed to leading order. After n steps, the standard deviation of the particle's distance from the origin grows of the order of $\sqrt{n}$, so most of its distribution is spread over an area of size the order of n. This means we expect $p_n(x)$ to go to zero at rate 1/n. Let's discard all terms which vanish faster than this.

The range of integration can be replaced by $[-\epsilon,\epsilon]$ for any $0<\epsilon<1/2$, creating an error only of order $e^{-cn}$ (some positive c). On such a range, $$ \frac15+\frac25\cos(2\pi u_1)+\frac25\cos(2\pi u_2)=\exp\left(-\frac{4\pi^2}{5}u^2+O(u^4)\right) $$ Substituting into the integral and changing variables, $$ \begin{align} p_n(x)&=\int_{[-\epsilon,\epsilon]^2}\exp\left(-\frac{4n\pi^2}{5}\vert u\vert^2+O(n\vert u\vert^4)+2\pi ix\cdot u\right)\\,du+O(e^{-cn})\\\\ &=\frac{5}{8n\pi^2}\int_{[-\epsilon'\sqrt{n},\epsilon'\sqrt{n}]^2}\exp\left(-\frac12\vert v\vert^2+O(\vert v\vert^4/n)+i\sqrt{\frac{5}{2n}}x\cdot v\right)\\,dv+O(e^{-cn}) \end{align} $$ Where $\epsilon'=\epsilon\sqrt{8\pi^2/5}$. Rearranging the v^4/n term in the exponential, it is possible to show that it only contributes of the order of 1/n to the integral. Integrals like $\int_{\epsilon\sqrt{n}}^\infty\exp(-x^2/2)\\,dx$ vanish at rate $\exp(-\epsilon^2n/2)/n$, giving $$ \begin{align} p_n(x)&=\frac{5}{8n\pi^2}\int_{\mathbb{R}^2}\exp\left(-\frac12\vert v\vert^2+i\sqrt{\frac{5}{2n}}x\cdot v\right)\\,dv+O(1/n^2)\\\\ &=\frac{5}{4n\pi}\exp\left(-\frac{5\vert x\vert^2}{4n}\right)+O(1/n^2) \end{align} $$ So, if $\vert x\vert\lt\vert y\vert$ then $p_n(x)\gt p_n(y)$ for large n. You can check that if $\vert x\vert\lt\vert y\vert-a$ for some constant a, $$ \begin{align} p_n(x)-p_n(y)&=\frac{5}{4n\pi}e^{-\frac{5\vert x\vert^2}{4n}}(1-e^{\frac{5(\vert x\vert^2-\vert y\vert^2)}{4n}})+O(1/n^2)\\\\ &\ge\frac{5}{4n\pi}e^{-\frac{5\vert x\vert^2}{4n}}(1-e^{-\frac{5a}{4n}})+O(1/n^2)\\\\ &\ge\frac{25}{4n^2\pi}\left(ae^{-\frac{5\vert x\vert^2}{4n}}+O(1)\right) \end{align} $$ To be certain that this is positive, we would need to evaluate the constant term in the final expression and compare to $a$.