Finite nuclei volume effect on ionization energy I am trying to understand following paper: Ionization Energy Difference
between Isotopes and its Effect on Isotope Abundance Measurement by Surface
Ionization Method.
I have question about section 3. Finite nucleus volume effect. I understand the main idea, but i dont understand the transformation of integral on quoted text:

Though this integral must be carried out over the all space, in fact,
  $\left(\varphi-\frac{Ze}{r}\right)$ is not zero only within the nucleus.
  On the other hand, the wave function of the s-state electron approaches to
  a finite value as $r\to 0$, and amounts to it already outside the nucleus.
  Consequently we can take $\psi^2$ out of the integrand and replace it by
  $\psi^2(0)$ calculated for Coulomb field of a point charge $Ze$.
  $$\begin{align}
\Delta E&=-e\psi^2(0)\int\left(\varphi(r) -\frac{Ze}{r}\right)\ dV \\
 &=-\frac{1}{6}e\psi^2(0)\int\left(\varphi-\frac{Ze}{r}\right)\Delta^2 r\ dV\\
 &=-\frac{1}{6}e\psi^2(0)\int r^2\Delta\left(\varphi-\frac{Ze}{r}\right)\ dV\ ,
\end{align}$$
  where we use the facts in the transformations of the integral that
  $\Delta^2 r= 6$, and the integrand is zero at infinitely remote surface.
  Now for all $r$ $\Delta\left(\frac{1}{r}\right)=-4\pi\delta(r)$ and
  $r^2\delta(r)=0$. According to Poisson's equation for electrostatic
  field, $\Delta\phi=-4\pi\rho$, where $\rho$ is the density of the
  charge distribution in the nucleus.

 A: I try to provide the solution of this problem. To be honest, an answer for energy shift, $\delta E$, depends on the model. For simplicity, one can consider that the nucleus is the ball of radius $R$. In this case we can write down Maxwell equation,
$$\partial\cdot{\bf E}=4\pi\rho,\quad \rho=3e/(4\pi R^3)\rightarrow E=\frac{er}{R^3}.$$
Then, for potential it easy to find
$$U(r)=\frac{e^2r^2}{2R^3}-\frac{3e^2}{2R},\quad r\leq R.$$
Note that the parameter of perturbative expansion is $R/a\ll 1$. Indeed, $R\sim 10^{-15}$ m and $a\sim 10^{-10}$ m. For area $r\geq R$, the potential is simply $U(r)=-e^2/r$. Perturbatio in our model is
$$V(r)=U(r)+\frac{e^2}{r}.$$
Now let us consider $1s$ state of hydrogen atom (as the simplest example). Its wavefunction is
$$\psi(r)=\frac{1}{\sqrt{\pi a^3}}e^{-r/a}.$$
Then, first order correction to ground state energy (=ionization energy) is
$$\delta E=\langle 100|V(r)|100\rangle$$
Performing integration (note that integration limit is $[0,R]$ because perturbation is non-zero only in this domain), you should obtain
$$\delta E=\frac{e^2}{2R^3}\left(3a^2-3R^2+\frac{2R^3}{a}-3e^{-2R/a}(a+R)^2\right).$$
Keeping in mind that $R/a\ll 1$ (expanding exp), you can find
$$\delta E= \frac{2e^2R^2}{5a^3}+\mathcal{O}(R^3).$$
Denoting $E_0=-e^2/(2a)$ (unperturbed ground state energy), you can see that
$$\delta E=-\frac{4}{5}\left(\frac{R}{a}\right)^2E_0\sim 10^{-10}.$$
In general, it is possible to consider different values $n$ & $l$. The energy shift for fixed $n$ & $l$ is
$$\delta E=\frac{e^2}{a}\int_{0}^{R}dr\,r^2\mathcal{R}_{nl}(r)\left[-\frac{3\beta}{2}+\frac{1}{a^2\beta}\left(\frac{r}{\beta}\right)^2\right], $$
where $\beta=R/a$ and $\mathcal{R}_{nr}$ is radial part of wave-function. For instance, in case of $n=2$, $l=1$, one can see that $\delta E_{2,1}\sim \beta^4$. In case of $n=2$ & $l=0$, it is relatively easy to show that $\delta E_{2,0}\sim \beta^2$. I hope that this will help.
