Does a random walk on a surface visit uniformly? Let $S$ be a smooth compact closed surface embedded in $\mathbb{R}^3$ of genus $g$.
Starting from a point $p$, define a random walk as taking discrete steps
in a uniformly random direction,
each step a geodesic segment of the same length $\delta$.
Assume $\delta$ is less than the injectivity radius and small with respect to the intrinsic diameter of $S$.

Q. Is the set of footprints of the random walk evenly distributed on $S$,
in the limit? By evenly distributed I mean the density of points per unit area of surface
is the same everywhere on $S$.

This is likely known, but I'm not finding it in the literature on random walks
on manifolds. I'm especially interested in genus $0$.
Thanks!

Update (6JUn2021). The answer to Q is Yes, going back
38 years to Toshikazu Sunada, as recounted in @RW's
answer.
 A: This random walk is known in the literature as the "geodesic random walk". For a manifold with positive curvature, theorems 1 and 4 of arXiv:1609.02901 prove that the uniform measure on the manifold is the unique stationary distribution of the geodesic walk.
A: Edit: As far as I understand it, this is the approach of T. Sunada, as described in the paper linked in R W's answer.
Let us fix $\delta>0$ such that any pair of points can be joined by a finite sequence of steps of size $\delta$. For instance, we may choose any $\delta$ less than the injectivity radius, as I discuss at the end of this answer.
It is a purely measure-theoretic result that the uniform distribution of the walk holds for almost all initial points. As I describe below, this implies, provided $\delta$ is at most the injectivity radius, that it actually holds for all initial points in the following sense. Call $\mathbb P_{x_0}$ the distribution of the random walk started at $x_0$, and $\mu$ the uniform measure on $S$.

For all $x_0\in S$, for all $f_0:S\to\mathbb R$ continuous, we have
$$ \lim_{n\to\infty}\frac1n\sum_{i<n}f_0(x_i) = \int f\mathrm d\mu $$
$\mathbb P_{x_0}$-almost surely.

I don't think the proof uses anything more than $S$ being a closed connected Riemannian manifold.
Ergodic argument
Let $\Omega\subset S^{\mathbb N}$ be the set of sequences such that two consecutive terms are at distance precisely $\delta$. This is a compact space (closed subset of a compact space, by Tychonoff or simply using a convenient metric), and it carries a natural probability $\mathbb P_\mu$, corresponding to $x_0$ being distributed according to $\mu$ and the following steps according to the random walk rules. The shift operator
$$ T:(x_0,x_1,\ldots)\mapsto(x_1,x_2,\ldots) $$
is such that $T_*\mu=\mu$ (let us just accept this until the end of the proof sketch), so $(\Omega,\mathbb P_\mu,T)$ is a dynamical system in the measure-theoretic sense. According to Birkhoff's ergodic theorem, for $\mathbb P_\mu$-almost every $x$, we have
$$ \lim_{n\to\infty}\frac1n\sum_{i<n}f(T^ix) = \mathbb E[f|\mathcal F_T](x) $$
for all $f:\Omega\to\mathbb R$ continuous, where $\mathcal F_T$ is the algebra of $T$-invariant sets. Now if $(\Omega,\mathbb P_\mu,T)$ is ergodic, then
$$\mathbb E[f|\mathcal F_T](x) = \int f\mathrm d\mathbb P_\mu = \int f_0\mathrm d\mu $$
for all $f$ depending only on $x_0$, i.e. $f:x\mapsto f_0(x)$. Thus my claim will follow by Fubini.
It remains to prove that $\mu$ is $T$ invariant, and that the resulting system is ergodic. According to the probabilistic interpretation in terms of Markov chains, it suffices to show that $x_1$ has the same distribution as $x_0$ (the uniform one) under $\mathbb P_\mu$. This is a consequence of the fact that the geodesic flow leaves the measure induced on $TM$ invariant, which itself is a consequence of Liouville's theorem because the geodesic flow is Hamiltonian.
Now let us show the system is ergodic. Let $f:\Omega\to\mathbb R$ be $T$-invariant, i.e. $f\circ T = f$, and suppose also that $f$ depends only on the first $k$ terms of the sequence. Let us show that it is constant by choosing $x$ and $y$ arbitrary and showing $f(x)=f(y)$. Because of the hypothesis that two points can be linked using steps of size $\delta$, we know that we can go from $x_k$ to $y_{k-1}$ in $N$ steps for some finite $N$. Killing the first $k$ terms of $x$ using $T$, then adding $N$ terms linking $x_k$ to $y_{k-1}$, and adding the first $k$ terms of $y$, we see that $f(x) = f(y_k)$, where $y_k$ is equal to $y$ up to the $k$th term. Because $f$ only depends on the first $k$ terms, we have in fact $f(x)=f(y)$. Now an approximation argument shows that the system is ergodic (see for instance a previous answer of mine here).
From almost all to all
Now this extends I think to all points in the manifold, provided $\delta$ is less than the injectivity radius. Let $A$ be the set of points $x_0$ such that the random walk started at $x_0$ is evenly distributed almost surely, and $\mathbf 1_A$ its indicator function. Then
$$ \mathbb P_{x_0}(\text{even distribution}) = \int \mathbb P_{x_2}(\text{even distribution})\mathrm d\mathbb P_{x_0} \geq \int\mathbf 1_A\mathrm d (T^2_*\delta_{x_0}) = 1, $$
because the distribution of $x_2$ under $\mathbb P_{x_0}$ is continuous with respect to the Lebesgue measure. Indeed, it is the image of the Lebesgue measure on the product of two sphere under a map that is a submersion almost everywhere. I have convinced myself of this last fact, but I can give more detail if people are interested.
All points can be reached
Define an equivalence relation so that $x\sim y$ if there is a chain of steps of length $\delta$ from $x$ to $y$. It is clearly an equivalence relation provided we allow for the trivial chain $(x)$ to be a witness for $x\sim x$. To prove that we can reach any given point from any other, it will suffice to show that the equivalence classes are open. We will show that if $\delta$ is less than the injectivity radius $r_\text{inj}$, then all points $y$ with $d(x,y)<\min(\delta,r_\text{inj}-\delta)$ can be reached from $x$ in two steps. Since $S$ is compact, this will actually show that there is an upper bound on the minimal number of steps needed to link two points of the manifold.
It will be enough to show that for any such pair $(x,y)$, the spheres $\partial B_x(\delta)$ and $\partial B_y(\delta)$ intersect, which in turn would follow from the existence of points $z_\pm$ on $\partial B_y(\delta)$ such that $\pm d(x,z_\pm)>\pm\delta$ (this sphere is topologically a sphere of dimension at least two, so it is path connected). But such points are easy to find. Follow the (unique) minimising geodesic from $y$ to $x$, continue until the geodesic has length $\delta$, and call this endpoint $z_-$. Since $x$ is on the minimising geodesic segment from $y$ to $z_-$, the distance between $z_-$ and $x$ is at most $\delta$. In the other direction, follow the geodesic from $x$ to $y$, continue until the geodesic has length $d(x,y)+\delta$, and call this endpoint $z_+$. Because the geodesic is minimising, $\delta=d(y,z_+)<d(x,z_+)$. This concludes the argument.
A: This problem was first considered and solved by Sunada, see his 1983 paper Mean-value theorems and ergodicity of certain geodesic random walks. Alas, the authors of the quoted arxiv paper were not aware of this. Any assumptions on curvature and dimension are not necessary - it is just enough to assume that the manifold is compact. As it has been pointed out by Pierre PC, the fact that the Riemannian volume is a stationary measure is an immediate consequence of the Liouville theorem. Its ergodicity with respect to the geodesic random walk is equivalent to the absence of invariant subsets of the manifold, which would follow, for instance, if any two points can be joined by a chain of geodesic segments of length $\delta$. Actually, geodesic random walks are always mixing for sufficiently small $\delta$ - this is a consequence of the fact that the cube of the transition operator has a density which is bounded away from 0 on the diagonal. The latter also implies the uniqueness of the stationary measure.
EDIT I had misattributed the term "geodesic random walk" to Sunada. Actually, it seems to be first introduced in 1975 by Jørgensen The central limit problem for geodesic random walks whose work is quoted by Sunada.
