Your version of the stress tensor matches what's currently on the Wikipedia link under
Electromagnetic stress–energy tensor
https://en.wikipedia.org/wiki/Electromagnetic_stress%E2%80%93energy_tensor
except for having $μ_0$ in place of your $4π$:
$$T^{μν} = \frac{1}{μ_0}\left(F^{μα}{F^ν}_α - \frac{1}{4} g^{μν} F_{αβ}F^{αβ}\right),$$
and with your metric $g_{μν}$ being written there as the Minkowski metric $η_{μν}$. A close examination of the context indicates that it is using the conventions (which it didn't spell out, explicitly):
$$
\left(x^0, x^1, x^2, x^3\right) = (ct, x, y, z), \hspace 1em F_{0i} = \frac {E_i} c, \hspace 1em F_{ij} = B^k, \\
η_{00} = -1, \hspace 1em η_{0i} = 0 = η_{i0}, \hspace 1em η_{ij} = δ_{ij}
$$
as $i = 1, 2, 3$ and $(i,j,k)$ range cyclically over $(1,2,3)$, $(2,3,1)$ and $(3,1,2)$; and with $F_{μν} = -F_{νμ}$. This yields the components listed in the article:
$$
T^{00} = \frac 1 2 \left(ε_0 E^2 + \frac{B^2}{μ_0}\right) = \frac {ε_0} 2 \left(E^2 + B^2 c^2\right), \\
\left(T^{10}, T^{20}, T^{30}\right) = \frac{×}{μ_0 c} = ε_0c× = \left(T^{01}, T^{02}, T^{03}\right), \\
T^{ij} = -ε_0 E_i E_j - \frac{B_i B_j}{μ_0} + δ^{ij} \left(ε_0 E^2 + \frac{B^2}{μ_0}\right) = -ε_0 \left(E_iE_j + B^iB^jc^2 - \frac 1 2 δ^{ij} \left(E^2 + B^2 c^2\right)\right),
$$
where $ε_0 = 1/(c^2 μ_0)$. Under your convention, $ε_0 = 1/(8πc^2)$. The last term can be written in tensor-dyad form as:
$$-ε_0 \left( + c^2 - \frac 1 2 \left(E^2 + B^2 c^2\right)\right).$$
The authors' analysis is, at best, sloppy, if that's what they really wrote. They're saying, just look at the $$ and $$ parts separately:
$$-ε_0 \left( - \frac 1 2 E^2\right), \hspace 1em -\frac{1}{μ_0} \left( - \frac 1 2 B^2\right).$$
More genrally, consider for any vector $$, the corresponding expression $- + ½ V^2$. Suppose the vector is in the $x$ direction with $ = V$, where $(, , )$ will be used here to denote the unit vectors, respectively, for the $x$, $y$ and $z$ directions. Then
$$- + \frac 1 2 V^2 = -V^2 + \frac 1 2 (++)V^2 = \frac{V^2}{2} (-++).$$
The dyad has the following principal components: $-½ V^2$ in the direction parallel to $$ and $+½ V^2$ in all directions perpendicular to $$. Now, apply this separately to the $$ and $$ parts, above, to get the desired results. The eigenvalues for the $$ part will be $±B^2/(2μ_0)$ or, in your notation $±B^2/(8π)$, while the eigenvalues for the $$ part will be $±ε_0E^2/2$, or in your notation, $±E^2/(8πc^2)$.
In other words, they were taking eigenvectors/eigenvalues ... but only for the $$ and $$ parts separately. Instead, you should actually be looking at the eigenvalues and eigenvectors for the whole thing, not just for the parts. You should be perform the eigenvector/eigenvalue problem on either the 3×3 matrix of the spatial-coordinate components of the stress tensor, if you're looking for the decomposition of that, or else for the full 4×4 matrix of the stress tensor, itself.
For convenience, we'll write $ = /μ_0$. Consider instead the following transform:
$$\sqrt{ε_0} = \cos θ - \sin θ, \hspace 1em \sqrt{μ_0} = \cos θ + \sin θ$$
for some vectors $$, $$ and angle $θ$ later to be determined. Then:
$$
\frac{ε_0 E^2 + μ_0 H^2}2 = \frac{m^2 + n^2}2, \\
\frac{ε_0 E^2 - μ_0 H^2}2 = \frac{m^2 - n^2}2 \cos{2θ} - · \sin{2θ}, \\
\sqrt{ε_0 μ_0} · = · \cos{2θ} + \frac{m^2 - n^2}2 \sin{2θ}, \\
\sqrt{ε_0 μ_0} × = ×.
$$
Now impose the requirement that $$ and $$ be orthogonal: $· = 0$. Then
$$\sqrt{ε_0μ_0} · \cos{2θ} = \frac{ε_0 E^2 - μ_0 H^2}2 \sin{2θ}.$$
Therefore,
$$ \cos{2θ} = \frac{ε_0 E^2 - μ_0 H^2}2, \hspace 1em \sin{2θ} = \sqrt{ε_0 μ_0} ·,$$
for some factor $$ that's determined by the condition that ${\cos}^2{2θ} + {\sin}^2{2θ} = 1$. From this condition, we find:
$$^2 = {\left(\frac{ε_0 E^2 - μ_0 H^2}2\right)}^2 + {\left(\sqrt{ε_0μ_0} ·\right)}^2 = {\left(\frac{ε_0 E^2 + μ_0 H^2}2\right)}^2 - {\left|\sqrt{ε_0μ_0} ×\right|}^2.$$
Write
$$ = \frac{ε_0 E^2 + μ_0 H^2}2 = \frac{m^2 + n^2}2, \hspace 1em = \sqrt{ε_0μ_0} × = ×.$$
Then
$$ = \sqrt{^2 - ^2}.$$
Substituting into the equations for $m^2$ and $n^2$, we have:
$$\frac{m^2 + n^2}2 = , \hspace 1em \frac{m^2 - n^2}2 = .$$
Thus
$$m = \sqrt{ + }, \hspace 1em n = \sqrt{ - }.$$
Denote the unit vectors aligned with $$ and $$, respectively, $$ and $$, and define $ = ×$. Then, we have
$$ = \sqrt{ + } , \hspace 1em = \sqrt{ - } , \hspace 1em = × = .$$
In tensor-dyad form, we can write the components $\left(T^{ij}: i, j = 1, 2, 3\right)$ as:
$$-ε_0 - μ_0 + \frac{ε_0E^2 + μ_0H^2}2 = - - + ( + + ) ,$$
or, noting that
$$ + = ( + ) + ( - ) ,$$
as:
$$-ε_0 - μ_0 + \frac{ε_0E^2 + μ_0H^2}2 = - + - .$$
For the other components, we have:
$$T^{00} = -, \hspace 1em
\left(T^{10}, T^{20}, T^{30}\right) = \left(T^{01},T^{02},T^{03}\right) = × = .$$
For the null field (i.e. where $ε_0E^2 = μ_0H^2$ and $· = 0$), one has $ = 0$ and $ = $, this reduces further to:
$$T^{00} = -, \hspace 1em
\left(T^{10}, T^{20}, T^{30}\right) = \left(T^{01}, T^{02}, T^{03}\right) = , \hspace 1em \left(T^{ij}: i, j = 1, 2, 3\right) = - .$$
It has not passed my notice that this looks a lot like what you encounter in the quantized version of the field theory. Essentially, this is a directed continuum of photons, with energy density $$ oriented in the direction $$ - the classical version of a photon stream. The stress tensor has eigenvalues (meaning: zero pressure) for the transverse directions $$ and $$ and an $-$ eigenvalue for the longitudinal direction $$.
In your units $μ_0 = 4π$, $ε_0 = 1/{4πc^2}$, $ = μ_0 = 4π$, and:
$$ = \frac{ε_0 E^2 + μ_0 H^2}2 = \frac{E^2/c^2 + B^2}{8π}.$$
So, the eigenvalues are 0, for eigenvectors $$ and $$ (i.e. the directions along the plane formed by $$ and $$); and $-(E^2/c^2 + B^2)/{8π}$ for $$ (directions perpendicular to $$ and $$).
For non-null fields, the zero eigenvalues split into $-$ and $+$ along the directions of $$ and $$ respectively. These are axes rotated from the $(/c, )$ frame, with the relations:
$$ = \sqrt{ε_0} \cos θ + \sqrt{μ_0} \sin θ, \hspace 1em = \sqrt{μ_0} \cos θ - \sqrt{ε_0} \sin θ$$
with
$$\cos{2θ} = \frac{(ε_0 E^2 - μ_0 H^2)/2}{}, \hspace 1em \sin{2θ} = \frac{\sqrt{ε_0μ_0} ·}{}.$$
We can solve for $\cos θ$ and $\sin θ$:
$$\cos θ = ±\sqrt{\frac{1 + \cos{2θ}}2} = ±\sqrt{\frac{ + (ε_0 E^2 - μ_0 H^2)/2}{2}}, \\
\sin θ = \text{sgn}(±·) \sqrt{\frac{1 - \cos{2θ}}2} = \text{sgn}(±·)
\sqrt{\frac{ - (ε_0 E^2 - μ_0 H^2)/2}{2}}.$$
One can also carry out the eigenvalue/eigenvector analysis for the full 4×4 tensor, rather than just the space-like 3×3 sub-matrix, posing the problem as:
$$T^{μν} u_ν = λ g^{μν} u_ν$$
Then, the eigenvectors $$ and $$ become $(0,)$ and $(0,)$, again with respective eigenvalues $-$ and $+$. In place of the eigenvector $$ are two eigenvectors of the form $\left(u_0, \right)$, with:
$$- u_0 + = λ g^{00} u_0 = -λ, \hspace 1em u_0 - = λ g^{11} = λ ,$$
leading to the equations:
$$( - λ) u_0 = , \hspace 1em u_0 = + λ,$$
or
$$^2 - λ^2 = ^2 ⇒ λ = ±, \hspace 1em u_0 = \frac{ ± }{}.$$
The $$ eigenvector and $-$ eigenvalue is now $(( - )/, )$ with eigenvalue $-$, while we also have the fourth eigenvector $(( ± )/, )$ with eigenvalue $$.
Up to powers of $c$, the quantities $$, $$ and $$ are energy-density, momentum-density and mass-density. The null field, which is effectively a directed continuum of photons, therefore has zero mass density, in this sense of the term.
Parts of this analysis predates Einstein. He wasn't the first to pose the formula $E = mc^2$. Nothing in the above is actually specific to relativity, nor even to a vacuum. It can equally well be carried out in a non-relativistic setting, or for isotropic media other than a vacuum (either non-relativistic media, or even relativistic media). In that case, $ε_0$ and $μ_0$ generalize, respectively to $ε$ and $μ$, while $c = 1/\sqrt{ε_0μ_0}$ generalizes to the wave speed, denoted here by $V = 1/\sqrt{εμ}$. Making the various densities dimensionally consistent, $$ would be replaced by $V$, while $$ would be replaced by $V^2$, and their relation with the energy density would become $^2 - ^2 V^2 = ^2 V^4$.
The references of note are the following:
Poincaré's analysis (1900)
"La théorie de Lorentz la physique expérimentale et la physique mathématique”, Revue générale des sciences pures et appliquées 11: 1163-1175.
which I suspect might dovetail into the analysis just done here, and is almost certainly couched firmly in the non-relativistic paradigm (though I haven't looked at it closely yet), versus
Einstein's analysis (1905)
“Ist die Trägheit eines Körpers von seinem Energieinhalt abhängig?”, Annalen der Physik, 18, 639-641, September 1905.
which is more fundamental and bona fide relativistic.
