Non-linear hyperbolic PDE I have the following PDE in two dimensions
$$
2\partial_x\partial_y\sqrt{1-u^2}+\left(\partial^2_x-\partial^2_y \right)u=0,
$$
with $u=u(x,y)$ with values between $-1$ and $1$, or alternatively
$$
2\partial_x\partial_y\sin2\theta(x,y)+\left(\partial^2_x-\partial^2_y \right)\cos2\theta(x,y)=0,
$$
with real $\theta(x,y)\sim\theta(x,y)+2\pi$, on some domain of the plane. Now, numerically I can obtain the solutions very quickly specifying some domain and an initial Cauchy line (as the equation hyperbolic), but I wish to have a deeper understanding of the solutions, so I'd like to see if there's a way to obtain analytic solutions. For instance, I know that $u=\cos(2\arctan(y/x))$ and $\theta(x,y)=\arctan(y/x)\pm1/2\arccos(c_1+c_2/(x^2+y^2))$, with $c_1, c_2$ some reals constants, are analytic, particular solutions, which strongly suggests that some general solution with arbitrary constants is plausible.
The problem is encountered in the context of elasticity of thin sheets. A so-called director field is imprinted on a thin elastic sheet, and it generates curvature upon a process called activation [1]. The director field $\theta(x,y)$ will induce a Riemannian metric on the new, curved sheet
$$
g(x,y)=R[\theta(x,y)]diag(\lambda_1,\lambda_2)R[\theta(x,y)]^T,
$$
with $R[\theta(x,y)]$ a $2\times2$ rotation matrix and $\lambda_1,\lambda_2$ some positive, known constants. Now, the aforementioned metric has a Gaussian curvature proportional to the equation written before, and the question I'm addressing is, for which $\theta(x,y)$s the generated curvature is zero ?, except for possibly isolated points where it may diverge. Now, the solutions I wrote before correspond to cones, but there should be more analytic solutions.
Any ideas ? Have you seen this equation or someone similar before ?
Thank you so much.
[1] Mostajeran, Cyrus; Warner, Mark; Ware, Taylor H.; White, Timothy J., Encoding Gaussian curvature in glassy and elastomeric liquid crystal solids, Proc. R. Soc. Lond., A, Math. Phys. Eng. Sci. 472, No. 2189, Article ID 20160112, 16 p. (2016). ZBL1371.82141.

Edit:
Here is a summary of the solution from Robert Bryant's great answer, in a language more familiar to physicist.
Consider the (in general muliply-)connected domain $\mathscr{{W}}\subseteq\mathbb{R}^{2}$, with Cartesian coordinates $u$ and $v$, and the function $f:\mathscr{{W}}\rightarrow\mathbb{R}$ that solves $\frac{\partial^{2}f(u,v)}{\partial u\partial v}=f(u,v)$,
with $f(u,v)$ and $\partial f(u,v)/\partial v$ non-vanishing. This is a linear hyperbolic equation, so it's Cauchy problem is always well defined on $\mathscr{{W}}$, and the space of solutions is always non-empty.
The 1-forms
$$
\alpha_{1}  \equiv f\cos\left(u-v\right)\mathrm{d} u+\frac{\partial f}{\partial v}\sin\left(u-v\right)\mathrm{d}v
,\:\:\alpha_{2} \equiv f\sin\left(u-v\right)\mathrm{d} u-\frac{\partial f}{\partial v}\cos\left(u-v\right)\mathrm{d}v,
$$
are closed. Therefore we can write locally $\mathrm{d} x    =\alpha_{1},\:
\mathrm{d} y    =\alpha_{2},$ for some functions $x$ and $y$ on $\mathscr{{W}}$. We define the function $(x,y):\mathscr{{W}}\rightarrow\mathbb{R}^{2}$ and the domain $\mathscr{{Z}\mathbb{\subseteq R}}^{2}$ as the image of $(x,y)$, i.e., $(x,y)\left(\mathscr{{W}}\right)=\mathscr{{Z}}$. We can use $x$ and $y$ as coordinates of $\mathscr{{Z}}$, as by definition they cover the latter completely. The function $u-v:\mathscr{{W}}\rightarrow\mathbb{R}$ can be pulled through $\mathscr{{Z}}$, that is $u-v=\theta\circ(x,y)$, with $\theta:\mathscr{{Z}}\rightarrow\mathbb{R}$ a function defined by the previous relation.
Inverting the definitions for $u$ and $v$ as functions of $x$ and $y$ we can write
$$\frac{\partial}{\partial x}   =\frac{1}{f}\cos\left(u-v\right)\frac{\partial}{\partial u}+\frac{1}{\frac{\partial f}{\partial v}}\sin\left(u-v\right)
\frac{\partial}{\partial v}\\
\frac{\partial}{\partial y} =\frac{1}{f}\sin\left(u-v\right)\frac{\partial}{\partial u}-\frac{1}{\frac{\partial f}{\partial v}}\cos\left(u-v\right)\frac{\partial}{\partial v},
$$
and it's just a matter of patience to verify that
$$
\begin{align}
&2\frac{\partial^{2}}{\partial x\partial y}\sin2\theta+\left(\frac{\partial^{2}}{\partial x^{2}}-\frac{\partial^{2}}{\partial y^{2}}\right)\cos2\theta\\
    &=\Bigg\{ \left[\frac{1}{f}\cos\left(u-v\right)\frac{\partial}{\partial u}+\frac{1}{\frac{\partial f}{\partial v}}\sin\left(u-v\right)\frac{\partial}{\partial v}\right]^{2}\\
&\:\:\:\:-\left[\frac{1}{f}\sin\left(u-v\right)\frac{\partial}{\partial u}-\frac{1}{\frac{\partial f}{\partial v}}\cos\left(u-v\right)\frac{\partial}{\partial v}\right]^{2}\Bigg\} \cos2\left(u-v\right)\\
    &+2\left[\frac{1}{f}\cos\left(u-v\right)\frac{\partial}{\partial u}+\frac{1}{\frac{\partial f}{\partial v}}\sin\left(u-v\right)\frac{\partial}{\partial v}\right]\\
&\:\:\:\times\left[\frac{1}{f}\sin\left(u-v\right)\frac{\partial}{\partial u}-\frac{1}{\frac{\partial f}{\partial v}}\cos\left(u-v\right)\frac{\partial}{\partial v}\right]\sin2\left(u-v\right)\\
    &=\frac{4\cos^{2}(u-v)}{f\left(\frac{\partial f}{\partial v}\right)^{2}}\left(f(u,v)-\frac{\partial f}{\partial u\partial v}\right)\\
    &=0.
\end{align}
$$
The particular solutions mentioned before are obtained with $f(u,v)=e^{\alpha u+v/\alpha}$, with some constant $\alpha$. But of course any other $f(u,v)$ will generate a solution. The most general real, separable solution is
$$ 
f(u,v)=\int_{-\infty}^{\infty}\mathrm{d}\rho\:C(\rho)\:e^{\rho u+\rho^{-1}v},
$$
for some arbitrary kernel $C(\rho)$. So one can classify arbitrarily many solutions by $C(\rho)$. From a calculative standpoint, once chosen some $f(u,v)$ the problem is that in general it's difficult to solve the algebraic system to write down explicitly $u(x,y)$ and $v(x,y)$, so one can say that the non-linear differential equation in two variables was transformed into a problem of two non-linear algebraic equations.
 A: 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.)
A: 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.
