Consider the vector field $$ \mathbf{f}(x_1,x_2)=\begin{bmatrix}f_1(x_1,x_2)\\ f_2(x_1,x_2)\end{bmatrix}=\begin{bmatrix}\frac{1}{1+x_2}+x_1\\ \frac{1}{1+x_1}+x_2\end{bmatrix}. $$
I'm interested in finding a scalar function $\psi(x_1,x_2)$ such that: $$\tag{$\ast$}\label{ast} \langle\nabla\psi,\nabla\psi+\mathbf{f} \rangle= \frac{\partial\psi(x_1,x_2)}{\partial x_1} \left(f_1(x_1,x_2)+\frac{\partial\psi(x_1,x_2)}{\partial x_1}\right)+\frac{\partial\psi(x_1,x_2)}{\partial x_2} \left(f_2(x_1,x_2)+\frac{\partial\psi(x_1,x_2)}{\partial x_2}\right)=0 $$
My question. Does there exist a solution $\psi(x_1,x_2)$ to \eqref{ast}? If so, how to compute it?
I'm aware that my question could be a trivial or naive one, but since I'm new to this kind of problems I would really appreciate any help/comment/suggestion. Thank you.
On the one hand, the nullity of the scalar product $$ \nabla\psi\cdot (f+\nabla\psi) = \psi_{,1} (f_1+ \psi_{,1} ) + \psi_{,2} (f_2+ \psi_{,2} ) = 0 $$ leads to the orthogonality of $\nabla\psi$ and $f+\nabla\psi$. Thus, $f+\nabla\psi$ is proportional to the vector $\nabla\psi^\perp = (\psi_{,1}, -\psi_{,2})^\top$. There exists $\alpha(x_1,x_2)$ such that $f+\nabla\psi = \alpha \nabla\psi^\perp$, or equivalently $$ (1-\alpha)\psi_{,1} = -f_1 \qquad\text{and}\qquad (1+\alpha)\psi_{,2} = -f_2 . $$ By eliminating successively $\alpha$ and unity of the previous system, we get to \begin{aligned} 2\psi_{,1}\psi_{,2} &= -f_2\psi_{,1} - f_1\psi_{,2} \\ 2\alpha\psi_{,1}\psi_{,2} &= -f_2\psi_{,1} + f_1\psi_{,2} . \end{aligned}
On the other hand, the nullity of the scalar product rewrites as $$ (\psi_{,1} + \psi_{,2})^2 - 2\psi_{,1}\psi_{,2} + f_1\psi_{,1} + f_2\psi_{,2} = 0 , $$ so that previous identities lead to $$ (\psi_{,1} + \psi_{,2})^2 + (f_1+f_2)(\psi_{,1} + \psi_{,2}) = 0 . $$ Thus, the solutions satisfy $\psi_{,1} + \psi_{,2} = -(f_1+f_2)$ if they have nonzero divergence -- otherwise, they satisfy $\psi_{,1} + \psi_{,2} = 0$, which is a particular case of the previous non-homogeneous equation. Solutions to these linear first-order PDEs can be obtained by using the method of characteristics for the Lagrange-Charpit system $\text d x_1 = \text d x_2 = {\text d \psi}/{r}$ with $r = -(f_1+f_2)$ (see related posts on this site). Solutions are of the form $$ \psi(x_1, x_2) = \int^{x_1} r(\xi,\xi + x_2-x_1)\, \text d \xi + F(x_2-x_1) , $$ where $F$ is an arbitrary function. The latter is determined by injecting the previous expression in the orthogonality condition $\nabla\psi\cdot (f+\nabla\psi) = 0$.
If $f$ is constant, then linear solutions of the form \begin{aligned} \psi(x_1,x_2) &= -f_1x_1 - f_2x_2 + C \\ \text{or}\qquad \psi(x_1,x_2) &= -\tfrac12(f_1+f_2)(x_1+x_2)+ C \end{aligned} are obtained. In particular, one notes that the solutions are constant if $f \equiv 0$. In theory, such solutions can be obtained similarly in the case $$ f_1(x_1,x_2) = x_1 + (1+x_2)^{-1}, \qquad f_2(x_1,x_2) = x_2 + (1+x_1)^{-1} , $$ but computations are more involved.