A photon gas obeys the equation of state $\rho=P/3$ and hence $T^{\mu}_{\quad\mu}=3P-\rho=0$.
(Can also be seen by expressing the stress energy tensor in terms of of the electromagnetic tensor as follows $T^{\mu \nu}=(F^{\mu \lambda}g_{\lambda \sigma }F^{\nu \sigma}-g^{\mu \nu} F^{\lambda \sigma}F_{\lambda \sigma})/4\pi $ and taking the trace.)
What is the physical meaning of this tensor being traceless?
To answer that would require first answering the question of what it would take for it to not be traceless, and we'll answer that by first revoking the assumption that the Lagrangian is the Maxwell-Lorentz Lagrangian and returning back to fundamentals, to proceed from there.
The most general Lagrangian density $(,)$ that one can construct from the fields $$ and $$, that is also Lorentz invariant, is one that is a function $(ℑ,)$ of the Lorentz invariants $ℑ = ½(E^2 - B^2c^2)$, and $ = ·$ of the fields.
The Maxwell equations $∇· = ρ$ and $∇× - ∂/∂t = $ are the Euler-Lagrange equations for the response fields, defined by $$ = \frac{∂}{∂}, \hspace 1em = -\frac{∂}{∂}.$$ For Lorentz-invariant Lagrangian densities, they can be expressed in terms of the derivatives with respect to the Lorentz invariants: $$ε = \frac{∂}{∂ℑ}, \hspace 1em θ = \frac{∂}{∂},$$ which might be termed the "response coefficients" or "constitutive coefficients" ($ε$ being none other than the permittivity), by $$ = ε + θ , \hspace 1em = ε c² - θ .$$
For Lagrangian densities of this form, the coefficients are generally functions $ε(ℑ,)$ and $θ(ℑ,)$ that satisfy the compatibility relation: $$\frac{∂ε}{∂} = \frac{∂θ}{∂ℑ},$$ and are otherwise arbitrary. For each Lagrangian density of this type, there is a corresponding set of constitutive coefficients, and each such set of functions defines a different Lagrangian and field theory; with the one corresponding to the Maxwell-Lorentz Lagrangian density given by $ε(ℑ,) = ε(0,0) ≠ 0$ and $θ(ℑ,) = θ(0,0) = 0$ as constants equal to their null field (i.e. (ℑ,) = (0,0)) values.
Under parity reversal $ℑ → ℑ$ and $ → -$, so a parity-symmetric Lagrangian density will have no linear or odd-order dependence on $$ ... which includes $θ = 0$ as a special case. Examples of parity-violating Lagrangian densities are those where the coefficient $θ$ is independent of $(ℑ,)$, but is a field in its own right (the axion, except authors seem to all be using the opposite sign for their $θ$'s).
Note that (barring the case where $θ$ a field in its own right) one can always set $θ(0,0) = 0$ by adjusting the definitions of response fields by $ → + θ(0,0) $ and $ → - θ(0,0) $. It won't affect the Euler-Lagrange equations, since one already has the Maxwell equations $∇· = 0$ and $∇× + ∂/∂t = $ for the field strengths $(,)$.
So, the only real empirical content of the Maxwell-Lorentz Lagrangian is that (1) $ε$ has a non-zero null-field limit and (2) the dynamics can be approximated in the presence of weak sources by their null-field limits. The assumption most definitely breaks down in the close vicinity of sources (e.g. "form factors" for fundamental particles) and Quantum Field Theory pays the price for it, if it tries to base quantization on the Maxwell-Lorentz Lagrangian. It is, to put it bluntly, the wrong classical limit to base quantization on. An example of a corrected classical limit is the Euler-Heisenberg Lagrangian density, which includes the second order terms of the Taylor expansion $(ℑ,)$ ... except there's no $ℑ$ term, since this would violate parity-symmetry.
The stress tensor, in its pre-geometric and "pre-metrical" form, is not a rank (2,0) or rank (0,2) tensor, but a rank (1,1) tensor density $^ρ_ν$ and is the one involved the continuity equation $∂_ρ ^ρ_ν = 0$, whose curved space-time form is $∂_ρ ^ρ_ν - Γ_{ρν}^μ ^ρ_μ = 0$. Here, and below I'm using the summation convention. The connection coefficients $Γ_{ρν}^μ$ are those associated with the metric $g_{μν}$ and the tensor that you're using to seeing $T_{μν}$ is metric-dependent and is related to the stress tensor density by $_{μν} ≡ g_{μρ} ^ρ_ν = \sqrt{|g|} T_{μν}$, where I also use the metric to define the tensor density version $_{μν}$ of it, and define $g$ (as usual) as the determinant of the matrix formed of the metric components $g_{μρ}$.
The stress tensor density is derived from the canonical stress tensor density of the Lagrangian density $$, which I'll denote $^ρ_ν$. For a first-order Lagrangian density $(q,v)$ that is a function of fields $q^A$ and their gradients $v^A_μ = ∂_μ q^A$, it is given by $$^ρ_ν = \frac{∂}{∂v^A_ρ} v^A_ν - δ^ρ_ν ,$$ where (as usual) $δ^ρ_ν$ is the Kronecker delta ($δ^ρ_ν = 1$ for $ρ = ν$ and $δ^ρ_ν = 0$ for $ρ ≠ ν$). The stress tensor density is derived from it by the inclusion of a gradient $$^ρ_ν = ^ρ_ν + ∂_μ ^{ρμ}_ν,$$ of a suitably-defined tensor density satisfying the condition $^{ρμ}_ν = -^{μρ}_ν.$
For the electromagnetic field, the role of the $q^A$'s is played by the potentials $A_μ$ (the $()_μ$ playing the role, here, of the $()^A$), with $A_0 = -φ$ and $ = \left(A_1, A_2, A_3\right)$, while the role of the gradients is (partially) played by the corresponding field strengths, which are given by $ = ∇×$ and $ = -∇φ - ∂/∂t$ or, in tensor form with $ = \left(F_{23}, F_{31}, F_{12}\right)$ and $ = \left(F_{10}, F_{20}, F_{30}\right)$ by $F_{μν} = ∂_μ A_ν - ∂_ν A_μ$, with $∂_0 = ∂/∂t$ and $\left(∂_1, ∂_2, ∂_3\right) = ∇$.
For a Lagrangian density $(A,F)$ whose dependence on the gradients of $A$ consists only of dependence on those contained in $F$, the response fields form the components $ = \left(^{01}, ^{02}, ^{03}\right)$ and $ = \left(^{23}, ^{31}, ^{12}\right)$ of the tensor density $$^{μν} = -\frac{∂}{∂F_{μν}}.$$
The relevant term $^{ρμ}_ν$ is that given by the Belinfante correction $^{ρμ}_ν = A_ν ^{ρμ}$. This will convert the explicitly gradient-dependent tensor density $^ρ_ν$ into one $^ρ_ν$ that has no explicit dependence on the potentials, and whose dependence on the gradients reduces to a dependence only on the field strengths and response fields (and the Lagrangian density). Explicitly, they are given by $$^ρ_ν = -^{ρμ} ∂_ν A_μ - δ^ρ_ν , \hspace 1em ^ρ_ν = ^ρ_ν + ∂_μ (A_ν ^{ρμ}) = ^{ρμ} F_{μν} - δ^ρ_ν ,$$ for a pure source-free field, i.e. one where $∂_μ ^{ρμ} = 0$.
In the presence of sources $∂_μ ^{ρμ} = ^ρ$, with $ρ = ^0$, and $ = \left(^1, ^2, ^3\right)$, you also have an extra contribution from $A_ν ^ρ$ to $^ρ_ν$. That's something J.D. Jackson includes and accounts for in his Electromagnetism monograph (https://en.wikipedia.org/wiki/Classical_Electrodynamics_(book)).
Interestingly (and something I just noticed), this does not exactly match the tensor density $_{μν}$ that arises from the equation for the Hilbert tensor density: $_{μν} = -2 ∂/∂g^{μν}$. The parity-violating invariant $$ has no metric-dependence. This can be seen directly when it is expressed in tensor form: $ = -⅛ ε^{μνρσ} F_{μν} F_{ρσ}$, where $ε^{μνρσ}$ is the anti-symmetric tensor density given by $ε^{0123} = 1$. So, all the $θ$ terms would drop out and the Hilbert tensor would be that which arose from the Belinfante-adjusted tensor density by taking $θ = 0$. Only the invariant $ℑ = -½ c √|g| g^{μρ} g^{νσ} F_{μν} F_{ρσ}$ has metric-dependence.
However, this won't have a much bearing on the answer this is all leading up to.
The tensor density derived from the Belinfante adjustment (again: for a pure source-free field) has the following components $$\begin{align} ^0_0 &= · - , & \left(^0_1, ^0_2, ^0_3\right) &= ×, \\ ^i_j &= D^i E_j + B^i H_j - δ^i_j \left(· + \right), & \left(^1_0, ^2_0, ^3_0\right) &= ×. \end{align}$$ The trace is $$^ρ_ρ = 2(· - ·) - 4 = 4(εℑ + θ - ) = 4\left(ℑ \frac{∂}{∂ℑ} + \frac{∂}{∂} - \right).$$ It expresses the extent to which $$ deviates from being homogeneous to the first degree in $ℑ$ and $$. In turn, given the relation to $ε$, it expresses the degree to which the field - with this Lagrangian - has a non-trivial permittivity; i.e. the degree to which the field, away from sources, behaves as a polarizable medium. It also expresses the degree to which the "axial" coefficient $θ$ is non-trivial.
For example, if $ = kℑ^{1+δ}$ for constant $k$, then $ε = kℑ^δ$ and $^ρ_ρ = 4δ$. So, an $$ with a fractal scaling dimension in $ℑ$ would produce a non-zero trace.
As we've already seen, for null fields, you already have $ε = ε(0,0)$ and (effectively) $θ = θ(0,0) = 0$. Null fields approximate the fields away from sources, where the only dominant contribution is the $1/r$ radiative part of the field. The expressions above, already, assume you're on the exterior of sources.
But when you get close to (meaning: at subatomic distances from) matter, things change. Quantum Field Theory is already giving you strong indication to that effect, which you have to take seriously even at the classical level; i.e. in the classical limit. There are literature references out there that treat the quantum field near particles as a polarizable medium and derive the relevant expressions for the constitutive coefficients, but I can't cite any off-hand. I might add them in a later edit (unless someone else adds them in first, hint, hint, editors).
