What does it mean that the Hessian is proportional to the metric? Let $(M,g)$ be a smooth manifold equipped with a metric tensor $g$, and $f\in C^\infty(M)$ a regular function (i.e., with nowhere vanishing differential).
Denote by $\mathrm{Hess}_g(f):=\nabla df$ the Hessian tensor of $f$ with respect to the metric $g$, and by $N_f:=f^{-1}(\{0\})$ the one-codimensional submanifold of $M$ determined by $f$.

QUESTION: how to characterise the functions $f$, with $\mathrm{Hess}_g(f)$ proportional to $g$? Do they belong to some well-known class? What about the submanifolds $N\subset M$, which can be written as $N=N_f$, for such an $f$?

By "proportional" I mean w.r.t. a conformal factor, and by "characterisation" I mean "what makes them special". The question may be restricted to the case when $g$ is flat.
Motivating example. Let $M=\mathbb{R}^2$ and $g$ the Euclidean metric. Then,
$$\mathrm{Hess}_g(f)=f_{xx}(dx)^2+f_{xy}dxdy+f_{yy}(dy)^2$$ 
is proportional to $g=(dx)^2+(dy)^2$ if and only if $f_{xy}=0$ and $f_{xx}=f_{yy}$, which means
$$
f=a(x^2+y^2)+bx+cy+d\, ,
$$
i.e., $f$ must be proportional to the square of the norm induced by $g$, plus some lower-order terms. I was wondering whether a similar result holds in general (of course, assuming $g$ to be flat).
 A: It is known (say, Y. Tashiro, Complete Riemannian manifolds and some vector fields, Trans.Amer.Math.Soc. 117(1965) 251–
275; I am not sure  that Tashiro is the first who proved it and there were many later papers which independently prove the same result later. The proof is pretty straightforward.)  that the existence of such a function for  a metric implies that the metric is a warped product metric, 
$$
g= dx_1^2 + \sigma(x_1) \sum_{i,j=2}^{n} h_{ij}(x_2,...,x_n)dx_i dx_j,    \ \ \ (*)
$$ 
and for such metrics a solution $f$ of your  equation    $\nabla\nabla f= \lambda(x) g$ is some  function of $\sigma$. 
The special case when the metric $g$ is flat is much more easy. In this case  the function $\sigma(x_1)$  from $(*)$
is automatically $\textrm{const_1} \ (x_1+ \textrm{const}_2)^2$ (this fact is also classically known), so in the nontrivial case, when $\textrm{const_1}\ne 0$, after the affine reparameterization  of $x_1$,   the form  $(*)$    is the cone form  $$ 
g= dx_1^2 + x_1^2 \sum_{i,j=2}^{n} h_{ij}(x_2,...,x_n)dx_i dx_j,   
$$ 
This implies that in the flat case  the function $\sigma(x_1)$, in  some euclidean  coordinate system $(y_1,...,y_n)$, is  (after addition of linear terms and constant which change nothing) simply the function 
$$f= y_1^2 +...+ y^2_n.$$ 
For flat metrics of other signatures the answer is essentially the same,  in this case $f= \varepsilon_1y_1^2 +...+ \varepsilon_ny^2_n,$ where
 $\varepsilon_i\in \{\pm 1\}$ are responcible for the signature
A: It is rare that manifold admits such functions.
Assume $Hess_g(f)=\lambda\cdot g$.
If $\lambda=0$ then it is so called affine function on $M$.
In this case $M$ is isometric to $\mathbb R\times M'$ and the function $f$ depends linearly on the first projection.
If $\lambda\ne 0$ the picture is similar, but you get so called warped product. 
The level sets of your function formed by manifolds with constant normal curvature.
