Suppose $\mathcal M\subseteq\mathbb R^n$ is a submanifold of $\mathbb R^n$. Take $d(\cdot,\cdot):\mathcal M\times\mathcal M\to\mathbb R_+$ to be the pairwise geodesic distance function, so $d(x,y)$ is the length of the shortest path from $x\in\mathcal M$ to $y\in\mathcal M$.
Holding one coordinate fixed and varying the other, wherever $d(\cdot,\cdot)$ is differentiable it satisfies the eikonal equation $\|\nabla f(x)\|_2=1$. In particular, we can write two conditions satisfied by $d(\cdot,\cdot)$ almost everywhere: $$\|\nabla_x d(x,y)\|_2=1\qquad\textrm{and}\qquad\|\nabla_y d(x,y)\|_2=1.$$
In a sense, these two conditions are redundant. If $d(\cdot,\cdot)$ satisfies the first condition, it "looks" like a distance function from $y$ to all the other points $x$ and the other condition follows by symmetry of $d$.
Is there a single, more symmetric PDE satisfied by $d(\cdot,\cdot)$ as a function on the product manifold $\mathcal M\times\mathcal M$?
That is, if the eikonal equation is the PDE behind the single-source-all-destinations geodesic distance problem, is there a different canonical PDE that governs the pairwise geodesic distance problem? I'm hoping to identify a condition that doesn't require enforcing the eikonal condition in the $x$ and $y$ coordinates individually.
[Cross post from Twitter:]
This is not rigorous, but it may help.
Consider the Euclidean case, where $d(x,y) = |x-y|$. Then
$$\nabla d = \frac{(x-y, y-x)}{|x-y|},$$
where $\nabla$ denotes the gradient on $\mathbb{R}^{2n}$. Hence,
$$|\nabla d|^2 = \frac{|x-y|^2 + |y-x|^2}{|x-y|^2} = 2.$$
So, a necessary condition is
$$ |\nabla d|^2 = 2, $$
subject to $d(x,x) = 0$.
In the Riemannian case, we likewise have
$$\nabla d = (\nabla_x d, \nabla_y d),$$
where $\nabla_x$ and $\nabla_y$ denote the gradient with respect to the first and second arguments, respectively. Hence,
$$|\nabla d|^2 =|\nabla_x d|^2+|\nabla_y d|^2 = 1+1 = 2.$$
So then, up to a factor 2, the distance function must satisfy the usual eikonal equation on the product space, subject to Dirichlet boundary conditions
$$d(x,x) = 0.$$
Since this equation has a unique (viscosity) solution, it should characterize the distance function.