The only rotation fields satisfying this PDE are constant $\newcommand{\div}{\operatorname{div}}$$\newcommand{\SO}{\operatorname{SO(2)}}$$\newcommand{\R}{\operatorname{\mathbb{R}}}$$\newcommand{\bdx}{\partial_x}$$\newcommand{\bdy}{\partial_y}$$\newcommand{\til}{\tilde}$
Let $\Omega \subseteq \R^2$ be an open connected domain, and let $U,V :\Omega \to \SO$ be smooth maps.

Claim: Suppose that $$ \div \bigg(U\begin{pmatrix} a & 0 \\\ 0 & b \end{pmatrix} V^T\bigg)=0 \tag{1}$$
for any $a,b \in \mathbb{R}$, where $\div$ is the standard divergence operator, acting row-by-row. Then $U,V$ are constant.

Comment: The assumption is equivalent to the following statement: For any  diagonal matrix $\begin{pmatrix} a & 0 \\\ 0 & b \end{pmatrix}$, the field $U\begin{pmatrix} a & 0 \\\ 0 & b \end{pmatrix}V^T$ equals $df$ for some smooth map $f:\Omega \to \mathbb{R}^2$.

I provide a proof below. I wonder whether there is a more conceptual or geometric proof. In particular, I don't "see a reason" for the appearance of complex analysis, which somehow arises in the proof.

Proof:
Set $$U=\begin{pmatrix} c & -s \\\ s & c \end{pmatrix}=R_{\theta}, V^T=\begin{pmatrix} \til c & - \til s \\\ \til s & \til c \end{pmatrix}=R_{\tilde \theta}.$$ Writing the system $(1)$ explicitly, we get
$$
 \begin{pmatrix} \bdx(c\til c) -\bdy(c\til s) & -\bdx(s\til s) -\bdy(\til c s) \\\ \bdx(s\til c)- \bdy(s \til s) & \bdx( \til s c) + \bdy(c \til c) \end{pmatrix} \begin{pmatrix} a   \\\ b \end{pmatrix}=0,
$$
so
$$
\begin{split}
&\partial_x(\tilde sc)=-\partial_y(c\tilde c) \\
 &\partial_y(\tilde s c)=\partial_x(c\tilde c) \\
\end{split} \tag{2}
$$
$$
\begin{split}
&\partial_x(s\tilde c)=\partial_y(s\tilde s) \\
&\partial_y(s\tilde c)=-\partial_x(s\tilde s).
\end{split} \tag{3}
$$
The systems $(2),(3)$ are equivalent to the holomorphicity of  $$f=-\tilde sc+i(c\tilde c)=i\cos(\theta)e^{i\tilde \theta},\,\,\,g=s\tilde c+i(s\tilde s)=\sin(\theta)e^{i\tilde \theta}.$$
Thus $(if)^2+g^2=e^{2i\tilde \theta}$ is holomorphic. Since its image lies on a circle, the open mapping theorem implies that it is constant, thus $e^{i\tilde \theta}$ is constant.
Together with the holomorphicity of $f,g$, we deduce that $c=\cos(\theta), s=\sin(\theta)$ are holomorphic, so are also constant.

Is there a simpler proof? Was this fact "obvious" -- could we somehow prove its without expanding $U,V$ explicitly in terms of $c,s,\tilde c,\tilde s$?
 A: NB: Evidently I used a different convention on divergence. In terms of matrix components my convention is $(\mathrm{div} M)_j = \sum_{i = 1}^2 \partial_{x^i} M_{ij}$. If your convention is $\sum_{i = 1}^2 \partial_{x^i} M_{ji}$, then everywhere I multiplied the divergence on the right by something, you should multiply on the left by the corresponding thing.

You don't even need all $a,b$. Let $I$ be the identity matrix and $R = \mathrm{diag}(1,-1) \in O(2)$ is the reflection matrix.
With the diagonal being $I$, you have
$$ \mathrm{div}(UV^T) = 0 $$
with the diagonal being $R$, you have
$$ \mathrm{div}(UR V^T) = 0 \iff \mathrm{div}(URV^T R) = 0 $$
the latter because you are just multiplying on the right with a constant, full rank, matrix so it commutes with taking the divergence.
Notice that $RV^TR = V$ by the conjugate action of the reflection on $SO(2)$.
Locally we can lift $U, V$ to mappings to the Lie algebra, and since $SO(2)$ is abelian we have that for some real valued functions $u,v$
$$ U = \exp(i u), \quad V = \exp(i v) $$
and that $UV^T = \exp(i(u-v))$ and $UV = \exp(i(u+v))$.
You have that
$$ 0 = \mathrm{div}(UV^T) VU^T = (id(u-v))^T $$
and similarly
$$ 0 = \mathrm{div}(UV) V^T U^T = (id(u+v))^T $$
adding and subtracting you find $du = dv = 0$.
