Euler-Lagrange equation in curved spacetime The action of a field $\phi^\mu$ in flat $n$-dimensional spacetime is
$$ S = \int \text{d}^n x \mathscr{L}(\phi^\mu(x),\partial_\alpha \phi^\mu(x)) $$
From an infinitesimal variation of field configuration, and discarding a boundary term, we have E-L equations
$$ \frac{\partial \mathscr{L}}{\partial \phi^\mu} - \partial_\alpha \frac{\partial \mathscr{L}}{\partial (\partial_\alpha \phi^\mu)} $$
In curved spacetime the action should now be
$$ S = \int \text{d}^n x \sqrt{|g|} \mathscr{L}(\phi^\mu(x),\nabla_\alpha \phi^\mu(x)) $$
but the covariant derivative and the element $\sqrt{|g|}$ make it hard to isolate a surface term. I tried various way round, but I can't obtain the expected, seemingly covariant, E-L equation
$$ \frac{\partial \mathscr{L}}{\partial \phi^\mu} - \nabla_\alpha \frac{\partial \mathscr{L}}{\partial (\nabla_\alpha \phi^\mu)} $$
even because this very formula may be wrong; I don't know in fact what is the actual form of E-L equation to be found.
I've looked a lot on this site and various sources, but I can't find a step-by-step demonstration of E-L equations in curved spacetime and a clear final statement, so what is the path to follow that brings to the result?
 A: I’ll consider the case of a scalar field.
You just need to apply the usual method of deriving the Euler-Lagrange equations by the variational principle. You get:
$$
\frac{\partial \mathcal L}{\partial \phi}-\frac{1}{\sqrt g}\partial_k\left(\sqrt g \frac{\partial\mathcal L}{\partial \phi_{,k}}\right)=0
$$
You can prove that this equation is covariant. This is because $\frac{\partial\mathcal L}{\partial \phi_{,k}}e_k$ is a well defined vector field and the second term is its divergence. You can express it in terms of the covariant derivative in order to make the covariance more explicit. This is done by the following general formula of the divergence for any vector field $v$:
$$
\frac{1}{\sqrt g}\partial_k(\sqrt g v^k)=(\nabla_k v)^k
$$
using the identity: $\Gamma_{ij}^j=\frac{1}{\sqrt g}\partial_i \sqrt g$
You can therefore rewrite the equation as:
$$
\frac{\partial \mathcal L}{\partial \phi}-\nabla \cdot \frac{\partial\mathcal L}{\partial\nabla\phi}=0
$$
Hope this helps.
A: Let us consider the action
$$
S = \int_\mathcal{M} \text{d}^n x \sqrt{|g|} \mathscr{L}(\phi^\mu(x),\nabla_\alpha \phi^\mu(x))
$$
where $\mathcal{M}$ is a region of spacetime with a fixed background metric $g$, and vary $\phi^\mu$.  To do this, we take $\phi^\mu$ and replace it with $\phi^\mu + \delta \phi^\mu$, and then discard any terms that are higher order than $\mathcal{O}(\delta \phi^\mu)$ to find $\delta S$.  The requirement that $\delta S = 0$ for arbitrary $\delta \phi^\mu$ will then allow us to find the Euler-Lagrange equations.
The variation of $\delta \phi^\mu$ yields
$$
S + \delta S = \int_\mathcal{M} \text{d}^n x \sqrt{|g|} \mathscr{L}(\phi^\mu(x)+ \delta \phi^\mu,\nabla_\alpha \phi^\mu(x) + \nabla_\alpha \delta \phi^\mu)
$$
and so the first-order variation $\delta S$ can be seen to be
$$
\delta S = \int_\mathcal{M} \text{d}^n x \sqrt{|g|} \left[ \frac{\partial \mathscr{L}}{\partial \phi^\mu} \delta \phi^\mu + \frac{\partial \mathscr{L}}{\partial (\nabla_\alpha \phi^\mu)}  \nabla_\alpha \delta \phi^\mu \right]
$$
We now need to integrate by parts, as we would in flat spacetime.  In the context of a curved background metric, Gauss's theorem is
$$
\int_\mathcal{M} \text{d}^n x \sqrt{|g|} \nabla_\alpha A^\alpha = \oint_{\partial \mathcal{M}} \text{d}^{n-1} x \sqrt{|h|} n_\alpha A^\alpha
$$
where $h_{\mu \nu}$ is the induced metric on the boundary $\partial \mathcal{M}$ and $n_\alpha$ is the unit normal on $\partial \mathcal{M}$.1  (See Appendix B of Wald's General Relativity  or Chapter 3 of Poisson's A Relativist's Toolkit for a derivation.)  In the present case, this means that we have
$$
\delta S = \int_\mathcal{M} \text{d}^n x \sqrt{|g|} \left[ \frac{\partial \mathscr{L}}{\partial \phi^\mu}  - \nabla_\alpha \frac{\partial \mathscr{L}}{\partial (\nabla_\alpha \phi^\mu)}  \right] \delta \phi^\mu + \oint_{\partial \mathcal{M}} \text{d}^{n-1} x \sqrt{|h|} n_\alpha \frac{\partial \mathscr{L}}{\partial (\nabla_\alpha \phi^\mu)} \delta \phi^\mu.
$$
If we are considering field variations $\delta \phi^\mu$ that vanish on the boundary but are otherwise arbitrary,2 then the boundary term vanishes and the quantity in brackets in the bulk integrand must vanish, leaving us with the result
$$
\frac{\partial \mathscr{L}}{\partial \phi^\mu}  - \nabla_\alpha \frac{\partial \mathscr{L}}{\partial (\nabla_\alpha \phi^\mu)}  = 0
$$
as expected.

1  This implicitly assumes that $\partial \mathcal{M}$ is not a null hypersurface.  See Wald and/or Poisson for details on how such cases can be treated.
2  Alternately, if one wishes to consider a wider class of variations, one can instead modify the action by adding an explicit surface term which cancels this term.  This turns out to sometimes be helpful if the Lagrangian contains higher-order derivatives.  The best-known example of this is the Gibbons-Hawking-York boundary term, an addition to the Einstein-Hilbert action that has various nice properties that the "plain-vanilla" Einstein-Hilbert action does not.
