Derivation of electromagnetic stress energy tensor in curved spacetime I would like to know how to derive the Electromagnetic Stress-Energy Tensor in curved spacetime.
I would like to arrive at the result 
$$T^{\mu\nu} = \frac{1}{\mu_0} \left[ F^{\mu \alpha}F^\nu{}_{\alpha} - \frac{1}{4} \eta^{\mu\nu}F_{\alpha\beta} F^{\alpha\beta}\right] \,.$$
 A: This is my derivation
\begin{equation}
 T_{\mu\nu} = \frac{-2 c}{\sqrt{-g}}   \frac{\delta S_{M}}{\delta g^{\mu\nu}} \; .
\end{equation}
\begin{equation}
 S_{EM}[g^{\mu\nu},A^\mu] = \frac{-1}{4 \mu_0}\int d^4x \sqrt{-g} F_{\alpha\beta} F^{\alpha\beta} \; ,
\end{equation}
\begin{eqnarray}
  \delta_g S_{EM} &=& \frac{-1}{4 \mu_0}\int d^4x \bigg[ \delta_g(\sqrt{-g})F_{\alpha\beta} F^{\alpha\beta} + \sqrt{-g} \delta_g (F_{\alpha\beta} F^{\alpha\beta} ) \bigg]  \; ,\\
  &=&\frac{-1}{4 \mu_0}\int d^4x \bigg[ - \frac 1 2 \sqrt{-g} g_{\mu\nu} \delta g^{\mu\nu}   F_{\alpha\beta} F^{\alpha\beta} +  \sqrt{-g} \delta_g (F_{\alpha\beta} F^{\alpha\beta} ) \bigg]\; .\\
\frac{\delta S_{EM}}{\delta g^{\mu\nu}}  &=&   \frac{-1}{4 \mu_0} \bigg[ - \frac 1 2 \sqrt{-g} g_{\mu\nu}   F_{\alpha\beta} F^{\alpha\beta} +  \sqrt{-g} \frac{\delta}{\delta g^{\mu\nu}} (F_{\alpha\beta} F^{\alpha\beta} ) \bigg]\;  .   (1)
\end{eqnarray}
Consider the last term in the vielbein form
\begin{equation}
 F_{\alpha\beta} F^{\alpha\beta} = e^I_\alpha e^J_\beta F_{IJ} e^\alpha_K e^\beta_L F^{KL}\; .(2)
\end{equation}
We have done here for isolated the flat structure ($g_{\mu\nu}$-independent) from curved structure. Next, we will use the chain rule
\begin{equation}
 \frac{\delta\,\,}{\delta g^{\mu\nu}} = \frac{\delta\,\,}{\delta e^\lambda_P} \;  \frac{\delta e^\lambda_P\,}{\delta g^{\mu\nu}}\;. (3)
\end{equation}
From $g^{\mu\nu} = \eta^{MP}e^\mu_M e^\nu_P \;$ we have
\begin{equation}
 \delta g^{\mu\nu} = 2 \eta^{MP}e^\mu_M \delta^\nu_\lambda \, \delta e^\lambda_P\; .(4)
\end{equation}
By using (2), (3) and (4) we can calculating the last term of (1)as
\begin{eqnarray}
 \frac{\delta (F_{\alpha\beta} F^{\alpha\beta} )}{\delta g^{\mu\nu}} &=& \frac{\delta}{\delta e^\lambda_P} (e^I_\alpha e^J_\beta F_{IJ} e^\alpha_K e^\beta_L F^{KL} ) \; \frac{\delta e^\lambda_P\,}{\delta g^{\mu\nu}}\;\\
 &=& 4 e^I_\alpha e^J_\beta e^\alpha_K \frac{\delta e^\beta_L}{\delta e^\lambda_P} F_{IJ} F^{KL} \;\frac{\delta e^\lambda_P\,}{\delta g^{\mu\nu}}\;\\
 &=& ( 4 e^I_\alpha e^J_\beta e^\alpha_K  \delta^\beta_\lambda \delta^P_L F_{IJ} F^{KL} )(\frac 1 2  \eta_{MP} e^M_\mu \delta^\lambda_\nu  )\;,\\
 &=& 2 e^I_\alpha e^J_\beta e^\alpha_K \delta^\beta_\lambda \delta^P_L \delta^\lambda_\nu \, e^M_\mu \eta_{MP} \, F_{IJ}F^{KL}\;,\\
 &=& 2 e^I_\alpha e^J_\nu e^\alpha_K e_{L \mu } F_{IJ} F^{KL}\;,\\
 &=& 2 F_{\alpha \nu} F^\alpha {}_\mu = 2 g^{\alpha\beta} F_{\alpha \mu} F_{\beta \nu} \; .
\end{eqnarray}
Then we obtain
\begin{equation}
 \frac{\delta S_{EM}}{\delta g^{\mu\nu}}  =  \frac{1}{8 \mu_0} \sqrt{-g} g_{\mu\nu}   F_{\alpha\beta} F^{\alpha\beta} - \frac{1}{4 \mu_0} \sqrt{-g} (2 g^{\alpha\beta} F_{\alpha\mu} F_{\beta\nu})\;,
\end{equation}
and the energy-momentum tensor of the electtromagnetic field reads
\begin{equation}
 T_{\mu\nu} = \frac{-2 c}{\sqrt{-g}}   \frac{\delta S_{EM}}{\delta g^{\mu\nu}} = \frac c {\mu_0}g^{\alpha\beta} F_{\alpha\mu} F_{\beta\nu} -\frac {c} {4 \mu_0} g_{\mu\nu}   F_{\alpha\beta} F^{\alpha\beta} \; . \label{TEM}
\end{equation}
A: Start with Hamiltonian density, the quantity in the integrand from the definition of the Hamiltonian:
$$H = \int d^3x \left( \psi_{,0} \frac{\partial \mathcal{L}}{\partial \psi_{,0}} - \mathcal{L} \right) \equiv \int d^3x \,\mathcal{H}$$
$\mathcal{L}$, of course, denotes the  Lagrangian density.
Since $\mathcal{H}$ corresponds to Hamiltonian density, it should be the $(00)$ component of the energy-momentum tensor, i.e. $T_{0}^{0}. $Instead of some generic field $\psi$, plug in the photon field $A_\mu$ and upgrade the equation to the full covariant form:
$$T_{\mu}^{\nu} = A_{\sigma ,\mu} \frac{\partial \mathcal{L}}{\partial A_{\sigma ,\nu}} - \delta^{\nu}_{\mu} \mathcal{L}_{ED}$$
Since $$\mathcal{L}_{ED} = -\frac{1}{4}F^{\mu \nu} F_{\mu \nu}$$
$$F^{\mu \nu} = A^{\nu,\mu} - A^{\mu,\nu}$$
a straightforward calculation gives
$$T_{\mu}^{\, \, \nu} = \frac{1}{4} \left( -A^{\sigma}_{,\mu}F^{\nu}_{\,\, \sigma} + \frac{1}{4} \delta^{\nu}_{\mu} F^{\sigma \rho} F_{\sigma \rho} \right).$$
This is the canonical energy-momentum tensor, which is generally not symmetric nor gauge invariant. To fix that, you simply make the tensor symmetric by adding a suitable (basically irrelevant) term $S_{\sigma \mu \nu}$ such that:
$$S_{\sigma \mu \nu} = - S_{\mu \sigma  \nu}$$
$$\bar{T}_{\mu}^{\, \, \nu} = \bar{T}_{\nu}^{\, \, \mu} = T_{\mu}^{\, \, \nu} + \partial^\sigma S_{\sigma \mu \nu}$$
which should give you the right expression. Incidentally, the term you should get is $S_{\mu \nu \sigma} =  A_\sigma F_{\mu\nu}$. You can simply plug it in as you would do with any other ansantz and see what it does.
I haven't done the complete calculation myself for the purposes of writing this answer, so I might be off by a minus sign, multiplicative constant or up to a permutation/relabeling of indices. Tell me if you see something wrong (or edit the answer yourself). But this should be enough to give you an idea how to derive it.
A: Saksith Jaksri's answer is correct and there is nothing wrong with using vielbeins to calculate $\delta(F^{\alpha\beta}F_{\alpha\beta})$. But in this case, a more direct calculation is actually easier.
$$
\begin{equation} \label{eq1}
\begin{split}
&& \frac{\delta \left(g^{\alpha\gamma} g^{\beta\delta} F_{\alpha\beta} F_{\gamma\delta}\right)}{\delta g^{\mu\nu}}  \\
&=& g^{\beta\delta} F_{\mu\beta} F_{\nu\delta} + g^{\alpha\gamma} F_{\alpha\mu} F_{\gamma\nu}  \\
&=& 2 {F_{\mu}}^\rho F_{\nu\rho}.
\end{split}
\end{equation}
$$
