A particle starts in a cell in an infinite hexagonal grid, and every second, jumps to an adjacent cell uniformly randomly. What is the expected amount of time until the particle is $N$ cell-jumps away from its starting point? With some linear algebra, for example, one finds values of $1$, then $10/3$, then $213/29$, for cases $N=1,2,3$ respectively. Computer simulation shows growth is approximately $4N^2/5$.
I expected to be able to solve this problem with similar methods (using polynomials in barycentric coordinates, constrained by dihedral symmetries) as to my recent Puzzling question, but to no avail so far. Curiously, by a coupling argument, this problem is equivalent to computing the expected value of the variable $\text{min}\{X_1,X_2\}$ where $X_i$ are iid variables representing the escape time of the honeybee from the center of its triangle in the linked problem, but that observation doesn't seem to help much.
Some rambling about my current attempts: In barycentric coordinates $(\alpha, \beta, \gamma)$ whereby we always have $\alpha + \beta + \gamma = 3N$, it would seem reasonable to demand that—in order to find the average escape time at $(\alpha, \beta, \gamma)$ from the $N-1$-hexagon centered at $(N,N,N)$—we find a function $H(\alpha, \beta, \gamma)$ algebraically satisfying the "average-of-6-neighbors-plus-1" property everywhere, which also satisfies $H = 0$ whenever $\alpha = 0, 2N$ or $\beta = 0, 2N$ or $\gamma = 0, 2N$.
After all, this approach is exactly how the triangular escape time problem is solved, just leaving out the $2N$ constraints. In that case, we think of the elementary symmetric polynomials in $\alpha, \beta, \gamma$, and realize $\alpha\beta\gamma$ is a good candidate. It doesn't quite satisfy the averaging-plus-one law—its difference from its nearby-average function is $3N$ and not $1$—so we tweak it to $\frac{3\alpha\beta\gamma}{\alpha+\beta+\gamma}$ to solve the problem.
So that is how I proceeded here, examining the obvious candidate $H=\alpha \beta \gamma (\alpha-2\beta-2\gamma)(\beta-2\alpha-2\gamma)(\gamma-2\alpha-2\beta)$. But its difference from its nearby-average function is gnarly, and not susceptible to obvious tweaks. With some thought, one realizes the field of rational functions invariant up to angular and mirror symmetry are generated by $H$ as well as $e_1 = \alpha+\beta+\gamma$ and $e_2 = \alpha \beta + \alpha \gamma + \beta\gamma$. Especially considering the empirical evidence that our formula will be of degree $2$, one might try candidate tweaks like $\frac{H}{e_1^4}$ or $\frac{H}{e_1^2 e_2}$ or $\frac{H}{e_2^2}$ or $\frac{H^2}{e_1^4 e_2^3}$… but some time spent in Mathematica proved fruitless.
It's become clear to me now that no rational function of the form $\frac{F}{e_1^n e_2^m}$ will satisfy the criteria of the first paragraph, because such a function will still be defined on and inside the full triangular region, thus restricting to a solution of the honeybee escape time problem. By standard Markov chain reasoning, this solution is unique, and obviously not the solution to the problem at hand. So, either an even more complicated denominator is needed (one giving poles outside the hexagon but inside the triangle), or we need to allow possibilities like $H \neq 0$ even if $\alpha = 0$ as long as we're outside the hexagonal boundary, or we need some even more radical change to our techniques.
Best Answer
Let us encode the hexagonal grid using the hexagonal lattice
$$ \mathsf{G} = \{ a + b \omega : a, b \in \mathbb{Z} \}, \qquad \omega = e^{i\pi/3},$$
where each $z \in \mathsf{G}$ represents the center of a hexagonal cell. Then two cells $z_1$ and $z_2$ are adjacent precisely when $\left| z_1 - z_2 \right| = 1$.
We also write $\mathsf{C}_n$ for the set of all cells with are precisely $n$ cells away from the origin.
Now let $(X_n)_{n\geq0}$ denote the simple random walk on $\mathsf{G}$, started at $X_0 = 0$. Denote by $\tau_n$ the hitting time of $\mathsf{C}_n$. Then by the second Wald's identity, the expectation of $\tau_n$ is
$$ \mathbb{E}[\tau_n] = \mathbb{E}\bigl[ \left| X_{\tau_n} \right|^2 \bigr]. $$
Now, if we define the continuous-time process $\tilde{X}^{(n)}_t = \frac{1}{n} X_{\lfloor n^2 t\rfloor}$ by the diffusive scaling of $X$, then by the invariance principle, $\tilde{X}^{(n)}$ converges to the complex Brownian motion $W$ started at $0$. So, if $\ell$ denotes the constant factor appearing in the asymptotic formula for $\mathbb{E}[\tau_n]$, then
$$ \ell = \lim_{n\to\infty} \frac{1}{n^2}\mathbb{E}[\tau_n] = \mathbb{E}\bigl[ \left| W_{\tau} \right|^2 \bigr] = \int_{\mathsf{C}} \left| z \right|^2 \, \mathbb{P}(W_{\tau_{\mathsf{C}}} \in \mathrm{d}z), $$
where $\mathsf{C}$ is the regular hexagon with vertices $e^{ik\pi/3}$ for $k = 0, 1, \dots, 5$, which arises as the "limit" of the rescaled set $n^{-1}\mathsf{C}_n$, and $\tau_{\mathsf{C}}$ is the hitting time of $\mathsf{C}$.
In order to compute the last integral, consider the Schwarz–Christoffel mapping
$$ \phi(z) = K \int_{0}^{z} \frac{1}{(1-\zeta^6)^{1/3}} \, \mathrm{d}\zeta $$
over the unit disc $\mathbb{D}$, and the normalizing factor $K$ is chosen as
$$ K = \left( \int_{0}^{1} \frac{1}{(1-x^6)^{1/3}} \, \mathrm{d}x \right)^{-1} = \frac{6 \cdot 2^{1/3} \pi^{1/2}}{\Gamma(\frac{1}{6})\Gamma(\frac{1}{3})} $$
so that $\phi(1) = 1$ holds. It is well-known that $\phi$ maps $\partial\mathbb{D}$ to $\mathsf{C}$, and $\phi$ is a conformal mapping from $\mathbb{D}$ to the interior of $\mathsf{C}$. So by the conformal invariance of $W$, we obtain
\begin{align*} \ell &= \int_{\partial\mathbb{D}} \left| \phi(w) \right|^2 \, \mathbb{P}(W_{\tau_{\partial\mathbb{D}}} \in \mathrm{d}w) = \frac{1}{2\pi} \int_{0}^{2\pi} \bigl| \phi(e^{i\theta}) \bigr|^2 \, \mathrm{d}\theta \\ &= K^2 \sum_{n=0}^{\infty} \binom{-1/3}{n}^2 \frac{1}{(6n+1)^2} \approx 0.80957626278006891494. \end{align*}