I'm reading Sect1.7 of Jackson's classical electrodynamics but I have trouble following his argument. Could someone help explain how exactly the Laplacian is evaluated in 1.30? Is it calculated with respect to the spherical coordinate? If so, where did the $\theta, \phi$ terms go? And why "it has a volume integral equal to 4$\pi$ for arbitrary a"? Also, I'd really appreciate it if someone could help me walk through the taylor expansion and the integration.
FWIW, $$\lim_{\varepsilon\to0^+}\frac{3}{4\pi}\frac{\varepsilon}{(r^2+\varepsilon)^{5/2}}~=~\delta^3(\vec{r})$$ is a well-known representation of the 3D Dirac delta distribution. It is straightforward to derive using test functions, cf. e.g. my related Phys.SE answer here.
The Laplacian is 1.30 is calculated in spherical coordinates. In the absence of $\theta, \phi$ dependence, the Laplacian is just $ \frac{1}{r} \frac{\partial^{2}}{\partial r^{2}}r$ and this does evaluate to the result given.
The volume integral is $4\pi$ because there is no $\theta, \phi$ dependence in the integrand and hence $\int_{0}^{2\pi} \int_{0}^{\pi}\sin(\theta)\, d\theta\, d\phi = 4\pi$.
The expression for $\nabla^2 f$ in spherical polar coordinates is
$$ \begin{equation} \frac{1}{r^{2}} \frac{\partial}{\partial r}\left(r^{2} \frac{\partial}{\partial r}f\right)+\frac{1}{r^{2} \sin ^{2} \phi} \frac{\partial^{2}}{\partial \theta^{2}}f+\frac{1}{r^{2} \sin \phi} \frac{\partial}{\partial \phi}\left(\sin \phi \frac{\partial}{\partial \phi}f\right) \end{equation}$$
If $f(r)$ depends only on the radius, as is true for the bracketed term in (1.30), then this reduces to
$$ \frac{1}{r^{2}} \frac{\partial}{\partial r}\left(r^{2} \frac{\partial}{\partial r}f(r)\right) $$
since partial derivatives of $f$ with respect to variables upon which $f$ does not depend vanish. In this case
$$\frac{1}{r^{2}} \frac{\partial}{\partial r}\left(r^{2} \frac{\partial}{\partial r}\frac{1}{\sqrt{r^2 + a^2}}\right) = - \frac{1}{r^{2}} \frac{\partial}{\partial r}\left(\frac{r^3}{(r^2 + a^2)^{3/2}}\right) = - \frac{3 a^2}{(a^2 + r^2)^{5/2}}$$
which is what Jackson has.
Next consider the volume integral of the bracketed expression, with the volume measure again in spherical polars:
$$\int_0^{2\pi} d\phi \int_{-1}^1 d\cos(\theta) \int_0^\infty dr\, r^2 \frac{3 a^2}{(a^2 + r^2)^{5/2}} $$
The integrand has no $\phi$ or $\theta$ dependence so those integrals may be trivially carried out to give $4\pi$, leaving
$$ 4\pi \int_0^\infty dr\, \frac{3 a^2 r^2}{(a^2 + r^2)^{5/2}} = 4\pi \int_0^\infty d\tau \, \frac{3 \sinh(\tau)^2}{\cosh(\tau)^4}$$
where the $a$ dependence has been eliminated by an $r = a \sinh(\tau)$ substitution. From here one can proceed by rewriting the integral as
$$\DeclareMathOperator{\sech}{sech} 4\pi \int_0^\infty d\tau\, 3 \tanh(\tau)^2\sech(\tau)^2 = 4\pi \int_0^1 d\tanh(\tau) \, 3 \tanh(\tau)^2 = 4\pi \left.\frac{3 \tanh(\tau)^3}{3} \right\rvert_0^1 = 4\pi$$
The Taylor expansion appears to have already been asked about on this site, with an answer available at https://physics.stackexchange.com/a/58107/250824. The general idea is that you should carry out a normal multivariate Taylor expansion and then adapt it to polar coordinates. You will find that the linear term vanishes on integration, leaving the constant and quadratic terms that Jackson has.
