Let $(M,g)$ be a smooth Riemannian manifold, and let $p \in M$.
Is there an open neighbourhood $U$ of $p$ that admits an orthonormal frame of divergence-free vector fields?
Edit: At least for dimension $2$, $M$ admits such a local frame if and only if it's flat (See a proof below). I am not sure about higher dimensions.
Edit 2:
It seems this question is addressed here. However, I don't understand the "counting" of the number of equations and variables:
Robert Bryant writes that this is a system of $n$ first-order PDE for $\tfrac12n(n{-}1)$ unknowns. I think that the following argument explains this:
The orthonormality conditions $\langle X_i,X_j \rangle=\delta_{ij}$ form $\tfrac12n(n{+}1)$ equations, and $\text{div}(X_i)=0$ form additional $n$ equations. So, together we have $\tfrac12n(n{+}1)+n$ equations in $n^2$ variables, when we express each $X_i$ in terms of some fixed arbitrary frame.
However, we can reduce the system by restricting the discussion in advance to orthonormal frames:
Indeed, given some fixe orthonormal frame $E_i$, we can represent any other orthonormal frame on $U$ as $QE_i$ where $Q:U \to \text{O}(n)$ is a smooth map. So, we have an "$ \text{O}(n)$" freedom to choose orthonormal frames; More explicitly, since $ \dim(\text{O}(n))=\tfrac12n(n{-}1)$, we can take $Q(p)=\text{Id}$ and for $q \in U$, $Q(q)$ can be expressed in terms of $\tfrac12n(n{-}1)$ functions (If $U$ is small enough so $Q(U)$ is contained in a coordinate chart around $\text{Id}$).
A proof for the $2$D case:
Suppose that $X,Y$ are orthonormal and divergence-free. Then $$ \langle X,X \rangle=1 \Rightarrow \langle \nabla_YX,X \rangle=\langle \nabla_XX,X \rangle=0. \tag{1}$$
The divergence-free condition means that $$ 0=\text{div}X=\text{trace}(\nabla X)=\langle \nabla_XX,X \rangle+\langle \nabla_YX,Y \rangle=\langle \nabla_YX,Y \rangle. \tag{2}$$
Combining equations $(1),(2)$ we see that $$ \langle \nabla_YX,X \rangle=\langle \nabla_YX,Y \rangle=0$$
so $\nabla_YX=0$. By symmetry, we also have $\nabla_XY=0$, so the symmetry of the connection implies $[X,Y]=0$. This in turn implies $X,Y$ can be realized as coordinate vector fields, but since they are orthonormal this means the metric is flat.
Alternatively, we can proceed from $\nabla_YX=\nabla_XY=0$ as follows:
Differentiating $\langle X,Y \rangle=0$ via $X$, we get
$$ 0=\langle \nabla_XX,Y \rangle+\langle X,\nabla_XY \rangle=\langle \nabla_XX,Y \rangle. \tag{3}$$
Combining this with $\langle \nabla_XX,X \rangle=0$ (see equation $(1)$ again) we deduce $\nabla_XX=0$, which together with $\nabla_YX=0$ implies $X$ is parallel. By symmetry, $Y$ is also parallel, so we have a parallel frame for $(TM,\nabla)$ which implies $g$ is flat.
Comment: We always have a local orthonormal frame;
Furthermore, there are always divergence-free frames: Indeed, every volume form can be locally written as $dx^1 \wedge \dots \wedge dx^n$, for some coordinates $x_i$. The divergence w.r.t this form is the standard one, i.e. if $V=v^i\partial_i$, then $\text{div}V=\partial_i v^i$, so in particular the coordinate frame $\partial_i$ form a divergence-free frame.
We can apply the Gram-Schmidt process on $\partial_i$, but I see now reason why the "divergence-free" property should be preserved.
By divergence of a vector field $X$, I refer to the Riemannian notion:
$\text{div} X= \text{trace}(\nabla X)$, where $\nabla$ is the Levi-Civita connection of $(M,g)$. Alternatively, $\text{div} X=0 \iff L_X\text{Vol}_g=0$ where $\text{Vol}_g$ is the Riemannian volume form of $(M,g)$.