A geometric law of large numbers Question: Choose independent uniformly distributed random variables $X$ and $Y$ on a closed Riemannian manifold $(M,g).$ The geodesic midpoint of $X$ and $Y$ is well-defined a.e., and is distributed according to some probability measure $\Pi_1$ on $M.$ We can then repeat the process and choose two independent random variables distributed according to $\Pi_1,$ and the geodesic midpoint of these will again be defined a.e. and will be distributed according to a new probability measure $\Pi_2$ on $M.$
If we iterate this process, does it converge and can we characterize the limit?
Hypothesized solution for convex surfaces: For convex hypersurfaces $M$ in $\mathbb R^d$ I hypothesize the sequence of measures converges to a volume form associated with a round spherical metric on $M.$
Further detail for motivation
Law of large numbers in Euclidean space
Let $X$ and $Y$ be independent identically distributed random variables on the real line with probability distribution $\Pi_1.$ The midpoint $(X+Y)/2$ of the line joining $X$ to $Y$ is a new random variable with distribution $\Pi_2$ (which can be computed from $\Pi_1$ by convolution and scaling). Now pick two independent random variables distributed according to $\Pi_2$ and iterate the procedure.
The law of large numbers says precisely that the sequence of distributions $\Pi_1,\Pi_2,\dots$ so defined will converge to a point mass at the mean of $\Pi_1.$ Indeed, $\Pi_N$ is just the distribution of $(X_1+\cdots + X_{2^N})/2^N,$ where $X_1,\dots, X_{2^N}$ are IID with distribution $\Pi_1.$
We may even characterize the rate of convergence through the central limit theorem, which can be proved using Fourier analysis for instance.
Set-up on a general manifold
The same arguments are fine in Euclidean spaces of all dimensions. But this procedure may in fact be implemented in any Riemannian manifold $(M,g),$ using the metric to define the midpoint between two points. Let us keep it simple and assume that $M$ is a closed manifold, and let $\Pi$ be the uniform probability measure on $M$ (the measure obtained by rescaling the volume element coming from $g$ to have unit mass). If $X$ and $Y$ are independent uniformly distributed random variables on $M,$ then we may define $Z$ to be the midpoint of the shortest geodesic segment between $X$ and $Y.$ The point $Z$ is well-defined a.e. (likely some application of Sard's theorem to the exponential function will prove this), and determines a new probability distribution $\Pi_1$ on $M$ that is absolutely continuous with respect to the uniform measure $\Pi.$
We may now iterate this process, choosing two independent random variables distributed according to $\Pi_1$ and letting $\Pi_2$ be the distribution of their geodesic midpoint (which is defined a.e. by absolute continuity of $\Pi_1$ with respect to $\Pi$), and so on. The result is a sequence $\Pi_1,\Pi_2,\Pi_3,\dots$ of probability measures on $M,$ each absolutely continuous with respect to volume. The question is whether this sequence converges, to what, and if possible how fast.
Easy examples
Any curve On a closed $1$-dimensional manifold, there is no intrinsic geometry and the evolution plainly results in the uniform measure for each iteration. In other words, $\Pi = \Pi_1 = \Pi_2 = \cdots,$ where $\Pi$ is the uniform measure on $M.$ The reason is that the operations preserve the translational symmetry under which all closed curves are invariant.
Constant positive curvature: spheres of any dimension In case of a sphere $S^d,$ the procedure also plainly results in the uniform measure for each iteration. This must be true because the operation of picking the midpoint is invariant under the sphere's rotational symmetries and so the resulting measure must be as well.
Constant zero curvature: flat torus of any dimension On the flat torus, I expect the story to be the same based on symmetry but have not reasoned through it carefully.
Non-trivial example 1: sphere with a spike Let $M$ be a $2$-dimensional sphere with a large thin spike protruding from it. Assume for simplicity that $M$ has unit volume, and let $\epsilon$ be the total area of the spike. Using the notation from the problem statement, we can see that the first measure $\Pi_1$ assigns smaller area to the spike, of the order $\epsilon^2$: speaking approximately, the midpoint of $X$ and $Y$ will live on the spike only when $X$ and $Y$ both live on the spike. This happens with probability $\epsilon^2$ given $X$ and $Y$ are independent.
Non-trivial example 2: dumbbell with a long neck Let $M^2$ be the combination of two $2$-dimensionsal spheres that are very far apart connected by a very thin tubular neck. Assume $M$ has unit volume and the neck has total area $\epsilon.$ We will see that the area of the neck grows under our procedure to have area approximately $1/2+\epsilon - \epsilon^2.$
Let $X,Y,$ and $\Pi_1$ be as in the problem statement. Then, roughly speaking, the geodesic midpoint of $X$ and $Y$ is contained on the neck whenever $X$ and $Y$ are not contained on the same sphere. Each sphere has volume roughly $(1-\epsilon)/2$ and the chance that $X$ and $Y$ are contained together in one of them is roughly $(1-\epsilon)^2/2.$ Therefore the complementary probability is $1-(1-\epsilon)^2/2 = 1/2 +\epsilon - \epsilon^2,$ which is roughly the chance that the midpoint of $X$ and $Y$ is found on the neck.
Some possible extensions
Non-compact manifolds The same procedure makes sense on any Riemannian manifold equipped with a probability measure that is absolutely continuous with respect to the volume induced by the metric. (Perhaps all that is needed is a metric measure space.) Therefore we may ask the same questions on a non-compact manifold provided we first choose an absolutely continuous probability measure, as we do in Euclidean space
A geometric evolution based on Minkowski's theorem for convex bodies The Minkowski problem in geometry assures us that a convex body in $\mathbb R^d$ is uniquely determined by its Gaussian curvature (or by the push-forward of volume to the sphere under the Gauss map). Now begin our iteration with a convex manifold $M$ in $\mathbb R^d.$ Choose $X$ and $Y$ on the sphere $S^d$ to be independent and identically distributed according to the pushforward of $M$'s volume measure under the Gauss map. We may now carry out our midpoint procedure on the sphere resulting in a new measure $\Pi_1$ on the sphere that is absolutely continuous with respect to volume. The first question is whether this new measure meets the conditions of Minkowski's theorem and determines a convex body (its centroid must be at the origin and it must not be concentrated on a great sub-sphere). Assuming it does, we obtain a new convex body with each iteration of the procedure, and we may ask two questions

*

*Does a convex body converge to the sphere under this procedure?

*Is there a 'natural' continuous geometric evolution equation which 'extends' this discrete evolution?

 A: Updated more quantitatively:
This repeated averaging can also lead to a distribution concentrated at a single point, even if $M$ is topologically a 2-sphere.
Consider the unit sphere (or the earth); let $N$ be the northern hemisphere, and let $C$ be the convex hull of the equator. So $N$ is a positively curved surface and $C$ is a flat surface going through the center of the sphere. Let $M=N\cup C$, or some smoothed version.
Since $\text{area}(N)=2 \pi$ and $\text{area}(C)=\pi$, the uniform distribution on $M$ assigns a weight of $\frac23$ to $N$ and $\frac13$ to $C$. But repeated averaging would lead to almost all the weight going to the center of $C$.
The key insight is that in this manifold, $N$ is not convex. For some pairs of points just north of the equator, the shortest path connecting them in $M$ will go through $C$. (In the earthly $M$, the shortest path from Aceh to Bogotá would not be the great-circle path through Europe, but the path going down to the equator and then close to the center of the earth.)

Similarly, take a blue strip and a green strip, each of area $\epsilon$ with points roughly $\pi/6$ away from the equator in $N$ and $C$ respectively. What is the chance that a random blue-green pair of points will have a midpoint in $C$? Points with similar longitudes will have midpoints roughly evenly split between $N$ and $C$. But for points with different longitudes, the angular distance will be traversed more quickly in $C$ than in $N$ (since the path will stay closer to the axis of the earth), and the midpoints will be in $C$ more often than in $N$.
We can analyze this more quantitatively using the area-preserving map
$f:[-\frac12,1] \times [0, 2\pi] \to M$ given by
$$f(u,v)=\left(\sqrt{1+\min(2u,-u^2)},\ v,\ \max(u,0)\right)$$
We can identify several regions where the midpoint of $f(u,v)$ and $f(u',v')$ is in $C$:
In $11$% (or $1/9$) of the cases, both points are in $C$, so the midpoint is obviously in $C$.
In ~$34$% of the cases, one point is in $N$, one point is in $C$, and the midpoint is in $C$. (This came from a numerical test: given $n\in N$ and $c\in C$, I checked many points $e$ on the equator, to see whether $d_N(n,e)<d_C(e,c)$ for the $e$ that minimizes $d_N(n,e)+d_C(e,c)$.)
In ~$2$% of the cases, both points are in $N$ and (following the earthly argument)
$$d_N(f(u,v),f(u',v')) > \arcsin(u)+\arcsin(u')+2\sin(\frac{v-v'}{2})$$
so the midpoint is in $C$.
In total at least $47$% of the midpoints are in $C$ after one averaging, even though $C$ started with only $33$% of the weight.
Repeating this will lead to a distribution with almost all the weight in $C$, and eventually with almost all the weight near the center of $C$, the center of the sphere.
For another example, consider a cone with its base. This has a single edge, and is flat elsewhere, so the calculations are (difficult with the cut locus but) as simple as possible. Repeated averaging on this surface will also concentrate a uniform distribution on the base and then in the center of the base.
