I am trying to understand some derivations surrounding the Riemann Invariants for the system: $$ \begin{pmatrix} u_t \\ \eta_t \end{pmatrix} + \begin{pmatrix} u & 1 \\ \eta & u \end{pmatrix} \begin{pmatrix} u_x \\ \eta_x \end{pmatrix} = 0$$
So far, all I have self-derived is that the eigenvalues are $\lambda_{1,2}=u\pm \sqrt{\eta}$ with corresponding eigenvectors:$ \begin{pmatrix} 1 \\ \pm \sqrt{\eta} \end{pmatrix}$
I know that the Riemann invariants are $R_{1,2}(\eta,u)=u\pm 2\sqrt{\eta}$ and I noticed I can arrive at the same invariants from another stackexchange answer which stated: "The corresponding Riemann invariants $R_\mp$ have their gradient orthogonal to the eigenvectors". Using this info I found that: $$ \begin{pmatrix} 1 \\ \pm \sqrt{\eta} \end{pmatrix} \cdot\nabla_{\bf U} R_\mp = 0 ,\qquad\text{e.g.}\qquad \nabla_{\bf U} R_\mp = \begin{pmatrix} 1\\ \mp \dfrac{1}{\sqrt{\eta}}\end{pmatrix} . $$
Where ${\bf U} = (u,\eta)^\top$. Which, when I solved via integration (excluding the arbitrary constant), gave me the correct Riemann invariants above. However, eigenvectors are determined up to a multiplicative constant. This is also true for the orthogonal vector, and if I replace $\begin{pmatrix} 1\\ \mp \dfrac{1}{\sqrt{\eta}}\end{pmatrix}$ with the equally orthogonal $\begin{pmatrix} \eta\\ \mp \sqrt{\eta}\end{pmatrix}$ (or the various other products I can make), I get entirely different values for the invariant via integration!
Am I doing this correctly? If so, how do I justify the issue I demonstrated? I feel pretty lost here, so any guidance would be appreciated.
Riemann invariants are not unique. They are functions of $(u,\eta)$ that are constant on characteristics, but if $R$ is constant on a curve then so is $2+R^2$, or any function of $R$.
Where you wrote "e.g." is a mistake because your example is not a gradient, by the mixed partial derivative test. It needs a multiple of that, an integrating factor.