Energy-Momentum Tensor for Electromagnetism in Curved Space $\newcommand{\l}{\mathcal L} \newcommand{\g}{\sqrt{-g}}$$\newcommand{\fdv}[2]{\frac{\delta #1}{\delta #2}}$I want to calculate the energy-momentum tensor in curved free space by functional differentiation with respect to the metric. The Lagrangian density that I have with units $c=1$ is the following
$$\l = \g \left(- \frac 1 4 F^{\mu \nu} F_{\mu \nu}  \right).$$
I have calculated the variation of the action $\delta S$ to be
$$\delta S = \frac 1 2 \int\g \left(\frac 1 4 F^{\mu \nu} F_{\mu \nu} g_{\rho \sigma}  - F^\tau_{\;\;\rho} F_{\tau \sigma} \right) \delta g^{\rho \sigma} - \g F^{\mu \nu} (\delta F_{\mu \nu})\; \mathrm d^4 x. \tag 1  $$
However I cannot get rid of the $\delta F_{\mu \nu}$ term. I've looked Wikipedia and seen that the energy momentum tensor with lower indices is the term in the brackets, which makes me think that the following association is correct
$$F^{\mu \nu} (\delta F_{\mu \nu}) \leftrightarrow T_{\rho \sigma} \delta g^{\rho\sigma} .\tag 2 $$
Using the association $(2)$, I can write equation $(1)$ as 
$$\delta S = \frac 1 2 \int\g \left[ \left(\frac 1 4 F^{\mu \nu} F_{\mu \nu} g_{\rho \sigma}  - F^\tau_{\;\;\rho} F_{\tau \sigma} \right)  - T_{\rho \sigma}  \right] \delta g^{\rho\sigma} \; \mathrm d^4 x \quad ,$$
so that 
$$\frac 2 \g \fdv \l {g^{\rho \sigma}} = \left(\frac 1 4 F^{\mu \nu} F_{\mu \nu} g_{\rho \sigma}  - F^\tau_{\;\;\rho} F_{\tau \sigma} \right) - T_{\rho \sigma} = 0 $$
$$ \implies T_{\rho \sigma} = \frac 1 4 F^{\mu \nu} F_{\mu \nu} g_{\rho \sigma}  - F^\tau_{\;\;\rho} F_{\tau \sigma}  \quad , $$
which is the correct result but it is unclear to me as to why I the association $(2)$ should be true.
 A: The energy momentum tensor is found by varying the metric and holding all other fields constant. Since clearly $$\frac{\partial F}{\partial g}=0\longleftrightarrow
\delta_gF=0$$
we end up with
$$\delta_g S=\frac{1}{2}\int\mathrm{d}v\,\left(F^2g_{\mu\nu}/4-F^\tau{}_\mu F_{\tau\nu}\right)\delta g^{\mu\nu}$$
and comparison with $$\delta_gS:=\frac{1}{2}\int\mathrm{d}v\,T_{\mu\nu}\delta g^{\mu\nu}$$
leads to the correct result.
Note that what we are really varying here is just the matter action. The full general relativistic action contains the Einstein-Hilbert gravitational action, the cosmological constant term and a matter term. Put together, we have
$$ S_\mathrm{GR}=S_\mathrm{EH}+S_\Lambda+S_\mathrm{M}$$
It is true that the variation of this vanishes wrt. the metric, i.e. $\delta_g S_\mathrm{GR}=0$. However, the energy-momentum tensor is calculated by examining only $\delta_g S_\text{M}$, which is not generally zero.
A: The Einstein-Hilbert lagrangian coupled to a matter action 
$$
S_m[\varphi,g] = \int d^Dx\, \sqrt{-g}\mathcal L_m(\varphi,\partial_\mu\varphi), 
$$
i.e.
$$
S[g,\varphi] = \frac{1}{16\pi G}\int d^Dx\,\sqrt{-g} R + \int d^Dx\,\sqrt{-g} \mathcal L_m(\varphi,\partial_\mu\varphi),
$$
satisfies
$$
\delta S=-\frac{1}{16\pi G} \int d^Dx\, \sqrt{-g} \mathcal E^{\mu\nu}\delta g_{\mu\nu} 
+ \int d^Dx\, \sqrt{-g} \frac{1}{\sqrt{-g}}\frac{\delta\left( \sqrt{-g}\mathcal  L_m\right)}{\delta g_{\mu\nu}}\delta g_{\mu\nu}
$$
where $\mathcal E_{\mu\nu}$ is the Einstein tensor $R_{\mu\nu}-g_{\mu\nu}R/2$ and hence we define the Rosenfeld energy-momentum tensor
$$\boxed{
T^{\mu\nu} = \frac{2}{\sqrt{-g}}\frac{\delta\left( \sqrt{-g}\mathcal  L_m\right)}{\delta g_{\mu\nu}}=
\frac{2}{\sqrt{-g}}\frac{\delta S_m[\varphi, g]}{\delta g_{\mu\nu}}}
$$
getting
$$
dS=-\frac{1}{16\pi G}\int d^Dx\,\sqrt{-g}\left(\mathcal E^{\mu\nu}- 8\pi G T^{\mu\nu} \right)\delta g_{\mu\nu}.
$$
The action of a massless vector $A_\mu$, where $F_{\mu\nu}=\partial_\mu A_\nu-\partial_\nu A_\mu=\nabla_\mu A_\nu-\nabla\nu A_\mu$, can be written as
\begin{align*}
S_m =&\ -\frac{1}{4}\int d^Dx\, \sqrt{-g}F_{\mu\nu}F^{\mu\nu} = 
-\frac{1}{2} \int d^Dx\, \sqrt{-g}  A_\mu \left(-g^{\mu\nu}\Box+\nabla^{\nu}\nabla^\mu\right)A_\nu.
\end{align*}
So the equations of motion for $A_\mu$ are
$$
\left(-g^{\mu\nu}\Box+\nabla^{\nu}\nabla^\mu\right)A_\nu=0.
$$
The energy-momentum tensor is easier to express from the first form of the action:
$$
\delta S_m = - \frac{1}{4}\int d^Dx\, \frac{1}{2} \sqrt{-g} g^{\mu\nu}\delta g_{\mu\nu} F_{\mu\nu}F_{\alpha\beta}g^{\mu\alpha}g^{\nu\beta}  
+\frac{1}{2}\int d^Dx\, \sqrt{-g} F_{\mu\nu}F_{\alpha\beta}g^{\mu\alpha}g^{\nu\sigma} g^{\beta\tau}\delta g_{\sigma\tau}
$$
from which
$$\boxed{
T^{\mu\nu} = F^{\alpha\mu} F^{\beta\nu}g_{\alpha\beta}-\frac{1}{4}g^{\mu\nu}F_{\rho\sigma}F^{\rho\sigma}}=F^{\alpha\mu} F^{\beta\nu}g_{\alpha\beta}+g^{\mu\nu} \mathcal L_m,
$$
and the Einstein equations
$$
R^{\mu\nu}-\frac{1}{2}g^{\mu\nu}R+\Lambda g^{\mu\nu} =8\pi G T^{\mu\nu}. 
$$
A: Comment to the question (v2): The association (2) is not correct. To find the Hilbert SEM tensor, one varies the action wrt. the metric $g_{\mu\nu}$; not wrt. the gauge potential $A_{\mu}$ (or the field strength $F_{\mu\nu}$). 
