Random Walk on Pentagonal Tiling I’ve recently been looking at closed walks on tilings of the plane in which the “player” can move from one tile to any of its edge-adjacent neighbors. In particular, I’m trying to find asymptotic formulae for various tilings giving the number of n-step closed paths starting and ending at an arbitrary tile. For an infinite square checkerboard we have
$$\sim \frac{4^n}{n}\cdot\frac{2}{\pi}$$
and for an infinite hexagonal honeycomb we have
$$\sim \frac{6^n}{n}\cdot \frac{\sqrt{3}}{2\pi}$$
and for an infinite tiling of equilateral triangles:
$$\sim \frac{3^n}{n}\cdot\frac{3}{\pi}$$
where the first and third formulae are applicable only for even $n$, since there are no closed paths of odd length on the square and triangular tesselations. All of these can be derived using the method employed in this answer.
QUESTION: Can anybody derive analogous formulae for less “nice” tilings, like the cairo pentagonal tiling, the floret pentagonal tiling, the tetrakis square tiling, the snub square tiling, or any others? Clearly they should be in the form of a constant times $c^n/n$, where $c$ is a constant depending on the tiling that is more easily calculated ($c=5$ for cairo and floret, $c=3$ for tetrakis square, and $c=\frac{1+\sqrt{33}}{2}$ for the snub square).
WHY I CAN’T SOLVE IT: The method I used to derive these 3 formulae relied heavily on the fact that all tiles in the tiling had the same orientation, and that the same “moves” could be made no matter what tile is occupied (for the square tiling, the player can always move N,E,S,W, and for the hexagonal tiling, the player can always move in the same 6 directions). This isn’t actually true for the triangular tiling, but luckily the set of moves that the player can make alternates from one move to the next, so I was able to use the same method with a slight modification. This is because there are two orientations of equilateral triangles in this tiling, and all triangles edge-adjacent to a triangle of one orientation are of the opposite orientation.
This problem has been slowly killing me for a long time, partly because it seems like it should be easy but isn’t (for me, at least). Help?

NOTE: Since posting the question, I have come up with the formula for another tiling that uses only regular hexagons and equilateral triangles. The formula applies to walks starting in a hexagon and is valid only for even $n$:
$$\sim \frac{(3\sqrt{2})^n}{n}\cdot \frac{9}{4\pi}$$
 A: Another standard approach is to use a multivariate local central limit theorem, which can be found e.g. in Petrov's book "sums of independent random variables" and the paper "Local Limit Theorems for Lattice Random Variables" by A. B. Mukhin https://doi.org/10.1137/1136086
I especially like "An elementary proof of the local central limit theorem B Davis, D McDonald - Journal of Theoretical Probability, 1995 " (see http://www.stat.purdue.edu/research/technical_reports/pdfs/1993/tr93-41.pdf 
(These references are from a recent discussion at Fair partitioning of a set - Weighted sums of Bernoullis).
If all the tiles have the same orientation, the local CLT can be applied directly.
Otherwise, by examining only the times when a translate of the original tile is visited, a random walk is obtained so the local CLT applies- the remaining issue is to compute the covariance matrix. 
In the tiling http://www.math.cmu.edu/~bkell/21110-2010s/3-6-3-6.png  one just has to go in skips of two. In the general case, there will be k tile types, and the random walk defines a Markov chain on these tile types. The covariance matrix can then be obtained by averaging on the starting tile (according to the stationary distribution) the covariance matrix of the step from that tile.  
See also Keilson, Julian, and D. M. G. Wishart. "A central limit theorem for processes defined on a finite Markov chain." In Mathematical Proceedings of the Cambridge Philosophical Society, vol. 60, no. 3, pp. 547-567. Cambridge University Press, 1964.
A: Let me outline a possible approach, which is very standard but I don't know a reference for this particular problem. Start with any doubly periodic embedded planar graph $\Gamma$, that is, a graph with $\mathbb{Z}^2$ acting on it by shifts. Let $\Gamma_1:=\Gamma/\mathbb{Z}^2$, a fundamental domain wrapped on a torus. Let us also take a large $N$ and let $\Gamma_N=\Gamma/(N\cdot \mathbb{Z}^2)$, an $N\times N$ torus. We will count walks on $\Gamma_N$ and then translate the results to walks on $\Gamma$. 
If $A_N$ denotes the adjacency matrix of $\Gamma_N$, then the number of walks of length $n$ from $z$ to $w$ is $(A^n_N)_{zw}$. Since $A_N$ is symmetric, it can be diagonalized; so, let $\lambda_1,\dots,\lambda_M$ be the eigenvalues and  $f_1,\dots,f_M:V(\Gamma_N)\to \mathbb{C}$ be the corresponding eigenfunctions, normalized so that $\sum_z |f_i(z)|^2=1$, where $M=|\Gamma_1|N^2$. Then, we have
$$
(A^n_N)_{zw}=\sum_i \lambda_i^nf_i(z)f_i(w).
$$
Note that $A_N$ commutes with the action of $\mathbb{Z}^2$, so it must have the eigenfunctions that are also eigenfunctions of the shifts by $1$, both in horizontal and vertical direction. The general form of such an eigenfunction is 
$$
f(z)=e^{2\pi i p x /N }e^{2\pi i q y/N}g(\hat{z}),
$$
where $p,q\in\{0,\dots,N-1\}$ and $\hat{z}$ is a projection of $z$ onto $\Gamma_1$, and $z=x+iy$. Looking for the eigenfunctions $f$ of $A_N$ of this form boils down to the computation of the eigenvalues of $A_1^{\varphi\theta}$, where $\varphi,\theta \in [0,2\pi]$ are parameters and $A^{\varphi\theta}_1$ is the weighted adjacency matrix of $\Gamma_1$ defined by
$$
(A^{\varphi\theta}_1)_{\hat{z}\hat{w}}= e^{i\varphi(\Re w-\Re z)}e^{i\theta(\Im w-\Im z)}(A_1)_{\hat{z}\hat{w}},
$$
and the representatives $z$ and $w$ are chosen in such a way that they are adjacent whenever $\hat{z}$ and $\hat{w}$ are adjacent. Note that alternatively, you could move all the dependence of $\varphi$ and $\theta$ into the edges that intersect the boundary of the fundamental domain, and then the resulting matrix algebraically depends on $e^{i\varphi}$ and $e^{i\theta}$.
Anyway, assuming that $\lambda^{\varphi\theta}_1,\dots,\lambda^{\varphi\theta}_{|\Gamma_1|}$ are the eigenvalues of $A^{\varphi\theta}_1$, and $g^{\varphi\theta}_i$ the corresponding (normalized) eigenfunctions; these are algebraic functions of $e^{i\varphi}$ and $e^{i\theta}$; concretely, lambdas are the roots of the polynomial $P(\lambda,{e^{i\varphi},e^{i\theta}}):=\det(\lambda I-A^{\varphi\theta}_1
)$ that is also a Laurent polynomial in the last two arguments. Then, we can approximate the sum for $A^n_N$ by an integral, to get
$$
\lim_{N\to\infty} (A^n_N)_{zz} = \frac{1}{(2\pi)^2}\int_{\varphi=0}^{2\pi}\int_{\theta=0}^{2\pi}\sum_{i=1}^{|\Gamma_1|}(\lambda^{\varphi\theta}_i)^n |g^{\varphi\theta}_i(z)|^2d\varphi d\theta.
$$
The asymptotics of this expression as $n\to\infty$ is given by the Laplace method. So, our next task is to find $\varphi$ and $\theta$ (and $i$) that maximize $|\lambda^{\varphi\theta}_i|$. I believe that it should be possible to prove that the unique maximum is always given by $(\varphi,\theta)=(0,0)$ (basically, because $(A^{\varphi\theta}_1)^n$ counts n-step walks on $\Gamma_1$ with complex weights of absolute value $1$, and that count always gives a smaller number than the unweighted count); there might be several maximising $i$'s in the case of periodicity phenomena. Also, from the explicit form of the polynomial $P$, it is possible to calculate the second derivative of $\lambda_{\text{max}}^{\varphi\theta}$. And that is all you need to run the Laplace method. Note also that if you are fine with averaging the quantity of interest over the fundamental domain, then you don't even need to calculate the eigenfunction $g^{00}_{\text{max}}$. 
This also explains nicely several general features of the resulting formulae. In the generic case (when the second derivative of $\lambda^{\theta\varphi}_{\text{max}}$ is non-degenerate), the Laplace method outputs $\sim C\cdot c^n/n$. The number $c$ is always algebraic, and the number $C$ is always algebraic divided by $\pi$ (because the second derivative is algebraic, $1/\pi^2$ is in front of the integral, and $\pi$ comes from the Laplace method)
