The purpose of this answer is to provide a proof for the following result which Sergei mentioned in the accepted answer:

**Proposition**. Let $M$ be a complete Riemannian manifold and $x,y \in M\times M$, $x\neq y$. Then the following are equivalent:

The Riemannian distance $d:M\times M\rightarrow[0,\infty)$ is smooth in a neighbourhood of $(x,y)$.

There is only one length minimising geodesic connecting the points $x$ and $y$ and they are not conjugate along that geodesic.

*Proof of Proposition.* The Proposition follows from the three Lemmas below which freely use some properties of the so called
*segment domains* $\Sigma_x=\{w\in T_xM: d(x,\exp_x(w))=\vert w \vert\}$:

- $\exp_x: \mathrm{int} \Sigma_x\rightarrow M$ is a diffeomorphism onto its image [Gallot-Lafontaine, Corollary 3.77 or Petersen, Lemma 5.7.8 and Proposition 5.7.10]
- $M= \exp_x(\mathrm{int}\Sigma_x)\cup\exp_x(\partial \Sigma_x)$ and the union is disjoint [Gallot-Lafontaine, Proposition 2.113].
- Denote $\partial^1\Sigma_x = \{w\in \partial \Sigma_x: \exp_x(w)=\exp_x(w')$ for some $w'\in \partial \Sigma_x\backslash\{w\}\}$. Then:
- If $w\in \partial \Sigma_x\backslash \partial^1\Sigma_x$, then $D\exp_x\vert_w$ is singular. [Gallot-Lafontaine, Scholium 3.78 or Petersen, Lemma 5.78 ]
- $\exp_x(\partial^1\Sigma_x) \subset \exp_x(\partial \Sigma_x)$ is dense [Sakai, Remark 4.9].

**Warning.** The denseness of $\exp_x(\partial^1\Sigma_x) \subset \exp_x(\partial \Sigma_x)$ is crucial for the proof of Lemma 1, however Sakai does not prove this result in Remark 4.9, but rather states that *It is known that...*. To me the statement is not trivial, so I would appreciate if someone could explain this result or show how one can work without the density argument.

**Lemma 1.** $d^2(x,\cdot):M\rightarrow [0,\infty)$ is smooth in a neighbourhood of $y$ if and only if $y\in \exp_x(\mathrm{int} \Sigma_x).$

*Proof.* On (the open set) $\exp_x(\mathrm{int}\Sigma_x)$ we have $d(x,y) = \vert \exp_x^{-1}(y)\vert^2$ which is clearly a smooth function in $y$.
For the converse assume that $d(x,\cdot)^2$ is smooth on some open set $U\subset M$ and note that it suffices to prove
$$U \cap \exp_x(\partial \Sigma_x) = \emptyset \tag{1}.$$
Without loss of generality we may assume that $x\notin U$, then also $d(x,\cdot)$ is smooth in $U$ and has a gradient $G\in C^\infty(U;TM)$. Let $\gamma:[0,l]\rightarrow M$ be a length minimising unit-speed geodesic with $\gamma(0)=x$ and $y=\gamma(1)\in U$. Then
$d(x,\gamma(t))=t$ and differentiation yields $\langle G_y, \dot \gamma(l)\rangle = 1$. In particular $d(x,\cdot)$ has non-vanising gradient at $\gamma(l)$ and thus it is a submersion in a neighbourhood of $\gamma(l)$. This implies that the orthogonal complement of $\dot\gamma(l)\in T_yM$ is spanned by vectors $\dot c(0)$, where $c:(-\epsilon,\epsilon)\rightarrow M$ are curves with $c(0)=y$ and $d(x,c(s))=\mathrm{const}$. Since $\langle G_y , \dot c(0)\rangle = 0$ for such curves, we have
$
G_y = \dot \gamma(l).
$
We conclude:
$$\text{Length minimising geodesics which start in $x$ don't intersect in $U$.} \tag{2}$$
Now we can prove (1). Assume to the contrary that $U \cap \exp_x(\partial \Sigma_x) \neq \emptyset$. Then by densitity of $\exp_x(\partial^1 \Sigma_x)$ we also have $U \cap \exp_x(\partial^1 \Sigma_x) \neq \emptyset$ and there are $w,w'\in \partial \Sigma_x$ with $w\neq w'$ and $\exp_x(w) = \exp_x(w')\in U$, which is in contradiction to (2).

**Lemma 2.** $d^2:M\times M\rightarrow [0,\infty)$ is smooth in a neighbourhood of $(x,y)$ if and only if $y \in \exp_x(\mathrm{int} \Sigma_x)$.

*Proof.* If $d^2$ is smooth near $(x,y)$, then $d^2(x,\cdot)$ is smooth near $y$ and the previous Lemma implies that $y\in \exp_x(\mathrm{int}\Sigma_x)$. For the converse define $\Sigma = \bigcup_x \Sigma_x \subset TM$ and note that
$$
\Sigma \text{ is closed }\quad \text{ and } \quad \mathrm{int} \Sigma \cap T_xM = \mathrm{int} \Sigma_x \tag{3}.
$$
Define
$$
F: \mathrm{int} \Sigma \rightarrow M\times M, (x,w) \mapsto (x,\exp_x(w))
$$
and note that
$$
DF\vert_{(x,w)} = \begin{bmatrix}
\mathrm{id} & 0\\
\ast& D \exp_x\vert_w
\end{bmatrix}
$$
is invertible for all $(x,w)\in \mathrm{int} \Sigma$. Further $F$ is easily seen to be injective and thus it has a smooth inverse $F^{-1}:F(\mathrm{int} \Sigma) \rightarrow \mathrm{int} \Sigma$. Hence $d^2(x,y)= \vert F^{-1}(x,y)\vert ^2$ is smooth in a neighbourhood of every $(x,y) \in F(\mathrm{int} \Sigma)$, which concludes the proof.

**Lemma 3.** $y\in \exp_x(\mathrm{int} \Sigma_x)$ if and only if there exists a unique distance minimising geodesic between $x$ and $y$ and along this geodesic they are not conjugate.

*Proof.* Let $y=\exp_x(w)$ with $w \in \mathrm{int}\Sigma_x$. Then $t\mapsto \exp_x(tw)$, $0\le t\le 1$ is length minimising (because $w \in \Sigma_x$) and $x$ and $y$ are not conjugate along this geodesic ($D\exp_x\vert_w$ is invertible because $\exp_x$ is a diffeomorphism on $\mathrm{int}\Sigma_x$). If there was another length minimising geodesic from $x$ to $y$, then $y=\exp_x(w')$ for some $w'\in \Sigma_x \backslash \{w\}$. Since $\exp_x(\mathrm{int}\Sigma_x)\cap\exp_x(\partial \Sigma_x)=\emptyset$ we must have $w'\in \mathrm{int} \Sigma_x$, but this is false (since $\exp_x$ is injective on $\mathrm{int} \Sigma_x$).

Conversely assume that there is a unique distance minimising geodesic from $x$ to $y$ and that they are not conjugate along that geodesic. Then $y=\exp_x(w)$ for some $w\in \Sigma_x$. If we had $w\in \partial \Sigma_x$, then either there would be two length minimising geodesic between $x$ and $y$ (corresponding to $w\in \partial^1 \Sigma_x$) or $x$ and $y$ would be conjugate (corresponding to $D \exp_x\vert_w$ being singular).