Conservation of stress-energy and the fluid equation in cosmology I'm trying to derive the equation for the cosmological fluid:
$$\dot \rho + 3 \frac{\dot a}{a}(\rho +P)=0$$
by starting from the conservation of the stress-energy tensor:
$$\nabla^\mu T_{\mu \nu} = 0$$
with the stress-energy for a perfect fluid in its own frame being:
$$ T_{\mu \nu} = \text{diag} (\rho, a(t)^2 P,a(t)^2 P,a(t)^2 P) $$
in a spatially flat FLRW metric:
$$g_{\mu \nu} = \text{diag}(1,-a(t)^2,-a(t)^2,-a(t)^2)$$
But I keep getting a bogus answer! Consider the equation you get from $ \nabla^\mu T_{\mu \nu} = 0$ when $\nu =0$:
$$
\begin{align*}
\nabla^\mu T_{\mu 0} &= 0 \\
g^{\mu \alpha}\nabla_\alpha T_{\mu 0} &= 0
\end{align*}
$$
$T$ is diagonal, so $\mu$ must be zero, but $g$ is diagonal as well, so if $\mu$ is zero, then so is $\alpha$. This gives:
$$
\begin{align*}
g^{0 0}\nabla_0 T_{0 0} &= 0 \\
\nabla_0 \rho &= 0 \\
\dot \rho &= 0
\end{align*}
$$
Because $\rho$ is just a scalar, so the covariant derivative is the partial derivative. Except this answer is wrong.
 A: It is easier to make this calculation covariantly. 
Given a time-like vector field $n^\mu$ representing the direction normal to the hyper-surfaces (in the Friedmann case the homogeneous and isotropic hyper-surfaces), you can always decompose the energy momentum tensor $T_{\mu\nu}$ as $$T_{\mu\nu} = \rho n_\mu n_\nu + 2n_{(\mu}q_{\nu)} + p \gamma_{\mu\nu}+ \Pi_{\mu\nu},$$ where $\rho \equiv T_{\mu\nu}n^\mu n^\nu$ is the energy density, $q_\alpha \equiv - T_{\mu\nu} n^\mu \gamma^\nu{}_\alpha$ the energy (heat) flow, $p \equiv T_{\mu\nu}\gamma^{\mu\nu} / 3$ the pressure and $\Pi_{\mu\nu} \equiv T_{\alpha\beta}\gamma^{\alpha}{}_\mu\gamma^{\beta}{}_{\nu} - p\gamma_{\mu\nu}$ all measured in the hyper-surfaces described by $n^\mu$, $\gamma_{\mu\nu} \equiv g_{\mu\nu} + n_\mu n_\nu$ is the spatial projector/metric and the signature is $(-1,1,1,1)$, i.e., $n_\mu n^\mu = -1$. 
The energy momentum tensor "conservation" gives $\nabla^\mu T_{\mu\nu} = 0$. To project this result in the time direction (that is want you want when you put $\nu = 0$) you just multiply this expression by $n^\nu$, i.e., 
$$n^\nu\nabla^\mu T_{\mu\nu} = \nabla^\mu (T_{\mu\nu}n^\nu) - T_{\mu\nu}\nabla^\mu n^\nu = \nabla^\mu(-\rho n_\mu - q_\mu) - \mathcal{K}^{\mu\nu}T_{\mu\nu},$$ where $\mathcal{K}_{\mu\nu} \equiv \nabla_\mu n_\nu$ is the extrinsic curvature.
Now for a Friedmann metric the extrinsic curvature is simply $\mathcal{K}_{\mu\nu} = \Theta \gamma_{\mu\nu} / 3$ where $\Theta$ is the expansion factor and the heat flow $q_\mu$ is zero, hence, we obtain
$$n^\nu\nabla^\mu T_{\mu\nu} = -\left[\dot{\rho}+\Theta(\rho + p) + \nabla_\mu q^\mu + \sigma^{\mu\nu}\Pi_{\mu\nu}\right],$$ where $\dot{\rho} \equiv n^\mu\nabla_\mu \rho$ and $\sigma^{\mu\nu}$ is the traceless symmetric part of the extrinsic curvature which is called shear tensor. At this point one can also introduces a scale factor satisfying $3\dot{a}/a = \Theta$.
The key ingredient to make this calculation correctly is the explicit introduction of the time vector $n^\mu$, note that it is its commutation with the covariant derivative and its derivative in terms of the expansion factor that introduces the missing terms in your result. If you want to make this calculation in the old fashion coordinate way you can just calculate everything in a specific coordinate system with all partial derivatives, Christoffel symbols, etc, and then project in the chosen direction.
A: Consider the (1,1)-tensor $T^{\mu}_{\phantom{\mu}\nu}=g^{\mu\lambda}T_{\lambda\nu}=\text{diag}(\rho, -P/c^2, -P/c^2, -P/c^2)$. Then
$$
\nabla_\mu T^{\mu}_{\phantom{\mu}\nu} = \partial_\mu T^{\mu}_{\phantom{\mu}\nu} +
\Gamma^{\mu}_{\mu\lambda}T^\lambda_{\phantom{\mu}\nu} - \Gamma^{\lambda}_{\mu\nu}T^\mu_{\phantom{\mu}\lambda} = 0.
$$
For $\nu=0$, we find
$$
\nabla_\mu T^{\mu}_{\phantom{\mu}0} = \partial_0 T^{0}_{\phantom{\mu}0} +
\Gamma^{\mu}_{\mu 0}T^0_{\phantom{\mu}0} - \Gamma^{\lambda}_{\mu 0}T^\mu_{\phantom{\mu}\lambda} = 0.
$$
Now, the only non-zero Christoffel symbols in this equation are
$$
\Gamma^{1}_{10} = \frac{1}{2}g^{11}\frac{\partial g_{11}}{\partial x^{0}} = \frac{\dot{a}}{ca} = \Gamma^{2}_{20} = \Gamma^{3}_{30},
$$
so that
$$
\nabla_\mu T^{\mu}_{\phantom{\mu}0} = \frac{\dot{\rho}}{c} + \frac{3\dot{a}}{ca}\rho + \frac{3\dot{a}}{ca}\frac{P}{c^2} = 0,
$$
which is indeed the conservation law
$$
\dot{\rho} + 3\frac{\dot{a}}{a}\left(\rho + \frac{P}{c^2}\right)=0.
$$
You can also verify that $\nabla_\mu T^{\mu}_{\phantom{\mu}\nu} \equiv0$ for $\nu\neq0$.
