to $$\label{eik3}
1=\sum_{i}\lambda_{i}(x)d_{x_i}d_{x_i}$$ The latter equation is easier but there is no Taylor expansion around $y$ as can be seen in the case of constant coefficients (and hence constant eigenvalues $\lambda$), where the solution is $$d(x,y)=\sqrt{\sum_{i=1}^n \frac{\Delta x_i^2}{\lambda_i}}$$

We use equation mainly for the theoretical purposes of this corollary. In general it cannot be in general used for numerical purposes since this would imply that we have an efficient procedure to compute the eigenvalue functions of a space dependent matrix. Since we are looking for high precision in this paper, this is not possible in general. An exception is the case of dimension $n=2$ where we have $$\lambda_{1,2}(x)=\frac{\mbox{tr}(A)(x)}{2}\pm\sqrt{\left( \frac{\mbox{tr}(A)(x)}{2}\right) ^2-\mbox{det}(A)(x)}$$ where $A(x)=(a_{ij}(x))$.

Next we plug in the power series expansion $$d^2(x,y)=\sum_{i=1}^n\lambda^i_0\Delta x_i^2 +\sum_{|\beta |\geq 3 } d^2_{\beta}(y) \Delta x^{\beta}$$ We have $$d^2_{x_i}=2\lambda_0^i(y)\Delta x_i +\sum_{|\beta |\geq 3 } d^2_{\beta}(y) \beta_i \Delta x^{\beta \dot -1_i},$$ where for any multiindex $\beta$ we define $$\beta\dot-1_i=(\beta_1,\cdots,\beta_i,\cdots \beta_n)\dot-1_i:=\left\lbrace \begin{array}{ll} (\beta_1,\cdots,\beta_i-1,\cdots \beta_n) \mbox{ if } \beta_i\geq 1\\ (\beta_1,\cdots,0,\cdots \beta_n)~~\mbox{ else }  \end{array}\right.$$ The term $\beta -2_i$ is defined analogously. Plugging in the power series ansatz and using the relation $\lambda^0_i\left(  \lambda^i_0\right)^2=\lambda^0_i$, this leads to

$$\begin{array}{ll}
&\left( \sum_{|\beta |\geq 3} d^2_{\beta}(y)\Delta x^{\beta}\right) \left(1-\sum_i\beta_i \lambda_0^i\lambda^0_i\right) \\
\\
=&\left( \sum_{|\beta 