Riemannian surfaces with an explicit distance function?

I'll briefly spell out what others have pointed to concerning geodesics on surfaces of revolution (or more generally, surfaces with a 1-parameter group of symmetries), because it's nice and not as widely understood as it should be.

Geodesics on surfaces of revolution conserve angular momentum about the central axis, so the geodesic flow splits into 2-dimensional surfaces having constant energy (~length) and angular momentum (The more general principle is that the inner product of the tangent to a geodesic with any infinitesimal isometry of a Riemannian manifold is constant). The surfaces are generically toruses. The shadow of these toruses on the surface of revolution is an annulus, a component of a set of $r \ge r_0$, where on each point with $r > r_0$ there are two vectors having the given angular momentum, but they merge at the boundary, both becoming tangent to the boundary of the annulus. If you sketch the picture, you will see the torus. The geodesics correspond to the physical phenomenon of the pattern of string or thread mechnically but passively wound around a cylinder. As string builds up in the middle, geodesics start to oscillate back and forth in a sinusoidal pattern, further amplifying the bulge in the middle.

To find the geodesic from point x to point y, you need to know which angular momentum will take you from x to y. For any two meridian circles and any choice of angular momentum, the geodesics of given angular momentum map one circle to the other by a rotation. Both the angle of rotation of the map and the length of the particular family of geodesics traversing the annulus is given by an integral over an interval cutting across the annulus, since the the slope of the vector field at all intervening points is known. I have an aversion to actual symbolic computation so I won't give you example formulas, but I believe this should meet your criterion for explicitness.

But to take a step back: this question, asking for an explicit formula, has an unstated (and probably unintended) connotation that is worth examining: this use of language implicitly suggests that non-symbolic forms are less worthy. I don't know the background motivation for the question, but an alternative question for some purposes would be to give example of surfaces where you can exhibit the distance function. Communication of mathematics is biassed toward symbolic forms. However, for many people and many purposes, some kind of graphical representation of the distance function, and/or diagrams or explanations of why it is what it is as well as a striaghtforward method for computing it, would often be better than a symbolic answer.

The geodesic flow of course is an ordinary differential equation. It is a vector field on the 3-manifold of unit-length tangent vectors to the surface, defined by very easy equations: the vectors are tangent to the surface, and their derivative (= the 2nd derivative of a geodesic arc) is normal to the surface. The solutions may not always have a nice symbolic form, but they always have a nice and easy-to-compute geometric form. Finding the distance involves the implicit function theorem, but this is easy and intuitive. One could, for instance, easily draw a parametric surface that is the graph of distance as a function of position directly from solutions to the ODE (which no doubt sometimes even have reasonable symbolic representations). Both the ODE for the geodesic flow and the inverse function to give distance as a function of position are easy to compute numerically, and easy to understand qualitatively.

NB (3/1/13): I revised this answer to make it more complete (and, to be frank, more accurate). My original answer did not take into account the difference between the cut locus and the conjugate locus, and, of course, this affects the formula for the distance between points.

I'm aware of a few metrics with non-constant curvature for which one can write the distance function explicitly in terms of the coordinates. The simplest such metric I know is the (incomplete) metric $ds^2 = y\ (dx^2+dy^2)$ on the upper half plane $y>0$. The Gauss curvature of this metric is $K = 1/(2y^3)>0$, so it's not constant.

Every geodesic of this metric in the upper half plane can be parametrized in the form $$ x = a + b\ t\qquad\qquad y = b^2 + \frac{t^2}{4} $$ for some constants $a$ and $b$, and, for such a geodesic, the arclength function along the curve is $$ s = c + b^2\ t + \frac{t^3}{12}\ . $$ for some constant $c$.

Using these formulae, one finds that two points $(x_1,y_1)$ and $(x_2,y_2)$ are joinable by a geodesic segment if and only if $4y_1y_2 \ge (x_1{-}x_2)^2$. In the case of strict inequality, there are two geodesic segments joining the two points, and the length of the shorter segment is $$ L_1\bigl((x_1,y_1),(x_2,y_2)\bigr) = {1\over3}\sqrt{3(x_1{-}x_2)^2(y_1{+}y_2)+4(y_1^3{+}y_2^3) - (4y_1y_2-(x_1{-}x_2)^2)^{3/2}}\ . $$ Note that, in a sense, this is better than the constant curvature case. Here, the distance function is algebraic in suitable coordinates, whereas, in the constant nonzero curvature cases, the distance function is not.

However, the function $L_1$ does not necessarily give the actual distance between the two points (i.e., the infinimum of the lengths of curves joining the two points), and it's not only because not every pair of points can be joined by a geodesic. To see this, one should complete the upper half plane by adding a point that represents the 'boundary' $y=0$. The Riemannian metric does not extend smoothly across this 'point', of course (after all, the Gauss curvature blows up at you approach this point), but it does extend as a metric space. The vertical lines, which are geodesics, can then be used to join $(x_1,y_1)$ to $(x_2,y_2)$ by going through the singular point, and the total length of this geodesic is $$ L_2\bigl((x_1,y_1),(x_2,y_2)\bigr) = \frac{2}{3}\bigl({y_1}^{3/2}+{y_2}^{3/2}\bigr). $$ (Also, note that $L_2$ is defined for any pair of points in the upper half-plane.) If one doesn't like this path that goes through the singular point, one can easily perturb it slightly to avoid the singular point and not increase the length by much, so it's clear that the infimum of lengths of curves lying strictly in the upper half plane and joining the two points is no more than $L_2$.

This suggests that the true distance function $L$ should be the minimum of $L_1$ and $L_2$ where they are both defined, i.e., where $4y_1y_2 \ge (x_1{-}x_2)^2$, and $L_2$ on the set where $4y_1y_2 < (x_1{-}x_2)^2$.

To get a sense of how these two formulae interact, one can use the fact that $x$-translation preserves the metric while the scalings $(x,y)\mapsto (ax,ay)$ for $a>0$ preserve the metric up to a homothety (and hence preserve the geodesics and scale the distances). These two actions generate a transitive group on the upper half plane, so, it suffices to see how these two functions interact when $(x_1,y_1) = (0,1)$, i.e., to see the conjugate locus and cut locus of this point.

The conjugate locus is easy: It's just $y-x^2/4=0$, which is the boundary of the region $y-x^2/4\ge0$ consisting of the points that can be joined to $(0,1)$ by a geodesic segment. Meanwhile, the cut locus is given by points $(x,y)$ that satisfy $y-x^2/4\ge0$ and for which $L_1\bigl((0,1),(x,y)\bigr) = L_2\bigl((0,1),(x,y)\bigr)$. In fact, one has $L_1\bigl((0,1),(x,y)\bigr) < L_2\bigl((0,1),(x,y)\bigr)$ only when $y > f(x)$, where $f$ is a certain even algebraic function of $x$ that satisfies $f(x) \ge x^2/4$ (with equality only when $x=0$). Moreover, for $|x|$ small, one has $$ f(x) = \left({\frac{{\sqrt{3}}}{4}}x\right)^{4/3} + O(x^2) $$ while, for $|x|$ large, one has $$ f(x) = \left({\frac{\sqrt{3}}{4}}x\right)^{4} + o(x^4). $$

Thus, all of the geodesics leaving $(x,y)=(0,1)$, other than the vertical ones, meet the cut locus before they reach the conjugate locus (and they all do meet the conjugate locus).

Thus, the actual distance function for this metric is explicit (it's essentially the minimum of $L_1$ and $L_2$), but it is only semi-algebraic.

Remark [by Matt F]: The following graph shows the contour lines for distances from $(0,1)$. The conjugate locus is in white, and the cut locus goes through the corners in the contour lines.

Remark: The thing that makes this work is that, while the metric has only a 1-parameter family of symmetries, it has a 2-parameter family of homotheties (as described above), and this extra symmetry of the geodesics is critical for making this work. Of course, there are other such metrics, all the ones of the form $ds^2 = y^{a}\ (dx^2+dy^2)$ ($a$ is a constant) have this property and don't have constant curvature unless $a = 0$ or $a = -2$. You don't get algebraic answers for all values of $a$, of course, but there is a way to get $D$ implicitly defined in terms of a special function (depending on the value of $a$).

More generally, the metrics whose geodesics admit more symmetries than the metric itself does tend to have such formulae. I'm not aware of any other cases in which one can get $D$ so explicitly.

In the course of writing an answer to a related MO question, I realized that there is a surface with a complete Riemannian metric of non-constant negative curvature for which one can write down the distance function explicitly, so I thought I would record it here for those who might be interested.

Such metrics are quite rare; even when the geodesic flow is integrable (or even rotationally symmetric), one cannot generally compute the arc length along geodesics in a sufficiently explicit form that one can actually compute the geodesic distance between two given points in any explicit way. This is the first complete example with nonconstant curvature that I have seen. (There are many explicit but non-complete examples with non-constant curvature in the classical literature, c.f. Tome III of Darboux' monumental Leçons sur la théorie générale des surfaces et les applications géométriques du calcul infinitésimal.)

The surface is $\mathbb{R}^2$ and the metric in standard coordinates is the rotationally symmetric metric $$ g = (x^2+y^2+2)\,(\mathrm{d}x^2 + \mathrm{d}y^2). $$ The Gauss curvature of $g$ is $K = -4/(x^2+y^2+2)^3<0$. It is complete, since it dominates the standard flat metric. It follows from general theory that any two points lie on a unique geodesic and each geodesic segment minimizes $g$-distance between its endpoints.

The geodesics of $g$ are easy to describe as curves: For every pair of numbers $(a,b)$ with $a^2+b^2\ge 1$, consider the equation $$ (1+a)\,x^2 + 2b\,xy + (1-a)\,y^2 = a^2+b^2-1. $$ When $a^2+b^2>1$, this is a hyperbola, and each of the branches is a geodesic. When $a^2+b^2=1$, this is the equation of a line through the origin, which is also a geodesic. Conversely, every geodesic of $g$ is either a line through the origin or a branch of one of the hyperbolae listed above.

The geodesic distance along a line through the origin is not hard to write down: On the line $y=0$, the element of arc length is $$ \mathrm{d}s = \sqrt{x^2+2}\,\mathrm{d}x = \mathrm{d}\left(\sinh^{-1}\left(\frac{x}{\sqrt{2}}\right)+\frac{x\sqrt{x^2+2}}{2}\right). $$ Set $$ f(x) = \sinh^{-1}\left(\frac{x}{\sqrt{2}}\right)+\frac{x\sqrt{x^2+2}}{2} \approx \sqrt2\left(x + \frac{x^3}{12}-\frac{x^5}{160}+\cdots\right). $$

I will now show that the $g$-distance between any two points $p,q\in\mathbb{R}^2$ is given by the formula $$ \delta(p,q) = f\left(\frac{|p+q|+|p-q|}{2}\right)-f\left(\frac{|p+q|-|p-q|}{2}\right), $$ where the norms are the Euclidean norms, i.e., taken with respect to the standard Euclidean inner product on $\mathbb{R}^2$.

To see this, first, note that, while the distance function $\delta:\mathbb{R}^2\times\mathbb{R}^2\to\mathbb{R}$ is not smooth along the diagonal, its square $\sigma = \delta^2$ is a smooth function on $\mathbb{R}^2\times\mathbb{R}^2$ that vanishes along the diagonal. In fact, because $g$ is real-analytic, it follows that $\sigma$ is real-analytic. Because $g$ is invariant under (Euclidean) rotation about the origin and reflection across lines through the origin, it follows that $\delta$ and $\sigma$ are also invariant under these rotations and reflections, now acting diagonally on $\mathbb{R}^2\times\mathbb{R}^2$. Using this, one can show that $\sigma$ must be representable as $$ \sigma(p,q) = C\bigl(|p|^2,\,p{\cdot}q,\,|q|^2\bigr) \quad\text{for all}\ p,q\in\mathbb{R}^2, $$ where $C(a,b,c)$ is a smooth function on the cone $\mathcal{C}_+$ defined by $a,c\ge 0$ and $ac-b^2\ge0$.

Now, for fixed $q\in\mathbb{R}^2$ the function $\delta_q:\mathbb{R}^2\to\mathbb{R}$, defined by $\delta_q(p) = \delta(p,q)$, vanishes at $q$ and satisfies $|\mathrm{d}(\delta_q)|^2_g = 1$ except at $q$ (where it is not differentiable). This implies that the corresponding $\sigma_q = {\delta_q}^2$ attains its minimum value of $0$ at $q$ and satisfies the first-order PDE $|\mathrm{d}(\sigma_q)|^2_g = 4\sigma_q$. Interpreting this in terms of the above representation of $\sigma$, we find that $C$ must satisfy the first order PDE $$ 4aC_a^2 + 4bC_aC_b+cC_b^2 - 4(a+2)C^2 = 0. $$ Similarly, using the fact that $C(a,b,c) = C(c,b,a)$ (since $\sigma(p,q) = \sigma(q,p)$), we find that $$ 4cC_c^2 + 4bC_cC_b+aC_b^2 - 4(c+2)C^2 = 0. $$ This pair of first-order PDE for $C$ is singular at $(a,b,c) = (0,0,0)$, but, since $C$ must vanish when $a+c-2b = |p-q|^2 = 0$ but otherwise be positive in the cone $\mathcal{C}_+$, it is easy to show that $C$ has a Taylor expansion $$ C\simeq (a{-}2b{+}c)\left(2 + \frac{(a{+}b{+}c)}{3}-\frac{(4a{+}7b{+}4c)(a{-}2b{+}c)}{360} + \cdots\right). $$ In fact, examining the higher terms, it becomes apparent that $C$ should be a function of $u = a{+}c$ and $v = a{-}2b{+}c$. Indeed, if $$ C(a,b,c) = H(a{+}c,\,a{-}2b{+}c) = H(u,v) $$ were to hold for some smooth function $H$ on the $uv$-domain defined by $0\le v\le 2u$, then one finds that $H$ would have to satisfy $$ u\,{H_u}^2 + 2v\,(H_uH_v+{H_v}^2) - (u+4)\,H = 0. $$ with $H \simeq v\,\bigr(2-\tfrac1{6}(v-3u)-\tfrac1{720}(15u-7v)v+\cdots\bigr)$. Using the theory of singular analytic first-order PDE, it is not difficult to show that such an analytic solution $H(u,v)$ exists, is unique, and is a multiple of $v$. (It's easy to show that there is a unique power series solution whose lowest term is $2v$, but one needs to show that this series converges.)

As a consequence, $C(a,b,c) = H(a{+}c,\,a{-}2b{+}c)$ satisfies the pair of first-order singular analytic PDE listed above. Consequently, $$ \sigma(p,q) = C\bigl(|p|^2,\,p{\cdot}q,\,|q|^2\bigr) = H\bigl(|p|^2{+}|q|^2,\,|p{-}q|^2\bigr), $$ Since $H$ is a multiple of $v = |p{-}q|^2$, it follows that $$ \delta(p,q) = |p{-}q|\,G\bigl(|p|^2{+}|q|^2,\,|p{-}q|^2\bigr) $$ for some smooth positive function $G(u,v)$. Meanwhile, for $b<a\in\mathbb{R}$, taking $p = (a,0)$ and $q=(b,0)$, we have $$ (a{-}b)\,G\bigl(a^2{+}b^2,\,(a{-}b)^2\bigr) = \delta(p,q) = f(a)-f(b). $$ The function $G(u,v)$ is determined in the wedge $0\le v\le 2u$ by this equation as $(a,b)$ vary over the half-plane $b<a$. It follows from this that $$ |p{-}q|\,G\bigl(|p|^2{+}|q|^2,\,|p{-}q|^2\bigr) = f\left(\frac{|p+q|+|p-q|}{2}\right)-f\left(\frac{|p+q|-|p-q|}{2}\right), $$ as desired.

Remark: The reader might be startled (as I was initially) upon realizing that the above formula implies a seemingly strange identity $$ f\left(\frac{|a+b|+|a-b|}{2}\right)-f\left(\frac{|a+b|-|a-b|}{2}\right) = |f(a)-f(b)| $$ for all real numbers $a$ and $b$, but, in fact, this identity holds for any increasing odd function $f$.

Added Remark (16 May 2020): A similar analysis, yielding an explicit distance function, can be done for the incomplete metric $$ g = (1-x^2-y^2)\bigl(\mathrm{d}x^2+\mathrm{d}y^2\bigr) $$ on the interior of the unit disk $D$ defined by $x^2+y^2<1$. This is a metric of positive curvature $K = 4/(1-x^2-y^2)^3$. What one finds is that, setting $$ s(x) = \tfrac12\arcsin(x) + \tfrac12x\sqrt{1-x^2} \quad\text{for}\ |x|\le 1, $$ the function $$ \delta(p,q) = s\left(\frac{|p+q|+|p-q|}{2}\right)-s\left(\frac{|p+q|-|p-q|}{2}\right) $$ gives the length of the shortest geodesic joining $p$ and $q$ when $|p+q|+|p-q|\le 2$. (This inequality is also the condition for the existence of a geodesic joining $p$ and $q$ within the interior of $D$.)

Meanwhile, regarding the boundary circle $x^2+y^2=1$ as a single point $z$ whose distance from $p\in D$ is $s(1) - s(|p|) = \tfrac14\pi - s(|p|)$, we see that there is always a path from $p$ to $q$ (through z) of length $L(p,q) = \tfrac12\pi - s(|p|)-s(|q|)$.

It is now not difficult to show that the actual distance from $p$ to $q$ is $L(p,q)$ when $|p+q|+|p-q|\ge 2$ and is the minimum of $\delta(p,q)$ and $L(p,q)$ when $|p+q|+|p-q|\le 2$.