Existence of orthogonal coordinates on a Riemannian manifold

This is probably a very naive question, but so far I could not find an answer:

Let $(M,g)$ be a Riemannian manifold. Can we always find "orthogonal coordinates" locally?

More precisely, I am asking if for every $p \in M$ there exists a neighbourhood $U$ and a diffeomorphism $\phi:\mathbb{R}^n \to U$, such that $g_{ij}=g(d\phi(e_i),d\phi(e_j))=0$ for $i \neq j$.

Clarification: Note that I want $g_{ij}=0$ on all $U$, not just at $p$. Also, I allow $g_{ii} \neq g_{jj}$ for $i \neq j$ (the special case where $g_{ii}$ is independent of $i$ is called isothermal coordinates-and corresponds to conformal flatness of $U$).

Of course, this is weaker than requiring $M$ to be conformally flat, since a (linear) map which maps an orthogonal basis to an orthogonal basis does not need to be conformal.

A Riemannian metric $g$ on an $n$-dimensional manifold is called locally diagonalizable if it is locally isometric to a Riemannian metric on a domain in $R^n$ with diagonal metric tensor. In dimension $n=2$ every Riemannian metric is locally diagonalizable due to existence of isothermal coordinates.

For $n\ge 3$ the problem of local diagonalizability was solved in

D. DeTurck, D. Yang, Existence of elastic deformations with prescribed principal strains and triply orthogonal systems. Duke Math. J. 51 (1984), no. 2, 243–260.

They proved that for $n=3$ every Riemannian metric is indeed locally diagonalizable while for all $n\ge 4$ there are obstructions to local diagonalizability. For instance, $W(e_i, e_j, e_k, e_l)=0$ for every orthonormal frame with distinct $i, j, k, l$, here $W$ is the Weyl tensor.