that given in Ref. [@huang] in a different form.

Throughout the paper we take $V_{\rm Ueh}(r)=V_{\rm Ueh}^{\rm
  step}(r)$ and $r_n=\sqrt{5/3}\, r_{\rm rms}$, where the nuclear root-mean-square radii $r_{\rm rms}$ are taken from Ref. [@rms]. We have found that there is agreement to all digits presented for binding-energy shifts calculated using the step-function density expression for the Uehling potential, Eq. (\[eq:step\]), compared to the full double integral Eq. (\[eq:double\]) with a two-parameter Fermi distribution used for the nuclear density $\rho_{\rm nuc}(r)$, provided the same root-mean-square nuclear radius is used.

In this work the integration of Eq. (\[eq:step\]) is performed using the GNU Scientific Library adaptive integration routine QAGI [@gsl].

First-order shifts {#sec:first}
==================

The first-order Uehling correction to the binding energies of the valence electron is given by $$\delta \epsilon_i^{(1)}=-\langle \varphi_i |V_{\rm Ueh}|\varphi_i\rangle \ ,$$ where $\varphi_i$ is the valence electron wave function of state $i$.

The zeroth-order valence energies $\epsilon_i$ and wave functions $\varphi_i$ are found by solving the relativistic equations $$\label{eq:rel}
( c\boldsymbol{\alpha}\cdot {\bf p} +(\beta -1)c^2 -V_{\rm nuc}-V_{\rm
  el}) \varphi_i= \epsilon_i \varphi_i \ ,$$ where $V_{\rm nuc}$ and $V_{\rm el}$ are the nuclear and electronic potentials, $\boldsymbol{\alpha}$ and $\beta$ are Dirac matrices, and ${\bf p}$ is the momentum 