I'm still not sure I see your reasoning for wanting to change the form, can you clarify your question.

I'll spit something out that has to do with energy density for the photon gas in hopes that it helps you out/ is right.

The approximate total number of photons is given by
[tex]
\langle N \rangle =\sum_{\omega} \langle n_\omega \rangle
[/tex]
with [itex]\langle n \rangle [/itex] as the approximate number of photons in a given frequency mode. Now if we take a continuum limit for the frequencies we write
[tex]
\langle N \rangle =\int_{0}^{\infty} g(\omega) \langle n \rangle d\omega
[/tex]
Now for bosons,
[tex]
\langle n \rangle =(\exp(\hbar \omega /kT)-1)^{-1}
[/tex]
For the density of states [itex]g(\omega)[/itex]
[tex]
g(\omega)d\omega=\frac{d^{3}p\, d^{3}x}{\hbar^3}=\frac{V}{\hbar^3}d^3 p
[/tex]
Now expressing the momentum in terms of frequency we have
[tex]
p=\frac{\hbar \omega}{c}\implies dp=\frac{\hbar}{c} d\omega
[/tex]
Now
[tex]
d^3 p=4\pi p^2 dp=4\pi \left(\frac{\hbar}{c}\right)^3 \omega^2 d\omega \implies g(\omega)d\omega=4\pi \frac{V}{\hbar^3}\left(\frac{\hbar}{c}\right)^3 \omega^2 d\omega
[/tex]
Now we have to remember to throw in the number 2 cause our math won't do the counting for us, there are two transverse modes to account for.
[tex]
g(\omega)d\omega=\frac{8 \pi V}{c^3}\omega^2 d\omega
[/tex]
then
[tex]
\langle N \rangle=\frac{8 \pi V}{c^3}\int_{0}^{\infty}\frac{\omega^2 d\omega}{e^{\frac{\hbar\omega}{kT}}-1}
[/tex]
to get [itex]\langle E \rangle[/itex] we just need to calculate
[tex]
\langle E \rangle =\sum_{\omega}\langle n_\omega \rangle \epsilon_\omega
[/tex]
with [itex]\epsilon=\hbar \omega [/itex].