An identity between the d'Alembertian and the covariant derivative Suppose $f$ is function which depends on $\phi$, $f = f(\phi)$; and $\phi$ is a scalar field. We define
$$\square \equiv g^{\mu\nu} \nabla_\mu \nabla_\nu$$ and $$\nabla_\mu f(\phi) \equiv  f_{;\nu}$$
Is this expression below correct?
$$ -f_{;\mu\nu} + \square f g_{\mu\nu} = 0  $$
I think it is correct since
$$  \nabla_\nu(\nabla_\mu f) =  f_{;\mu\nu} $$ and
$$(g^{\mu\nu} \nabla_\mu \nabla_\nu f) g_{\mu\nu} = \square f g_{\mu\nu} = (\nabla_\mu \nabla_\nu f)g^{\mu\nu}g_{\mu\nu} =\nabla_\mu \nabla_\nu f$$
but I doubt it, I need another point of view.
ps: for further, you can check arXiv:2009.11827. I try to prove equation (18).
 A: This expression:
$$(g^{\mu\nu} \nabla_\mu \nabla_\nu f) g_{\mu\nu} = \square f g_{\mu\nu} = (\nabla_\mu \nabla_\nu f)g^{\mu\nu}g_{\mu\nu} =\nabla_\mu \nabla_\nu f,$$
is wrong. This does not hold. The Box operator is a scalar quantity. It is defined as:
$$\Box = g^{\mu\nu}\nabla_{\mu}\nabla_{\nu}$$
and for a four dimensional diagonal metric this is:
$$\Box = g^{tt}\nabla_{t}\nabla_{t} + g^{xx}\nabla_{x}\nabla_{x} + g^{yy}\nabla_{y}\nabla_{y} + g^{zz}\nabla_{z}\nabla_{z}$$
You abused notation and thus you're not able to see this but the indices for the box term $\mu, \nu$ denote summation (dummy indices) while the indices for the metric denote components (free indices). A more correct way do write it would be:
$$\Box f g_{\mu\nu} =  (g^{\alpha\beta}\nabla_{\alpha}\nabla_{\beta} f) g_{\mu\nu} $$
Also $g_{\mu\nu}g^{\mu\nu} =d$ where $d$ is the dimensionality of spacetime for a diagonal metric. I hope this helps.
A: You are trying to prove eq (18) in 2009.11827. That has nothing to do with the identity you are asking about in the question (which is completely incorrect, btw). To derive (18) in the ref. we simply need to vary the action w.r.t. the metric ($G_4$ is a function only of $\varphi$)
$$
\delta \int \sqrt{-g} G_4 R = \int \left[ \delta  \sqrt{-g} G_4R + \sqrt{-g} G_4\delta  R \right]
$$
Next, we use
$$
\delta \sqrt{-g} = - \frac{1}{2} \sqrt{-g} g_{\mu\nu} \delta g^{\mu\nu} , \qquad \delta R = R_{\mu\nu} \delta g^{\mu\nu} - \delta g^{\mu\nu}{}_{;\mu\nu} + g_{\mu\nu} \Box \delta g^{\mu\nu} 
$$
Then,
\begin{align}
\delta \int \sqrt{-g} G_4R &= \int \sqrt{-g}  \left[ - \frac{1}{2} g_{\mu\nu} \delta g^{\mu\nu}  G_4R +  G_4 R_{\mu\nu} \delta g^{\mu\nu} - G_4 \delta g^{\mu\nu}{}_{;\mu\nu} + G_4 g_{\mu\nu} \Box \delta g^{\mu\nu}  \right]\\
&= \int \sqrt{-g}  \left[ G_4 E_{\mu\nu} - G_{4;\mu\nu}  + g_{\mu\nu}  \Box G_4 \right] \delta g^{\mu\nu} 
\end{align}
Thus, the equations of motion are
$$
{\cal E}_{\mu\nu}^{(4)} = G_4 E_{\mu\nu} - G_{4;\mu\nu}  + g_{\mu\nu}  \Box G_4
$$
The identity in the question is not correct and is also not relevant to the calculation.
A: Your expression is correct if you're using a derivative that is compatible with the metric, so if $\nabla_\mu g_{\rho\sigma}=0$, otherwise the d'Alembertian would bring other terms too.
Edit: The reason why I think this is correct is
$\Box f g_{\rho\sigma}=g^{\mu\nu}\nabla_\mu\nabla_\nu f g_{\rho\sigma}=\nabla^\nu\nabla_\nu f g_{\rho\sigma}$
$f_{,\rho\sigma}=\nabla_\rho\nabla_\sigma f =\nabla^\rho\nabla_\rho f g_{\rho\sigma}=\nabla^\nu\nabla_\nu f g_{\rho\nu}g^{\rho\nu}g_{\rho\sigma}=\nabla^\nu\nabla_\nu f g_{\rho\sigma}$
I may be wrong, as @Prahar says, but I can't see where.
