This a continuation of a previous question, which you can find here.
Question: Let $M$ be a Riemann surface with constant curvature, but assume $M$ is not flat. In fact, I would be satisfied with an analysis of the case $M = S^2$. Can there exist two pointwise linearly-independent vector fields $X$ and $Y$ defined on a small open subset of $M$ such that each is parallel along the other? Here, parallel means parallel with respect to the Levi-Cevita connection.
I haven't thought of an organized way to approach this problem, but I have ground out some calculations, which may be of use. Let $g$ be the Riemann metric and let $\nabla$ be its Levi-Cevita connection. Since $\nabla$ is torsion free, and the vector fields are each parallel along the other, we have $[X,Y] = \nabla_X Y - \nabla_Y X = 0$. Thus, the vector fields give us a with coordinates $(x,y)$ in $M$. In this coordinate system, the vector fields are just the standard ones $X = e_1$, $Y = e_2$. The Riemann metric $g = \begin{pmatrix} g_{11} & g_{12} \\ g_{21} & g_{22} \end{pmatrix} = \begin{pmatrix} u & \beta \\ \beta & v \end{pmatrix}$ is a smoothly varying positive definite matrix i.e. each of $u$, $v$ and $\Delta = uv - \beta^2$ is a positive function. One calculates that $$ \nabla_{e_1} e_2 = \frac{1}{2} g^{-1} \begin{pmatrix} u_y \\ v_x \end{pmatrix}$$ and so the condition that $e_2$ be parallel along $e_1$ (since these vector fields commute and the torsion vanishes, it follows that $e_1$ is parallel along $e_2$) is precisely that $u_y = v_x = 0$. In other words $u$ is only a function of $x$ and $v$ is only a function of $y$. Next, one calculates that \begin{align*} \nabla_{e_1} e_1 &= \frac{1}{2} g^{-1} \begin{pmatrix} u_x \\ 2 \beta_x \end{pmatrix} = \frac{1}{2\Delta} \begin{pmatrix} vu_x-2\beta\beta_x \\ -\beta u_x+2u\beta_x \end{pmatrix} = \frac{1}{2\Delta} \begin{pmatrix} \Delta_x \\ -\beta u_x+2u\beta_x \end{pmatrix} \\ \nabla_{e_2} e_2 &= \frac{1}{2} g^{-1} \begin{pmatrix} 2 \beta_y \\ v_y \end{pmatrix} =\frac{1}{2\Delta} \begin{pmatrix} 2v\beta_y - \beta v_y \\ -2 \beta \beta_y + u v_y \end{pmatrix} =\frac{1}{2\Delta} \begin{pmatrix} 2v\beta_y - \beta v_y \\ \Delta_y \end{pmatrix} \end{align*} The Christoffel symbols of the connection can be read off from here. These, in turn, can be used to express the curvature, which we could then equate to a constant, ending up with a complicated nonlinear PDE to ponder. I wasn't able to rule out solutions.
If one makes the extra assumption:
Simplying assumption: Suppose that $\Delta = \det(g)$ is identically equal to $1$.
the calculations become rather more tractable. In fact, earlier I thought I had proved the curvature needs to vanish under this extra hypothesis... but looking again at my calculations, I am unable to follows them myself, so I have not put them here.