$Ric-\lambda g$ equaling to Hessian of a function implies Einstein I am studying some basic knowledge of Riemannian Geometry and I am stucked in the following problem:

$(M,g)$ is a compact orientable Riemannian manifold and there exists $f \in C^{\infty}(M)$ and $\lambda \in \mathbb R_{\leq 0}$ such that 
  $\text{Ric} - \lambda g =\text{Hess}(f)$ as $(0,2)$-tensor over $M$. Prove that $f$ is constant hence $M$ is Einstein.

My thought in the problem is that the Second Bianchi identity may help as the left hand side involves $Ric$ and $g$, but the right hand side may be tedious after taking covariant derivatives. Also, I think maybe some identity like Bochner  formula may help as the condition $\lambda \leq 0$ may force something be zero after integral, but the left hand side will be not easy to handle.
In other words, I cannot find a good way to characterize Hessians among symmetric bilinear tensors which is discussed in this problem.
 A: We have that $$div(Ric-\lambda g)=div(Ric)=\dfrac12 dR,$$ where $R$ stands for the scalar curvature, and $$div(Hess(f))=Ric(\nabla f,\cdot)+g(\nabla \Delta f, \cdot).$$ Using that $$Ric-\lambda g=Hess(f)$$ in the last equality we get
$$div(Hess(f))=\lambda g(\nabla f,\cdot)+Hess(f)(\nabla f,\cdot)+g(\nabla \Delta f, \cdot).$$ Now $$Hess(f)(\nabla f,\cdot)=g(\nabla_{\nabla f}\nabla f,\cdot)=\dfrac12 d|\nabla f|^2.$$ So 
$$div(Hess(f))=\lambda g(\nabla f,\cdot)+\dfrac12 d|\nabla f|^2+g(\nabla \Delta f, \cdot).$$ Thus
$$\dfrac12 dR=\lambda df+\dfrac12 d|\nabla f|^2+d\Delta f.$$ That is, there exists a constant $c$ such that
$$R=2\lambda f+|\nabla f|^2+2\Delta f+c.$$ Since $R=n\lambda +\Delta f$ we get
$$2\lambda f+|\nabla f|^2+\Delta f=C,$$ where $C=n\lambda -c.$ Since $M$ is compact there exists $p,q\in M$ such that $f(p)\le f\le f(q).$ So
$$2\lambda f(p)+(\Delta f)(p)=2\lambda f(q)+(\Delta f)(q)$$ from where
$$(\Delta f)(p)-(\Delta f)(q)=2\lambda (f(q)-f(p)).$$ If $\lambda <0$ then we get 
$$f(q)-f(p)\le 0$$ since $$(\Delta f)(p)-(\Delta f)(q)\ge 0.$$ That is, $f$ is constant.
If $\lambda=0$ we have that $$|\nabla f|^2+\Delta f=C.$$ Thus
$$0\le (\Delta f)(p)=C, 0\ge (\Delta f)(q)=C,$$ from where $C=0.$ Integrating we get
$$\int_M|\nabla f|^2=-\int_M\Delta f=0$$ from where $f$ is constant.
