Why is $\nabla_\mu T^{\mu\nu} = 0$ valid for $f(R)$ theories? My question is in the context of $f(R)$ theories,  assuming that in the early universe $f(R)\sim R^n$ ($n>1$). 
Why is the co-variant conservation law for the energy-momentum tensor found in General Relativity still valid?
 A: Diffeomorphism invariance of the matter action $S_m$ leads (via Noether's 2nd theorem) to the identity$^1$
$$ \nabla_{\mu} T^{\mu\nu}~\stackrel{m}{\approx}~0, \qquad  T^{\mu\nu}~:=~\mp\frac{2}{\sqrt{|g|}}\frac{\delta S_m}{\delta g_{\mu\nu}},
\tag{1} $$
cf. e.g. Ref. 1. [Here the $\stackrel{m}{\approx}$ symbol means equality modulo matter eoms. The connection $\nabla$ is the Levi-Civita connection. The Minkowski sign convention is $(\pm,\mp,\mp,\mp)$.]
References:


*

*R.M. Wald, GR; Appendix E.1.


--
$^1$ Note that eq. (1) is not a conservation law by itself. To get a conservation law, we need a Killing vector field, cf. e.g. my Phys.SE answer here.
A: Consider you get the field equations by setting $\delta S=0$, where $S$ is the action taken. 
Remember that for f(R) theories, the action goes:
$S=\dfrac{16\pi G}{c^4}\int\,d^4\mathbf{x}\,\left(\sqrt{-g}\,f(R)+\mathcal{L}_m\right)$
and original, Einstein field equations are only modified in the left side of the equality ("geometric terms"):
$\dfrac{\partial f(R)}{\partial R}R_{\mu\nu}-\dfrac{f(R)}{2}g_{\mu\nu}+\left[g_{\mu\nu}g^{\alpha\beta}\nabla_{\alpha}\nabla_{\beta}-\nabla_{\mu}\nabla_{\nu}\right]\dfrac{\partial f(R)}{\partial R}=\dfrac{8\pi G}{c^4}T_{\mu\nu}$.
For General Relativity, the Einstein-Hilbert action goes:
$S=\dfrac{16\pi G}{c^4}\int\,d^4\mathbf{x}\,\left(\sqrt{-g}\,R+\mathcal{L}_m\right)$
and thus one get famous, well-known, Einstein field equations:
$R_{\mu\nu}-\dfrac{R}{2}g_{\mu\nu}=\dfrac{8\pi G}{c^4}T_{\mu\nu}$
where as you can see, we don't get yet into the cosmological constant (or the vacuum energy density).
You can take the definition of the stress-energy tensor as:
$T_{\mu\nu}=\,-\dfrac{2}{\sqrt{-g}}\dfrac{\,\delta\,\left(\sqrt{-g}\mathcal{L}_m\right)}{\delta\,g^{\mu\nu}}$
for any action given where $\mathcal{L}_m$ is the lagrangian of matter present in your spacetime configuration, corresponding in the end in the stress-energy tensor. 
From above you can see that the "energy-momentum terms" on the right side of the equations are the same for both theories and are added manually in the actions. Typically, we get Einstein-field equations taking E-H action negliging $\mathcal{L}_m$ for vacuum:
$R_{\mu\nu}-\dfrac{R}{2}g_{\mu\nu}=0\,$   for   $\,\mathcal{L}_m=0$
The conservation law will be still obeyed from the definition of the action for f(R).
