S. Carroll in his book tells that the geodesic equation (together with metric compatibility) implies that the quantity $$\epsilon =-g_{\mu \nu}\frac{dx^{\mu}}{d\lambda}\frac{dx^{\nu}}{d\lambda}\tag{5.55}$$ is constant along the geodesic. I couldn't prove it. I would really appreciate if someone could help me about it.
My try: We must show that $\bigtriangledown_{\rho}(\epsilon)=0$. By metric compatibility, we have $\bigtriangledown_{\rho}(\epsilon)=-g_{\mu \nu}\bigtriangledown_{\rho} (\frac{dx^{\mu}}{d\lambda}\frac{dx^{\nu}}{d\lambda})$. Now by definition of $\bigtriangledown_{\rho}$, we get $\bigtriangledown_{\rho}(\epsilon)=-g_{\mu \nu}\big[ \partial_{\rho}(\frac{dx^{\mu}}{d\lambda}\frac{dx^{\nu}}{d\lambda})+\Gamma^{\mu}_{\rho \sigma}\frac{dx^{\sigma}}{d\lambda}\frac{dx^{\nu}}{d\lambda}+\Gamma^{\nu}_{\rho \sigma}\frac{dx^{\mu}}{d\lambda}\frac{dx^{\sigma}}{d\lambda}\Big]$. On the one hand, we have $$\partial_{\rho}(\frac{dx^{\mu}}{d\lambda}\frac{dx^{\nu}}{d\lambda})=\partial_{\rho}(\frac{dx^{\mu}}{d\lambda})\frac{dx^{\nu}}{d\lambda}+\frac{dx^{\mu}}{d\lambda}\partial_{\rho}(\frac{dx^{\nu}}{d\lambda})=\frac{\partial}{\partial x^{\rho}}(\frac{dx^{\mu}}{d\lambda})\frac{dx^{\nu}}{d\lambda}+\frac{dx^{\mu}}{d\lambda}\frac{\partial}{\partial x^{\rho}}(\frac{dx^{\nu}}{d\lambda})$$ $$=\delta^{\mu}_{\rho}\frac{d}{d\lambda}\frac{dx^{\nu}}{d\lambda}+\frac{dx^{\mu}}{d\lambda}\delta^{\nu}_{\rho}\frac{d}{d\lambda}=\delta^{\mu}_{\rho}\frac{d^2 x^{\nu}}{d\lambda^2}+\delta^{\nu}_{\rho}\frac{d^2 x^{\mu}}{d\lambda^2}.$$
On the other hand, from the geodesic equation we have $\frac{d^2 x^{\mu}}{d\lambda^2}=-\Gamma^{\mu}_{\rho \sigma}\frac{dx^{\rho}}{d\lambda}\frac{dx^{\sigma}}{d\lambda}$. But now I don't know how get the result.