Existence and uniqueness of a quasi-linear pde system on a surface I have the following system of first order quasi-linear pde:
$$ -(\Delta+1) a^{\alpha\beta} [b_{\beta\rho} I_{\alpha;\sigma}+b_{\beta\sigma} I_{\alpha;\rho}]
 + a^{\alpha\beta}  [(\Delta+1) b_{\beta\rho;\sigma}+b_{\beta\rho}\Delta_{;\sigma}+b_{\beta\sigma}\Delta_{;\rho}]I_\alpha
 +a^{\alpha\beta} b_{\rho\sigma}I_\alpha I_\beta=0. \,\,\,\,\,\,\,(1)$$
Here, the subscript $(\cdot)_{;}$ denotes the covariant derivative with respect to the usual surface Christoffel symbols $\Gamma^\mu_{\alpha\beta}\equiv\frac{1}{2}a^{\mu\nu}(a_{\nu\alpha,\beta}+a_{\nu\beta,\alpha}-a_{\alpha\beta,\nu})$ on a non-compact surface $S$, embedded in $\mathbb{R}^3$, with first fundamental form $a_{\alpha\beta}$ and second fundamental form $b_{\alpha\beta}$. $\Delta(\ne -1)$ is a sufficiently smooth known function defined on $S$.
It can be observed that $I_\alpha=0$ clearly satisfies the above pde, and in my research problem, I want only this solution to occur. But I have to argue about the uniqueness of the above solution. Otherwise, if there are other solutions, I have to rule them out on `physical' grounds.
I wonder if anything could be said about the uniqueness of this pde, namely conditions on the coefficients, and, if possible, existence.
A brief background of the problem:
On an embedded smooth enough non-compact surface $s$ in $\mathbb{R}^3$, with first and second fundamental form $A_{\alpha\beta}$ and $B_{\alpha\beta}$ respectively, given are five real fields: $[E_{\alpha\beta}]$ (a symmetric positive definite matrix), $[\Lambda_{\alpha\beta}]$ (a non-necessary symmetric invertible matrix), ${\Lambda_\alpha}$, $\Delta_\alpha$ and $\Delta (\ne -1)$. We are required to construct, out of this information, a surface $S$ given by $\mathbf{r}(\theta^\alpha)$, with a transverse vector field $\mathbf{d}(\theta^\alpha)$ attached to every point on $S$, and having first fundamental form $a_{\alpha\beta}=A_{\alpha\beta}+2E_{\alpha\beta}$ and a suitably constructed second fundamental form $b_{\alpha\beta}$ out of the given fields, so that the following are true:
$$\Lambda_{\alpha\beta}=\mathbf{d}_{,\beta}\cdot\mathbf{a}_\alpha+B_{\alpha\beta}$$
$$\Lambda_\alpha=\mathbf{d}_{,\alpha}\cdot\mathbf{n}$$
$$\Delta_\alpha=\mathbf{d}\cdot\mathbf{a}_{\alpha}$$
$$\Delta=\mathbf{d}\cdot\mathbf{n}-1$$
Here, $\mathbf{n}$ is the unit normal field on $S$ and $\mathbf{a}_\alpha=\mathbf{r}_{,\alpha}$.
The necessary and sufficient conditions for the above to happen are given by
$$\Lambda_{[\alpha\beta]}-\Delta_{[\alpha;\beta]}=0,$$
$$I_\alpha\equiv \Lambda_\alpha-\Delta_{,\alpha}+\Delta_{\beta}a^{\beta\gamma}b_{\gamma\alpha}=0,$$
and the Gauss and Codazzi-Mainardi relations of $a_{\alpha\beta}$ and $b_{\alpha\beta}$ where $b_{\alpha\beta}\equiv\dfrac{\Lambda_{\gamma\alpha}-B_{\gamma\alpha}-\Delta_{\gamma;\alpha}}{\Delta+1}$.
Here, the subscript $(\cdot)_{[.]}$ denotes the skew part of a matrix. 
This has already been proved in M. Epstein, A Note on Nonlinear Compatibility Equations for Sandwich Shells and Cosserat Surfaces, Acta Mechanica 31, 285-289 (1979). In order to put this result in a different light, I arrived at equation (1).
 A: This is not a full solution, but it indicates the main points.  (A full solution would involve an EDS analysis that would be too long to include here; see below.)
First, setting $I_\alpha = (\Delta + 1)a_{\alpha\beta} J^\beta$ and noting that $\Delta+1$ is nonvanishing, the above equation becomes
$$ 
b_{\beta\rho} J^\beta_{;\sigma} + b_{\beta\sigma} J^\beta_{;\rho} = 
   b_{\beta\rho;\sigma}J^\beta +b_{\rho\sigma} \, a_{\alpha\beta} J^\alpha J^\beta .
$$
Note, by the way, that the Codazzi equation for $b$, the second fundamental form, is $b_{\beta\rho;\sigma}=b_{\beta\sigma;\rho}$, so both sides of the equation are symmetric in $\rho$ and $\sigma$, and this is three equations for two unknowns.
Second, as has already been remarked, if $b\equiv0$, then any vector field $J^\beta$ will furnish a solution.  The OP says that this case can be discarded, though.  
Third, in the 'generic' case that $b_{\rho\sigma}$ is a nondegenerate quadratic form, the situation is different.  In this case, this overdetermined system is of finite type.  In fact, the usual EDS analysis shows that knowing $J$ and its full first covariant derivative at one point determines $J$ uniquely in a neighborhood.  Thus, the space of local solutions on a connected open set has dimension at most $3$.  There may be conditions on the surface $S$ that will ensure that this local maximum is actually reached, but deriving those conditions will require differentiating $b$ at least two times.  For the generic surface, though, the prolonged system will not be compatible, and the only solution (even locally) will be $J=0$.  (This last statement may have some bearing on your problem, in the sense that, if the 'physical problem' you are considering is supposed to apply to the 'generic' surface, then, possibly, you can justify ignoring the nonzero solutions to your equation by simply saying that the nonzero solutions are not physical because they don't persist under arbitrarily small generic perturbations.) 
An EDS analysis that is too long to include here shows that, in fact, the congruence classes of real-analytic surfaces $S$ in $\mathbb{E}^3$ with nonvanishing Gauss curvature that admit at least one nonzero solution $J$ to the above equation depend on $4$ real-analytic functions of $1$ variable.  Thus, such surfaces do exist, at least locally, but it might be hard to write down any specific one since the equations are complicated.  They might simplify for surfaces of rotation, though, so it's not out of the question that explicit solutions could be found.
Fourth, when $b$ is a quadratic form of rank $1$, the induced metric is flat and the equations degenerate, so the solutions can be understood fairly simply.  In the generic case, in which the surface is not a generalized cylinder, i.e., if the principal curves of zero curvature (which are lines in 3-space) are not parallel, then it turns out that there are exactly two solutions to the above equations, one is $J=0$ and the other is a certain nonzero vector field tangent to the linear principal curves (see below).  In the special case that the surface is a generalized cylinder, the solutions depend on one arbitrary function of one variable.
Explicitly, if one diagonalizes the first and second fundamental forms in the form
$$
\mathrm{I} = {\omega_1}^2+{\omega_2}^2,
\qquad \mathrm{I\!I} = b_{11}\,{\omega_1}^2+0\,{\omega_2}^2 = b_{11}\,{\omega_1}^2\not=0,
$$
then one finds that $\mathrm{d}b_{11} = b_{11}(c_0\,\omega_1 -b_0\,\omega_2)$ for some functions $c_0$ and $b_0$, where $b_0$ vanishes if and only if the surface is a generalized cylinder.  The vector field $J$ is now described by the functions $J_i = \omega_i(J)$.
When $b_0$ is nonvanishing, a vector field $J$ on the surface satisfies the given equation if and only if either $(J_1,J_2)=(0,0)$ or $(J_1,J_2) = (0,3b_0)$.  Thus, there are two solutions.
When $b_0$ vanishes identically, one finds that there are coordinates $x$ and $y$ so that $\omega_1=\mathrm{d}x$ and $\omega_2=\mathrm{d}y$ and $J$ satisfies the given equation if and only if $J_1$ and $J_2$ are functions of $x$ that satisfy the single ODE
$$
J_1' = \tfrac12({J_1}^2+{J_2}^2 + c_0 J_1)
$$
(Of course, in this case, $c_0$ is also a function of $x$ alone.)
