Non-linear hyperbolic PDE

As I understand it, the equation you are imposing on the function $\theta(x,y)$, defined on a domain $D\subset\mathbb{R}^2$ in the $xy$-plane is that, for some positive constants $\lambda_1\not=\lambda_2$, the metric $$ g = \lambda_1\,(\cos\theta(x,y)\,\mathrm{d}x+\sin\theta(x,y)\,\mathrm{d}y)^2 + \lambda_2\,(\sin\theta(x,y)\,\mathrm{d}x-\cos\theta(x,y)\,\mathrm{d}y)^2 $$ should be flat, i.e., that there should exist functions $p,q:D\to\mathbb{R}^2$ so that $g = \mathrm{d}p^2 + \mathrm{d}q^2$. In other words, the Jacobian matrix of the mapping $f = (p,q):D\to\mathbb{R}^2$ should have distinct, constant singular values.

Now, this is exactly (a local version of) the question posed in Are all maps $\mathbb{R}^2 \to \mathbb{R}^2$ with fixed singular values affine? In my answer to that question, I showed that any sufficiently differentiable solution $f$ defined on all of $D = \mathbb{R}^2$ must, in fact, be affine (equivalently, that $\theta$ must be constant), and I gave a formula for the local solutions that satisfy a non-degeneracy condition that shows how to reduce this problem locally to the hyperbolic linear equation $f_{uv} = f$ on an auxilliary domain $D'$ in the $uv$-plane plus a couple of 'quadratures' (i.e., writing an explicit closed $1$-form as the differential of a function).

I won't reproduce the analysis here, I'll just give the solution described there: Let $D'$ be a simply-connected domain in the $uv$-plane and let $f$ be a function on $D'$ such that $f_{uv} = f$ while $f$ and $f_v$ are nonvanishing. (It's easy to write many explicit solutions to this equation by separation of variables.) One easily sees that, when $f_{uv}=f$, the $1$-forms $$ \begin{aligned} \alpha_1 &= \cos(u{-}v)\,f\,\mathrm{d}u +\sin(u{-}v)\,f_v\,\mathrm{d}v\\ \alpha_2 &= \sin(u{-}v)\,f\,\mathrm{d}u -\cos(u{-}v)\,f_v\,\mathrm{d}v \end{aligned} $$ are closed, and hence one can write $\alpha_1 = \mathrm{d}x$ and $\alpha_2 = \mathrm{d}y$ for some functions $x$ and $y$ on $D'$. Suppose that $(x,y):D'\to \mathbb{R}^2$ is one-to-one and, hence, a diffeomorphism, and let $\theta:(x,y)(D')=D\to\mathbb{R}$ be the function that satisfies $\theta(x,y) = u-v$. Then $\theta$ satisfies the given equation, as follows from the Chain Rule. This is the 'general' local smooth solution, i.e., every smooth solution for which $\mathrm{d}(\cos\theta(x,y)\,\mathrm{d}x+\sin\theta(x,y)\,\mathrm{d}y)$ and $\mathrm{d}(\sin\theta(x,y)\,\mathrm{d}x-\cos\theta(x,y)\,\mathrm{d}y)$ are nonvanishing can be written in the above form locally. (Of course, if one imposes the condition that one of these $2$-forms vanishes identically, the problem is even more easily solved, and one can see the formula for those solutions in the solution that I referenced above.)

Write $z = e^{i2\theta}$ where $\theta$ is as in your second formulation, you have that the equation is equivalent to

$$ -2i \partial^2_{xy} (z - \bar{z})+ (\partial^2_{xx} - \partial^2_{yy})(z + \bar{z}) = 0 $$

which you can rewrite as

$$ \Re (\partial^2_{xx} - 2i \partial^2_{xy} + i^2 \partial^2_{yy}) z = 0 $$

factoring you get

$$ \Re (\partial_x - i \partial_y)^2 z = 0 \tag{*}$$

with the constraint that $|z| = 1$. A special class of (*) would be when $(\partial_x - i \partial_y)^2 z = 0$. It is impossible for a non-trivial solution to satisfy $(\partial_x - i \partial_y)z \equiv 0$: this would mean $z$ is (anti)holomorphic and the modulus constraint would imply that $z$ is constant. However, it is possible for $z$ to be linear in $x + i y$. To maintain the constant modulus condition this would require $$ z(x,y) = \frac{x + i y + c}{x - i y + \bar{c}} \tag{**} $$ for some fixed complex number $c$.

(**) includes the particular solutions that you mentioned and is special in that it comes from solving (completely) the stronger equation $(\partial_x - i \partial_y)^2 z = 0$. Note that the original equation (*) reads essentially $(\partial_x - i \partial_y)^2 z \in i \mathbb{R}$. A slight generalization of the special solution can be taken by postulating that $z$ takes the form

$$ z = \frac{w}{\bar{w}} $$

with some holomorphic function $w$. Applying the equation we find this requires

$$ w'' \cdot w \in i \mathbb{R}$$

and hence $w$ must solve the (complex) ODE

$$ w w'' = \lambda i, \quad \lambda\in\mathbb{R}. $$

Near a point this can be solved using power series methods, and Wolfram Alpha only returns a function that is expressible in terms of the inverse error function, so I don't think general solutions can be found using very nice closed-form formulae in terms of elementary functions.