George Green's derivation of Poisson's equation I was reading George Green's An Essay on the Application of Mathematical Analysis to the Theories of Electricity and Magnetism, and I got confused on one step in his derivation of Poisson's Equation. Specifically, how does Green obtain conclude that:
$$\delta\left(2\pi a^2\varrho-\frac{2}{3}\pi b^2\varrho\right)=-4\pi\varrho.$$
Here are two pages to provide context; I understand everything except for the equality above.


 A: Let's first derive the value of $V$ inside the small sphere:
$$
V_\text{sphe} = \rho\int\frac{\text{d}x'\text{d}y'\text{d}z'}{r'},
$$
Where the sphere is sufficiently small such that $\rho$ can be considered constant. We can orientate the axes such that $p$ lies on the $z'$ axis. In spherical coordinates, the integral then has the form
$$
\begin{align}
V_\text{sphe} &= \rho\int_0^a\text{d}r'\int_0^{2\pi}\text{d}\varphi
\int_0^{\pi}\frac{r'^2\sin\theta}{\sqrt{r'^2 + b^2 - 2r'b\cos\theta}}\text{d}\theta\\
&=\frac{2\pi}{b}\rho\int_0^a r' \left(\sqrt{r'^2+b^2+2r'b}-\sqrt{r'^2+b^2-2r'b}\right)\text{d}r'\\
&= \frac{2\pi}{b}\rho\int_0^a r'(r'+b - |r'-b|)\text{d}r'\\
&=\frac{4\pi}{b}\rho\left[\int_0^b r'^2\text{d}r' + \int_b^a r'b\,\text{d}r'\right]\\
&= \frac{4\pi}{3}\rho b^2 + 2\pi\rho a^2 - 2\pi\rho b^2 = 2\pi\rho a^2 - \frac{2\pi}{3}\rho b^2.
\end{align}
$$
Since
$$
b^2 = (x-x_l)^2 + (y-y_l)^2 + (z-z_l)^2,
$$
we get
$$
\begin{align}
\frac{\partial b^2}{\partial x} &= 2(x-x_l),\qquad \frac{\partial^2 b^2}{\partial x^2} = 2 = \frac{\partial^2 b^2}{\partial y^2} = \frac{\partial^2 b^2}{\partial z^2}
\end{align}
$$
so that
$$
\delta b^2 = \frac{\partial^2 b^2}{\partial x^2} + \frac{\partial^2 b^2}{\partial y^2} + \frac{\partial^2 b^2}{\partial z^2} = 6
$$
and $\delta a^2 = 0$ since $a$ is a constant. Therefore,
$$
\delta V = \delta V_\text{sphe} = -\frac{2\pi}{3}\rho (\delta b^2) + \left(2\pi a^2 - \frac{2\pi}{3} b^2\right)\delta\rho = -4\pi\rho.
$$
The term with $\delta\rho$ disappears since $a$ and $b$ are exceedingly small.
A: Pulsar has already given a correct answer. In this answer we will use a slightly different (but equivalent) method, and we will keep a watchful eye on possible distributional contributions.

*

*George Green uses Gaussian units where Coulomb's constant $k_e=1$. He is considering a ball with radius $a$ and uniform charge density $\rho_0$, i.e. of total charge $$
Q~=~\frac{4\pi a^3}{3}\rho_0, \qquad \rho~=~\rho_0\theta(a-r) . \tag{1}$$
Here $\theta(a-r)$ is the Heaviside step function.


*Now since Coulomb's law has the same $1/r^2$ dependence as Newtonian gravity, we know from Newton's shell theorem that the electric field $\vec{E}$ is radial, and given by
$$ E_r~=~\frac{Q}{r^2}\theta(r-a)+\frac{Qr}{a^3} \theta(a-r) , \tag{2}$$
cf. e.g. this Phys.SE post and links therein.


*If we wish, we can integrate to find the potential
$$\begin{align} \Phi(r)~=~&\int_r^{\infty}\! dr^{\prime}~E_r(r^{\prime})\cr ~\stackrel{(2)}{=}~&
\frac{Q}{r}\theta(r-a)+\frac{Q}{2a}\left(3-\frac{r^2}{a^2}\right) \theta(a-r) \cr
~\stackrel{(1)}{=}~&2\pi \rho_0 \left(\frac{2a^3}{r}\theta(r-a)+\left(a^2-\frac{r^2}{3}\right) \theta(a-r)\right).\end{align}\tag{3}$$


*Instead of integrating and then applying the Laplacian, we can also differentiate the electric field (2) directly to deduce that
$$\begin{align} - \Delta\Phi~=~& \vec{\nabla}\cdot \vec{E}\cr
~=~&\frac{1}{r^2}\frac{\partial}{\partial r}  (r^2E_r)\cr
~\stackrel{(2)}{=}~&\frac{3Q}{a^3}\theta(a-r)\cr
~\stackrel{(1)}{=}~&4\pi\rho_0\theta(a-r)\cr~\stackrel{(1)}{=}~&4\pi\rho.\end{align}\tag{4}  $$
Eq. (4) is Poisson's equation in Gaussian's units. It can be viewed as a version of George Green's formula listed in the beginning of OP's post. The rhs. of eq. (4) unsurprisingly confirms that the charge density is $\rho_0$ inside the ball, and zero outside.


*Finally, note that when differentiating the Heaviside step functions from eq. (2), that the possible Dirac delta distribution terms $\delta(r-a)$ in eq. (4) precisely cancel, i.e. the ball has no surface charge.
