Solution (source code)

= Solution

Use the <electromagnetic field tensor> convention stated above, and lower indices with the <Minkowski metric> $g=\operatorname{diag}(-1,1,1,1)$. Lowering a single time index changes its sign. Summing both orders of each antisymmetric index pair gives the <Lorentz scalar>
$$
F_{\kappa\lambda}F^{\kappa\lambda}
=2\sum_iF_{0i}F^{0i}+\sum_{i,j}F_{ij}F^{ij}
=2(B^2-E^2/c^2).
$$
For $T^{00}$, $F_j{}^0=F^{j0}=-E_j/c$ and $g^{00}=-1$. Direct substitution in the printed <electromagnetic stress-energy tensor> yields
$$
\begin{aligned}
T^{00}&=\epsilon_0c^2\left[-\sum_j(E_j/c)(-E_j/c)
+\frac14\,2(B^2-E^2/c^2)\right]\\
&=\frac{\epsilon_0}{2}(E^2+c^2B^2).
\end{aligned}
$$
For a spatial index $i$, $g^{0i}=0$ and $F_j{}^i=F^{ji}=-\epsilon_{ijk}B_k$. Hence
$$
T^{0i}=-\epsilon_0c^2\sum_jF^{0j}F_j{}^i
=\epsilon_0c\epsilon_{ijk}E_jB_k.
$$
Thus the full top row is
$$
\boxed{T^{0\nu}=\left(u_{\rm EM},\frac{S_x}c,\frac{S_y}c,\frac{S_z}c\right),\qquad
u_{\rm EM}=\frac12\left(\epsilon_0E^2+\frac{B^2}{\mu_0}\right),\quad
\mathbf S=\frac1{\mu_0}\mathbf E\times\mathbf B.}
$$
Here $u_{\rm EM}$ is the <electromagnetic energy density> and $\mathbf S$ the <Poynting vector>, using $\epsilon_0\mu_0c^2=1$. By the symmetry $T^{0i}=T^{i0}$, these components describe both <energy> flux divided by $c$ and $c$ times the <electromagnetic momentum density> $\mathbf g=\mathbf S/c^2$. Keeping the factor $c$ is essential because the time coordinate is $ct$.