This confirms uniqueness. (Note that there is no loss of generality if we choose the normal coordinates for the second order terms). In general the solution is not globally analytic in the sense that $d^2$ is not representable by a globally converging power series.

A smooth solution $d^2$ of the eikonal equation has the representation $$\begin{array}{ll}
d^2(x,y)&=d(y,y) + \nabla d(y,y)\cdot (x-y)\\
\\
&+\sum_{|\gamma|=2}\int_0^1(1-\theta)^{1}\frac{\Delta x^{\gamma}}{\gamma !}\partial^{\gamma}d^2(x+\theta \Delta x,y)d\theta.
\end{array}$$ We abbreviate $R(x,y)=\sum_{|\gamma|=2}\int_0^1(1-\theta)^{1}\frac{\Delta x^{\gamma}}{\gamma !}\partial^{\gamma}d^2(x+\theta \Delta x,y)d\theta.$ Since $d(y,y)=0$ we have $$\begin{array}{ll}
d^2(x,y)=\nabla d^2(y,y)\cdot (x-y)+R(x,y)
\end{array}$$ The ’only if’-condition of the boundary condition leads to $\nabla d^2(y,y)=0$. To see this assume that $\nabla d^2(y,y)\neq 0$. Since $R(x,y)\leq C\|\Delta x\|^2$ there is a small $\Delta x$ such that $\nabla d^2(y,y)\cdot \mu\Delta x >C\|\Delta x\|^2$ and $\nabla d^2(y,y)\cdot (-\mu)\Delta x <-C\|\Delta x\|^2$ for some $\mu\in (0,1]$. Hence there exists some $\rho$ such that with $x':=y+(\rho\mu)\Delta x$ $$d^2(x',y)=\nabla d^2(y,y)\cdot (\rho\mu)\Delta x + R(x',y)=0,$$ contradicting one part of the boundary condition $d^2(x,y)=0~\mbox{ iff }~x=y$. Next one computes that $\sum_{ij}a^{ij}(y)\Delta x^i\Delta x^j$ satisfies the equation $$d^2(x,y)=\frac{1}{4}\sum_{ij}a_{ij}(y)d^2_{x_i}d^2_{x_j},$$ and the uniqueness of theorem 2.3. (which we established by arguing with uniqueness of related diffusions and Varadhan’s result, in the Atiyah-Singer spirit of short-range analytic expansions) identifies the coefficients $a^{ij}(y)$ 