Let $(M^n,g)$ be a Riemannian manifold with a fixed point $p$. Can we find a local coordinate system $(x_1,x_2,\cdots, x_n)$ around $p$ satisfying the following two at the same time:
(1) $(x_1,x_2,\cdots, x_n)$ are harmonic coordinates, i.e., $\Delta_g x_i=0$;
(2) $(x_1,x_2,\cdots, x_n)$ are normal coordinates at $p$, i.e., $\partial_k g_{ij}=0$ at $p$?
Yes. This follows from a standard fact about the Laplacian $\Delta_g$, (because it is a second-order elliptic operator): The fact is this:
If $u$ is a smooth function on an open neighborhood $U$ of $p$ such that $\Delta_gu$ vanishes to order $k$ at $p$, then there is a (possiblly smaller) $p$-neighborhood $V\subseteq U$ on which there exists a smooth function $v$ that satisfies $\Delta_gv=0$ and $u-v$ vanishes to order $k{+}2$ at $p$.
Now fix $p$ and choose normal coordinates $y=(y^1,\ldots,y^n)$ centered at $p$. Because these are normal coordinates, $\displaystyle\frac{\partial g_{ij}}{\partial y^k}$ vanishes at $p$ for all $i,j,k$, and hence $\Delta_g y^i$ vanishes to order $1$ at $p$. Hence there exist functions $x^i$ on an open neighborhood of $p$ such that $\Delta_g x^i=0$ and $y^i-x^i$ vanishes to order $3$ at $p$ for all $1\le i\le n$.
By shrinking $V$, one can arrange that $x = (x^1,\ldots,x^n)$ is a coordinate system on a neighborhood of $p$ as well. Since $y^i = x^i + O(|x|^3)$, and since $$ g = g_{ij}(y)\,\mathrm{d}y^i\mathrm{d}y^j $$ where $g_{ij}(y)-g_{ij}(0)$ and $\mathrm{d}y^i - \mathrm{d}x^i$ all vanish to order at least $2$ at $p$, it follows that $g - g_{ij}(0)\,\mathrm{d}x^i\mathrm{d}x^j$ vanishes to order at least $2$ at $p$.
Hence, if $g = {\bar g}_{ij}(x)\,\mathrm{d}x^i\mathrm{d}x^j$, then $\displaystyle\frac{\partial {\bar g}_{ij}}{\partial x^k}$ vanishes at $p$ for all $i,j,k$, as desired.
