There are three variables involved, and it's important not to mix them up, or you'll go astray. There's the distance from the origin to the field point, call this $r$. There's the distance from the point on the surface of the disc being integrated to the field point, call this $\mathscr{R}$. And there's the distance from the origin on the disc to the point being integrated, call this $r'$. Assuming $\sigma$ is not a function of $r'$ the last equation will then look like:
$\frac{ \sigma}{4\pi \epsilon_o}\frac{1}{r} \sum_{l = 0} ^{\infty} p_l(\cos \phi) \left( \frac{r'}{r} \right)^l dt$,
for a general surface or volume element $dt$. You are integrating with respect to $r'$, so the $r$ comes outside the integral and you get (in polar coordinates):
$\frac{ \sigma}{4\pi \epsilon_o} \sum_{l = 0} ^{\infty} \frac{1}{r^{1+l}}\int p_l(\cos \phi) \left( r' \right)^l r'dr'd\phi.$
It's then just a matter of "pulling out" as many terms as you like, like:
$\frac{ \sigma}{4\pi \epsilon_or}\int r'dr' + \frac{ \sigma}{4\pi \epsilon_o r^2}\int r'^2\cos(\phi)dr'd\phi...$
to get an approximation for the potential to any accuracy you desire.
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.
Best Answer