Derivation of Klein-Gordon equation in General Relativity I am trying to derive the Klein-Gordon equation for the case of GR using the action:
$$S\left[ {\varphi ,{g_{\mu \nu }}} \right] = \int {\sqrt g {d^4}x\left( { - {1 \over 2}{g^{\mu \nu }}{\nabla _\mu }\varphi {\nabla _\nu }\varphi  - {1 \over 2}{m^2}{\varphi ^2}} \right)}   \tag{1}$$
So in the Euler - Lagrange equation for the SR case which is:
$${\partial _\mu }\left( {{{\partial L} \over {\partial \left( {{\partial _\mu }\varphi } \right)}}} \right) = {{\partial L} \over {\partial \varphi }}\tag{2}$$
I use the correspondence 
$${\partial  \over {\partial \left( {{\partial _\mu }\varphi } \right)}} \to {\nabla _{{\nabla _\mu }\varphi }},{\partial  \over {{\partial _\mu }\varphi }} \to {\nabla _\mu } \tag{3}$$
and so,
$${\nabla _{{\nabla _\mu }\varphi }}L =  - {1 \over 2}\sqrt g {g^{\mu \nu }}{\nabla _\nu }\varphi  - {1 \over 2}\sqrt g {g^{\mu \nu }}{\nabla _\mu }\varphi {\nabla _{{\nabla _\mu }\varphi }}\left( {{\nabla _\nu }\varphi } \right) \tag{4}$$
and using this dubious relation that I am unable to prove (extending the SR case of partial derivatives),
$${\nabla _{{\nabla _\mu }\varphi }}\left( {{\nabla _\nu }\varphi } \right) = \delta _\nu ^\mu  \tag{5}$$
I finally get the desired KG equation,
$${\nabla _\mu }\left( { - \sqrt g {g^{\mu \nu }}{\nabla _\nu }\varphi } \right) =  - \sqrt g {m^2}\varphi $$
However, I am very uncomfortable with my assumption 5 that I am unable to prove. Is my analysis correct? I am a 60y old doing this as retirement fun so please don't cut me down too brutally :-)
 A: It's probably more straightforward to express your action in terms of partial derivatives:
$$ S\left[ {\varphi ,{g_{\mu \nu }}} \right] = \int {\sqrt g {d^4}x\left( { - {1 \over 2}{g^{\mu \nu }}{\partial_\mu }\varphi \,{\partial_\nu }\varphi  - {1 \over 2}{m^2}{\varphi ^2}} \right)}  \,, $$
since $\varphi$ is just a scalar. Then
$$ \frac{\partial L}{\partial (\partial_\mu\varphi)} = -\sqrt{g}\,g^{\mu \nu} \partial_\nu \varphi \quad\text{and}\quad \frac{\partial L}{\partial \varphi} = -\sqrt{g}\,m^2 \varphi^ \,.$$
Finally we need the useful expression for the covariant divergence:
$$ \nabla_\mu V^\mu = \frac{1}{\sqrt{g}} \partial_\mu(\sqrt{g}V^\mu) \,,$$
to reproduce the usual Klein-Gordon equation.

The correspondence you state in $(3)$ is not correct. The normal procedure for turning a theory in flat spacetime into one in curved spacetime involves replacing $\partial_\mu$ by $\nabla_\mu$, yes, but the Euler-Lagrange equations don't change. The E-L equations can be derived for any action, whether in flat spacetime, curved spacetime, or in some other context entirely. One clue that we shouldn't be replacing partial derivatives by covariant derivatives here is the fact that the derivatives are with respect to fields; in the minimal coupling prescription, $\partial_\mu \to \nabla_\mu$ is a statement about how we should take derivatives with respect to position in spacetime.
A: Note: technically, everywhere we've been using $\sqrt{g}$ it should be $\sqrt{|g|}$, but I'll let the former denote the latter in a slight abuse of notation like the rest of this page.
gj255 has already provided the textbook approach, in which we keep to partial derivatives. You can actually do it a different way; the same from-first-principles argument that derives Eq. (2) can be modified to obtain something that manifestly works with covariant derivatives. Write the Lagrangian density as $L=\sqrt{g}L_0$ with $L_0=\frac{1}{2}\nabla_\mu\varphi\nabla^\mu\varphi-\frac{1}{2}m^2\varphi^2$, so $$0=\delta S=\int d^4x \sqrt{g}\left(\delta\varphi \frac{\partial L_0}{\partial\varphi}+\nabla_\mu\delta\varphi \frac{\partial L_0}{\partial\nabla_\mu\varphi}\right)$$($\nabla_\mu\varphi$ is actually just $\partial_\mu\varphi$, because $\varphi$ is a scalar field), where we have written $\delta\nabla_\mu\varphi=\nabla_\mu\delta\varphi$. But$$\int d^4x\sqrt{g}\nabla_\mu\left(\delta\varphi \frac{\partial L_0}{\partial\nabla_\mu\varphi}\right)=\int d^4x\partial_\mu\left(\sqrt{g}\delta\varphi \frac{\partial L_0}{\partial\nabla_\mu\varphi}\right)$$is a boundary term, so we can rewrite $\delta S=0$ as $$0=\int d^4 x\sqrt{g}\delta\varphi \left(\frac{\partial L_0}{\partial\varphi}-\nabla_\mu\frac{\partial L_0}{\partial\nabla_\mu\varphi}\right),$$giving us an Euler-Lagrange function more in the spirit of what you wanted, viz.$$0=\frac{\partial L_0}{\partial\varphi}-\nabla_\mu\frac{\partial L_0}{\partial\nabla_\mu\varphi}=-m^2\varphi-\nabla_\mu\nabla^\mu\varphi.$$Note: I don't know what the standard symbol is for $L/\sqrt{g}$, but I do know it's sometimes called the scalar Lagrangian density.
A: In order to derive the Klein-Gordon equation you must vary with respect to the scalar field $\phi$. The action reads: 
$$S = \int d^{4}x\sqrt{-g}\left(-\cfrac{1}{2}g^{μν}\nabla_{μ}\phi\nabla_{ν}\phi - \cfrac{1}{2}m^{2}\phi ^{2}\right) $$
For the kinetic term you have: 
\begin{align}
δ(g^{μν}\nabla_{μ}\phi\nabla_{ν}\phi) &= g^{μν}δ(\nabla_{μ}\phi)\nabla_{ν}\phi + g^{μν}δ(\nabla_{ν}\phi)\nabla_{μ}\phi\\
& = δ(\nabla_{μ}\phi)\nabla^{μ}\phi + δ(\nabla_{ν}\phi)\nabla^{ν}\phi \\ 
& = 2δ(\nabla_{μ}\phi)\nabla^{μ}\phi \\
& = 2\nabla_{μ}(δ\phi)\nabla^{μ}\phi
\end{align}
where I've summed the terms with dummy indices. From Leibnitz rule we know that:  
$$ \nabla_{μ}(δ\phi\nabla^{μ}\phi) = \nabla_{μ}(δ\phi)\nabla^{μ}\phi + δ\phi\nabla_{μ}\nabla^{μ}\phi $$
The left hand side term in the above expression in spacetime takes the following form:
$$\int d^{4}x\sqrt{-g} \nabla_{μ}(δ\phi\nabla^{μ}\phi).$$
We can see that this term is a surface term so: 
$$\int d^{4}x\sqrt{-g} \nabla_{μ}(δ\phi\nabla^{μ}\phi) = 0$$
because $δ\phi = 0$ at the surface. So we have:  
$$\int d^{4}x \sqrt{-g}\nabla_{μ}(δ\phi)\nabla^{μ}\phi =  -\int d^{4}x \sqrt{-g}δ\phi\nabla_{μ}\nabla^{μ}\phi $$
The mass term is: 
$$δ(m^{2} \phi ^{2}) = 2\phi δ\phi m^{2}$$
Combining all the above relations one obtains the K-G equation for a massive scalar field: 
$$(\Box  - m^{2})\phi = 0$$
where $\Box = g^{μν}\nabla_{μ}\nabla_{ν}$ is the D'Alambert operator.
Note - Hint: You can derive K-G equation using the conservation of the energy momentum tensor: $\nabla^{μ}Τ_{μν} = 0$
