Why Photon gas's Equation of State Diverges? We know that, for Photon gas, the equation of state is given by $pV=\frac{1}{3}U$; where $p$ is the Pressure, $V$ is the Volume and $U$ is the Internal Energy of the Photon Gas. (see Equation of State of Photon gas (as done with Grand Canonical Partition Function & Grand Canonical Potential Function))
For Nonrelativistic Bosons (e.g. Higgs Boson), the Equation of State can be shown as $pV=\frac{2}{3}U$; where $p$ is the Pressure, $V$ is the Volume and $U$ is the Internal Energy of the Nonrelativistic Boson Gas (by considering $\varepsilon=P^2/2m$; where $\varepsilon$ is Kinetic Energy, $P$ is Linear Momentum & $m$ is Mass of Nonrelativistic Boson). (see Equation of State of Nonrelativistic Fermion & Boson gas (as done with Average Number of Nonrelativistic Bosons in Energy state $\varepsilon\ {\rm and}\ \varepsilon+d\varepsilon$))
But, we know that Photon is Ideal Relativistic Boson. So, if we apply the Relativistic Approach i.e. $\varepsilon=Pc$; where $c$ is Speed of Light, on the Average Number of Photons (Relativistic Bosons), we should get the Equation of State $pV=\frac{1}{3}U$ for Relativistic Bosons. But, the problem there is, it will be true if and only if it could be shown that
$$I=\int_0^\infty \frac{d \varepsilon}{e^{\varepsilon/k T}-1}=0,$$
but actually it is $\infty$ as shown in the Calculation.
Then, how can we show that, for Photon gas $pV=\frac{1}{3}U$ in this way?
 A: I couldn't quite follow your calculations for $U$, especially when the bracket appears line $3$. Normally you do a simple integration by parts to get the state equation. You have for non-interacting bosons, with density of state $D$ :
$$
\Omega = -T\int d\epsilon D(\epsilon)(-\ln(1-e^{-\beta(\epsilon-\mu)}))
$$
$$
U = \int d\epsilon D(\epsilon)\frac{\epsilon}{e^{\beta(\epsilon-\mu)}-1}
$$
So it is tempting to relate the two by an integration by parts. This is possible when in general: $D(\epsilon)\propto 1_{\mathbb R_+}(\epsilon)\epsilon^\alpha$ (for 3D photons, $\alpha=2$ and $\mu=0$). Both integrals are well-defined as long as $\alpha>-1$ (for low energy limit, high energy is exponentially supressed). You therefore get:
$$
U = \left[-TD(\epsilon)\epsilon(-\ln(1-e^{-\beta(\epsilon-\mu)})\right]_0^\infty-(\alpha+1)\Omega
$$
and the bracket is trivial since $\alpha>-1$, so the power law overcomes the $\ln$ in the low energy limit, and you still have exponential suppression in the high energy limit. If you further assume $D(\epsilon)\propto V$, then
$$
\Omega=-pV
$$
and you recover:
$$
U = \frac{1}{1+\alpha}pV
$$
and as a sanity check, you can see that the condition $\alpha>-1$ is important for the final result to make sense.
Hope this helps and tell me if you need more details.
