This is the differential form of the Reilly formula. It holds for a function on any pseudo-Riemannian manifold. (Robert C. Reilly. _Applications of the Hessian operator in a Riemannian manifold_, Indiana Univ. Math. J. **26** (1977), no. 3, 459–472, doi:[10.1512/iumj.1977.26.26036](https://doi.org/10.1512/iumj.1977.26.26036))

Use the product rule to say
$$(\Delta f)^2=\operatorname{div}(\Delta f\cdot\nabla f)-\langle\nabla f,\nabla\Delta f\rangle.$$
Use the commutation formula for covariant derivatives to replace the last term by
$$\langle\nabla f,\nabla\Delta f\rangle=\langle\nabla f,\Delta\nabla f\rangle-\operatorname{Ric}(\nabla f,\nabla f).$$
Use the product rule to replace the second to last term by
$$\langle\nabla f,\Delta\nabla f\rangle=\operatorname{div}\big(\nabla^2f(\nabla f,\cdot)\big)-|\nabla\nabla f|^2.$$
Finally $\nabla^2f(\nabla f,\cdot)=\frac{1}{2}\nabla|\nabla f|^2$. This gives your formula.

**Edit**. As pointed out by Jeffrey Case below, this also follows from the Bochner formula
$$\frac{1}{2}\Delta|\nabla f|^2=|\nabla\nabla f|^2+\langle\nabla\Delta f,\nabla f\rangle+\operatorname{Ric}(\nabla f,\nabla f),$$
where you just need to use the very first line above to replace the middle term on the RHS. The proof of the Bochner formula is by the other lines above