It is a purely measure-theoretic result that this holds for almost all initial points. As I describe below, I think this implies 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 fd\mu $$ $\mathbb P_{x_0}$-almost surely.

The fact that $\delta$ is small enough plays no role in the first part.

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 fd\mathbb P_\mu = \int f_0d\mu $$ for all $f$ depending only on $x_0$, i.e. $f:x\mapsto f_0(x)$. Thus my claim will follow by Fubini.

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)$. Let $N$ be a number large enough that $N\delta$ is larger than the diameter of the manifold. Killing the first $k+N$ 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).

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})d\mathbb P_{x_0} = \int 1d 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 (image of the Lebesgue measure on the product of two sphere under a submersion almost everywhere). In fact the result is true for any closed Riemannian manifold, as maybe we could have expected, although we have to look for $x_d$, where $d$ is the dimension of the manifold. The “I think” I use comes from the fact that although $x_d$ is the image of a map from a product of $d$ spheres to $S$, I cannot give a simple argument of the fact (natural to me) that this is a submersion at all points that form a basis of $\mathbb R^d$.

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.

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.

It is a purely measure-theoretic result that this holds for almost all initial points. As I describe below, I think this implies 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 fd\mu $$ $\mathbb P_{x_0}$-almost surely.

The fact that $\delta$ is small enough plays no role in the first part.

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 fd\mathbb P_\mu = \int f_0d\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)$. Let $N$ be a number large enough that $N\delta$ is larger than the diameter of the manifold. Killing the first $k+N$ 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).

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})d\mathbb P_{x_0} = \int 1d 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 (image of the Lebesgue measure on the product of two sphere under a submersion almost everywhere). In fact the result is true for any Riemannian manifold, as maybe we could have expected, although we have to look for $x_d$, where $d$ is the dimension of the manifold. The “I think” I use comes from the fact that although $x_d$ is the image of a map from a product of $d$ spheres to $S$, I cannot give a simple argument of the fact (natural to me) that this is a submersion at all points that form a basis of $\mathbb R^d$.

It is a purely measure-theoretic result that this holds for almost all initial points. As I describe below, I think this implies 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 fd\mu $$ $\mathbb P_{x_0}$-almost surely.

The fact that $\delta$ is small enough plays no role in the first part.

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 fd\mathbb P_\mu = \int f_0d\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)$. Let $N$ be a number large enough that $N\delta$ is larger than the diameter of the manifold. Killing the first $k+N$ 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).

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})d\mathbb P_{x_0} = \int 1d 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 (image of the Lebesgue measure on the product of two sphere under a submersion almost everywhere). In fact the result is true for any closed Riemannian manifold, as maybe we could have expected, although we have to look for $x_d$, where $d$ is the dimension of the manifold. The “I think” I use comes from the fact that although $x_d$ is the image of a map from a product of $d$ spheres to $S$, I cannot give a simple argument of the fact (natural to me) that this is a submersion at all points that form a basis of $\mathbb R^d$.