How do one rigorously prove that the electric potential energy of an conducting sphere with charge $Q$ is $\frac{Q^2}{8\pi\epsilon_0R}$
We will show that one need not calculate the energy stored in the electrostatic field by integration. To that end, we proceed.
We start with the electrostatic energy $W$ stored as given by
$$W=\frac{\epsilon_0}{2}\int_V \vec E\cdot \vec E \,dV \tag 1$$
where $V$ is all of space.
We now assume that the field is induced by a charge distribution residing on a conductor of volume $V$ bounded by surface $S$. Inasmuch as the electric field is zero inside the conductor, the volume of integration in $(1)$ can be taken as all of space outside of the conductor.
We can change the form of $(1)$ using $\vec E=-\nabla \Phi$ along with the vector identity
$$\nabla \cdot \Phi\,\vec E = \Phi \nabla \cdot \vec E+\nabla \Phi \cdot \vec E\tag2$$
Then, we have for any conducting surface $S$
$$\begin{align} W&=\frac{\epsilon_0}{2}\int_V \vec E\cdot \vec E \,dV \\\\ &=-\frac{\epsilon_0}{2}\int_V \nabla \Phi \cdot \vec E \,dV \tag 3\\\\ &=\frac{\epsilon_0}{2}\int_V \Phi \nabla \cdot \vec E-\nabla \cdot (\Phi \vec E) \,dV \tag 4\\\\ &=\frac{\epsilon_0}{2}\oint_S\Phi (\hat n \cdot \vec E)dS \tag 5\\\\ &=\frac{\epsilon_0}{2}\left. \Phi(\vec r)\right|_{\vec r\in S}\frac{Q}{\epsilon_0} \tag 6\\\\ &=\frac12Q\left. \Phi(\vec r)\right|_{\vec r\in S} \tag 7 \end{align}$$
NOTES:
In arriving at $(3)$, we used $\vec E = -\nabla \Phi$.
In going from $(3)$ to $(4)$, we used the vector identity $(2)$.
In going from $(4)$ to $(5)$, we used the Divergence Theorem, recognizing that the unit normal points out of $V$ and therefore into the surface of the conductor; we then absorbed the minus sign wherein the final result, the unit normal points out of the conductor's surface. We also used that Gauss's Law, $\nabla \cdot \vec E=\rho/\epsilon_0=0$ in $V$.
In going from $(5)$ to $(6)$, we used the fact that the potential is constant on the conducting surface and the integral form of Gauss's Law, $\oint \vec E\cdot \hat n dS=Q/\epsilon_0$.
Thus, the electrostatic energy stored is given by
$$\bbox[5px,border:2px solid #C0A000]{W=\frac12 Q\left. \Phi(\vec r)\right|_{\vec r\in S}} \tag 8$$
For a conducting sphere, the potential on the surface is $\Phi = \frac{Q}{4\pi \epsilon_0R}$ which from $(8)$, immediately gives the result
$$W=\frac{Q^2}{8\pi\epsilon_0R}$$
as expected!
So, we see that one does not need to integrate $\vec E\cdot \vec E$ to determine the stored energy! Rather, we can simply use $(8)$ directly.