Let us first recall some definitions and useful formulas for a surface in $\mathbb{R}^3$ given by an immersion $F\colon U \rightarrow \mathbb{R}^3$, where $U \subset \mathbb{R}^2$ is an open set.
Denote $F_x = \frac{\partial F}{\partial{x}}$, $F_{x x} = \frac{\partial{F_x}}{\partial{x}}$, and so on.
The components of the first fundamental form are $I_{x x} = F_x \cdot F_x$, $I_{y y}=F_y \cdot F_y$, $I_{x y} = I_{y x} = F_x \cdot F_y$
For the second fundamental form we may use the expression for the unit normal vector
$$
n = \frac{F_x \times F_y}{|F_x \times F_y|}
$$
so that $II_{x x} = n \cdot F_{x x}$, $II_{y y} = n \cdot F_{y y}$, $II_{x y} = II_{y x} = n \cdot F_{x y}$
The Gaussian curvature $K(x,y)$ has the following expression
$$
K(x,y) = \frac{II_{x x} II_{y y} - II_{x y}^2}{I_{x x} I_{y y} - I_{x y}^2} \tag{1}
$$
and the mean curvature can be computed by
$$
H(x,y) = \frac{I_{x x} II_{y y} - 2 I_{x y} II_{x y} + I_{y y} II_{x x}}{I_{x x} I_{y y } - I_{x y}^2} \tag{2}
$$
In the proposed problem we a given a surface represented by a graph of function $z = f(x) + g(y)$, so our immersion has the following form:
$$
F(x,y) = \begin{pmatrix}
x \\
y \\
f(x) + g(y)
\end{pmatrix}
$$
and we calculate $F_x = \begin{pmatrix} 1 \\ 0 \\ f' \end{pmatrix}$, $F_y = \begin{pmatrix} 0 \\ 1 \\ g' \end{pmatrix}$, $F_{x x} = \begin{pmatrix} 0 \\ 0 \\ f'' \end{pmatrix}$, $F_{y y} = \begin{pmatrix} 0 \\ 0 \\ g'' \end{pmatrix}$, $F_{x y} = \begin{pmatrix} 0 \\ 0 \\ 0 \end{pmatrix}$.
This is enough to find the unit normal
$$
n = \frac{(f', -g', 1)^T}{\sqrt{1 + (f')^2 + (g')^2}}
$$
so we find the components of the second fundamental form $II_{x x} = \frac{f''}{\sqrt{1 + (f')^2 + (g')^2}}$, $II_{y y} = \frac{g''}{\sqrt{1 + (f')^2 + (g')^2}}$, $II_{x y} = 0$
The components of the first fundamental form are, of course, $I_{x x} = 1 + (f')^2$, $I_{y y} = 1 + (g')^2$, and $I_{x y} = f'g'$.
I would let you to finish the job by substituting these quantities into equations (1) and (2).
Yes, you can compute all the coefficients $e,f,g,E,F,G$ and get the gaussian and mean curvature and yes, it's tedious.
Here's another way:
From the first step we get : $Y_u\times Y_v=(1-2Ha+Ka^2)(X_u\times X_v)$, ie if $N$ and $\overline N$ are the normal vectors of $X$ and $Y$ respectively, then $\overline N\circ Y$ and $N\circ X$ coincide, since they're parallel. If these functions coincide then we have the following relations :
$$d\overline N(Y_u)=(\overline N\circ Y)_u=(N\circ X)_u=dN(X_u) \tag1$$
$$d\overline N(Y_v)=(\overline N\circ Y)_v=(N\circ X)_v=dN(X_v) \tag2$$
Let $\overline B$ be the matrix of $d\overline N$ with respect to $\{Y_u,Y_v\}$ and $B$ the matrix of $dN$ with respect to $\{X_u,X_v\}$.
Now, to compute $\overline K$ and $\overline H$ we need to find the expression of $\overline B$.
Put $$B=\begin{bmatrix}b_{11} & b_{12}\\ b_{21} & b_{22}\\ \end{bmatrix}$$
From the definition of $Y$ we have:
$$Y_u=X_u+a\cdot N_u=(a\cdot b_{11}+1)\cdot X_u+a\cdot b_{21}\cdot X_v$$
$$Y_v=X_v+a\cdot N_v=a\cdot b_{12}\cdot X_u+(a\cdot b_{22}+1)\cdot X_v$$
From these equations we can get the "change of basis" matrix : $Q=\begin{bmatrix}a\cdot b_{11}+1 & a\cdot b_{12}\\ a\cdot b_{21} & a\cdot b_{22}+1\\ \end{bmatrix}$ from $\{X_u,X_v\}$ to $\{Y_u,Y_v\}$. Then from the initial relations $(1)$ and $(2)$, we have the following equation:
$$B=Q\cdot \overline B$$
Since $Q$ is invertible: $$ \overline B=Q^{-1}\cdot B$$ From this point you can compute the entries of $\overline B$ and calculate $\overline H $ and $ \overline K$.
You can also notice that, since $Q^{-1}=(I+a\cdot B)^{-1}$, you have $\overline B=(I+a\cdot B)^{-1}\cdot B $. So, if $B$ has eigenvalues $-\lambda_1$ and $-\lambda_2$, then the eigenvalues of $\overline B$ are $\frac{-\lambda_1}{1-a\cdot \lambda_1}$ and $\frac{-\lambda_2}{1-a\cdot \lambda_2}$ and you can easily compute $\overline H$ and $\overline K$.
Best Answer
We've that $dN_p (X_u) = \lambda X_u$, because the points are all umbilical. But you have just showed that $\langle dN_p (X_u) , X_u \rangle = 0 $. Thus, we have that $dN_p (X_u) = 0$. The same arguing shows that $dN_p (X_v) = 0$. This plainly implies that $N$ is constant, and we've that $S$ is contained in a plane, as desired.