Since the proof that such a function implies the metric is a warped product metric is fairly simple, I include a complete copy below.

### 1

Start with $\nabla^2 f = \lambda g$. For any vector field it holds

$$ \frac12 \nabla_a \nabla_b(X^c X_c) = (\nabla_a \nabla_b X^c) X_c - \nabla_b X^c \nabla_a X_c $$

which implies, after taking the antisymmetric part and using that Hessians of scalars are symmetric, that

$$ (\nabla_a \nabla_b X_c) X^c - (\nabla_b\nabla_a X_c) X^c = 0 $$

Now specializing to $X = \nabla f$ the gradient, we find

$$ \nabla_a \lambda X_b - \nabla_b \lambda X_a = 0 $$

which shows that $d\lambda$ is collinear with $df$; or that $\lambda$ is constant on the level sets of $f$.

### 2

$\nabla_X X = \lambda X$ and so $X$ is (pre-)geodesic.

$\nabla g(X,X) = 2 g(X, \nabla X) = 2 \lambda X$, hence $g(X,X)$ is constant along level sets of $f$.

Letting $Y,Z$ be orthogonal to $X$ (hence tangent to level sets of $f$), we find

$$ - g(II(Y,Z),X) = - g(X, \nabla_Y Z) = g(\nabla_Y X,Z) = \lambda g(Y,Z) $$

As $X$ is a normal to the level set of $f$ with $g(X,X)$ fixed along said level set, this shows that the level sets of $f$ is totally umbilic with constant proportionality factor.

### 3

The Lie derivative $\mathcal{L}_X g = 2 \lambda g$ by assumption (so $X$ is a homothety). Since $X$ is hypersurface orthogonal, we can write

$$ g = \frac{X^\flat X^\flat}{g(X,X)} + h $$

where $h$ is the induced metric on the level sets of $f$. Taking the Lie derivative we find, after a computation, that this implies

$$ \mathcal{L}_X h = 2 \lambda h $$

As $\lambda$ is constant along the leaves of $f$, we can find a function $F$ that is constant on the leaves of $f$ so that $X(F) = - 2 \lambda$. Then writing $\tilde{h} = e^F h$, we find that $\mathcal{L}_X \tilde{h} = 0$. So we can write

$$ g = \frac{X^\flat X^\flat}{g(X,X)} + e^{-2F} \tilde{h} $$

where $\mathcal{L}_X \tilde{h} = 0$. This is precisely the warped product form.