Solution (source code)

= Solution
{c}

A nonrelativistic $2p$ <wavefunction> vanishes linearly at the origin. Its probability density is therefore proportional to $r^2$, so its first-order <finite nuclear size correction> shift scales as $R^4$, rather than the $R^2$ result for $1s$. It is positive but much smaller.

More explicitly, the normalized <hydrogen-like ion> radial function obeys $R_{21}^2(r)\simeq Z^5r^2/(24a_0^5)$ near the origin. Since the <spherical harmonic> has unit angular norm,
$$
\Delta E_{2p}\simeq\frac{k_ZZ^5}{24a_0^5}\int_0^R\left(r^3-\frac{3r^4}{2R}+\frac{r^6}{2R^3}\right)dr
=\frac{k_ZZ^5R^4}{1120a_0^5}.
$$
Thus
$$
\boxed{\frac{\Delta E_{2p}}{\Delta E_{1s}}\simeq\frac1{448}\left(\frac{ZR}{a_0}\right)^2\simeq2.71\times10^{-7}.}
$$
This is about $2.4\,\mu\mathrm{eV}$ in the nonrelativistic model. At large $Z$, the relativistic $2p_{1/2}$ <Dirac spinor> has an $s$-type small component and a much stronger nuclear sensitivity than this nonrelativistic argument predicts; $2p_{3/2}$ remains much more suppressed. Neglecting <spin-orbit interaction> corrections would miss that distinction.