This is a way to prove convergence and also to understand why this integral represents potential energy between two charge distributions. Another way to look at this is
$$\iiint \frac{1}{2} \epsilon_0(\vec{E}_1+\vec{E}_2)^2d\tau $$ : $$\iiint \frac{1}{2} \epsilon_0(\vec{E}_1)^2d\tau +$$
$$\iiint \frac{1}{2} \epsilon_0(\vec{E}_2)^2d\tau +$$
$$\iiint \epsilon_0\vec{E}_1\cdot\vec{E}_2\,d\tau$$ The last component represents PD between charge distributions, $$=\iiint \epsilon_0\vec{E}_1\cdot(-\nabla V_{2})\,d\tau.$$ Rearranging the vector field identity $$\nabla \cdot (V_{2}\vec{E}_1) = V_{2} \nabla \cdot \vec{E}_1 + \vec{E}_{1} \cdot \nabla V_{2},$$ gives $$\vec{E}_{1} \cdot (-\nabla V_{2})= -\nabla\cdot(V_{2}\vec{E}_1) + V_{2} \nabla \cdot\vec{E}_1.$$ Substituting this in the integral gives $$\iiint\epsilon_0\left[-\nabla \cdot (V_{2}\vec{E}_1) + V_{2} \nabla \cdot \vec{E}_1\right]d\tau$$ $$=-\iiint\epsilon_0 \nabla \cdot (V_{2}\vec{E}_1)\,d\tau+\iiint \epsilon_0V_{2} \nabla \cdot \vec{E}_{1}\,d\tau.$$
Invoking Stokes' Theorem on the first integral gives $$-\iint\epsilon_0 (V_{2}\vec{E}_1) \cdot d\vec{S}+ \iiint\epsilon_0 V_{2} \nabla \cdot \vec{E}_1\,d\tau.$$ With the integration volume as all of space, the surface integral evaluates to zero for localized sources, since at large distances $(r\rightarrow\infty)$, the integrand falls off as $V_{2}\vec{E}_{1}\sim r^{-3}$. This leaves $$\iiint\epsilon_0V_{2} \nabla \cdot \vec{E}_1\,d\tau.$$ Invoking Gauss' Law, $$\iiint\epsilon_0 V_{2} \frac{\rho_1}{\epsilon_0}\,d\tau$$ $$=\iiint V_{2} \rho_1\,d\tau$$
In this form, it is clear to see why it represents potential energy, as you are building up charge distribution 1 in the presence of potential 2.
In this form you can also understand why it converges, as the $V\rho$ is always a finite number, because $V$ is never evaluated at any point that causes it to be infinity.
