There are two main questions here: a technical question about the "exact" eigenvalues of the Hamiltonian
$$
H=-\frac{1}{2m}\sum_{i=1}^D\frac{\partial^2}{\partial x_i^2}+U(\vec{r})+\lambda\delta(\vec{r})\equiv H^{(0)}+\lambda\delta(\vec{r}) \ ,
$$
and a conceptual one about incorporating higher relativistic corrections beyond first-order perturbation theory. The TL;DR version of the answer is: the eigenvalue problem of this $H$ is ill-defined and there are no well-behaving bound state eigenvalues (with a possible loophole mentioned in Addendum 4), but this is not a problem from a physical point of view because solving the Schrodinger equation with a Darwin-like effective potential would not be suitable to account for higher relativistic effects anyway. Let us discuss both issues.
Eigenvalue problem with a repulsive Dirac delta term
Consider the above Hamiltonian in $D$ spatial dimensions with an arbitrary potential $U(\vec{r})$ ($U(\vec{r})=-Z\alpha/r$ for the Coulomb potential in $D=3$), and $\lambda>0$ (for most $U(\vec{r})$, $H$ would not be bounded from below with $\lambda<0$, $D\geq2$, as shown here).
To get a formal solution of $H\psi(\vec{r})=E\psi(\vec{r})$, we expand $\psi(\vec{r})$ in the basis of the eigenfunctions of $H^{(0)}$ (the usual Coulomb eigenfunctions for $U(\vec{r})=-Z\alpha/r$):
$$
\psi(\vec{r})=\sum_k a_k\phi_k(\vec{r}) \ .
$$
Inserting this into $H\psi(\vec{r})=E\psi(\vec{r})$ and using $H^{(0)}\phi_k(\vec{r})=\varepsilon_k\phi_k(\vec{r})$ as well as the properties of the delta function lead to
$$
\psi(\vec{r})=\sum_k\left(\lambda\psi(0)\frac{\phi_k^*(0)}{E-\varepsilon_k}\right)\phi_k(\vec{r})=\lambda\psi(0){\cal{G}}(\vec{r},0;E) \ ,
\tag{$\star$}
$$
where ${\cal{G}}(\vec{r},\vec{r}{}';E)$ is the Green's function of $H^{(0)}$:
$$
{\cal{G}}(\vec{r},\vec{r}{}';E)=\sum_k\frac{\phi_k(\vec{r})\phi^*_k(\vec{r}{}')}{E-\varepsilon_k}=\left\langle \vec{r}\left |\frac{1}{E-H^{(0)}}\right|\vec{r}{}'\right\rangle \ .
$$
Note that I assumed $E$ to be a bound state energy distinct from any $\varepsilon_k$, so no need for any $i0^+$ prescription. An alternative way to derive $(\star)$ is to start from the homogeneous Lippmann-Schwinger equation for bound states (see e.g. notes 38 of Littlejohn):
$$
\psi(\vec{r})=\int\mathrm{d}^Dr \, {\cal{G}}(\vec{r},\vec{r}{}';E)W(\vec{r}{}')\psi(\vec{r}{}') \ ,
$$
and just substitute $W(\vec{r})=\lambda\delta(\vec{r})$. Setting $\vec{r}=0$ in $(\star)$ gives an algebraic equation from which the exact eigenvalues could in principle be found, since everything there except for $E$ is known by construction:
$$
\frac{1}{\lambda}={\cal{G}}(0,0;E) \ .
\tag{$\star\star$}
$$
Once a specific $E_i$ is found, it can be substituted back into $(\star)$, and $\psi_i(0)$ can be found from the requirement $\langle\psi_i|\psi_i\rangle=1$; the eigenfunction belonging to that $E_i$ is then completely determined. To make connection with perturbation theory, we can rearrange $(\star\star)$ as
$$
\frac{1}{\lambda}=\frac{|\phi_i(0)|^2}{E_i-\varepsilon_i}+\sum_{k\neq i}\frac{|\phi_k(0)|^2}{E_i-\varepsilon_k} \ ,
$$
$$
E_i=\varepsilon_i+\lambda|\phi_i(0)|^2+\lambda(E_i-\varepsilon_i)\sum_{k\neq i}\frac{|\phi_k(0)|^2}{E_i-\varepsilon_k}
$$
for some $\varepsilon_i$ that is assumed to be close to $E_i$, and for small $\lambda$ (and assuming no degeneracies), we can start iterating this equation to get
$$
E_i=\varepsilon_i+\lambda|\phi_i(0)|^2+\lambda^2|\phi_i(0)|^2\sum_{k\neq i}\frac{|\phi_k(0)|^2}{\varepsilon_i-\varepsilon_k}+{\cal{O}}(\lambda^3) \ .
$$
These are just the standard Rayleigh-Schrodinger perturbation corrections.
In Am. J. Phys. 43 301 (1975), Atkinson and Crater used the above formalism to find exact eigenvalues for various $U(\vec{r})$ in $D=1$; however, they also showed that the question becomes ill-defined for $D\geq2$, as the Green's function becomes divergent in the $\vec{r}\rightarrow0$ limit. A pedestrian way to see this for the Coulomb potential is to write
$$
{\cal{G}}(0,0;E)=\sum_{n=1}^{\infty}\frac{|\phi_{n00}(0)|^2}{E-\varepsilon_n}+\int_0^\infty\mathrm{d}\varepsilon\frac{|\phi_{\varepsilon 00}(0)|^2}{E-\varepsilon}
$$
for the bound and scattering part of the hydrogenic spectrum (only $S$ states contribute, so $l=m_l=0$), and use
$$
|\phi_{n00}(0)|^2=\frac{m(Z\alpha)^3}{\pi n^3} \ ,
$$
$$
\begin{aligned}
|\phi_{\varepsilon 00}(0)|^2
&=\frac{pm}{2\pi^2}\left|\Gamma\left(1+i\frac{m Z\alpha}{p}\right)\right|^2\exp\left(\frac{\pi m Z\alpha}{p}\right) \\
&=
\frac{m^2Z\alpha}{\pi}
\left[
1-\exp\left(-2\pi\frac{mZ\alpha}{p}\right)
\right]^{-1}
\ ,
\end{aligned}
$$
where $p=\sqrt{2m \varepsilon}$ and $\Gamma(z)$ is the Euler gamma function. This shows $|\phi_{n 00}(0)|^2\sim n^{-3}$ and $|\phi_{\varepsilon 00}(0)|^2\sim\sqrt{\varepsilon}$ for large $\varepsilon$. The sum over bound states converges, but the integral over scattering states is divergent, and no meaningful eigenvalues can be extracted from $(\star\star)$.
Alternatively, one can show that ${\cal{G}}(\vec{r},0,E)\sim1/r$ for small $r$ (see e.g. Sec. 4.3.1-4.3.2 of Jentschura & Adkins), which by $(\star)$ means $\psi(\vec{r})\sim1/r$ for small $r$. But then $\delta(\vec{r})\psi(\vec{r})$ is ill-defined in the Schrodinger equation, and we ran into a contradiction.
Atkinson and Crater tried to regularize the problem by treating the delta function term as the limit of a constant potential barrier, and found that the effect of the potential well vanishes when the delta function limit is performed. These results, however, contradict the ones obtained by first-order perturbation theory, and highlight that the delta function term is extremely sensitive to the method of regularization (see also Addendum 2).
Higher relativistic corrections
The fact that the eigenvalue problem of $H$ is not well-defined is at first rather uncomfortable: there are no well-defined exact eigenvalues, and all perturbative energy corrections beyond first-order are divergent (this follows from the large $\varepsilon$ asymptotics of $|\phi_{\varepsilon 00}(0)|^2$, just like in the exact case; see also Schwartz, Ann. Phys. 6 2 156 (1959), where the first-order wave function is analytically constructed for the delta potential, making the divergence of the second-order energy explicit). But at the same time, the first-order correction is finite, and gives a correct and important contribution to the fine structure of hydrogen.
The punchline is that spin-orbit and Darwin-like interactions cannot be used universally, because they are not "fundamental" as e.g. a Coulomb interaction is; they are only effective potentials specifically constructed to imitate the results of a more complete theory (in this case, the Dirac theory) up to some perturbative order. Higher-order effects are incorporated by constructing higher-order effective operators and using those in perturbation theory as well, not by trying to insert the same low-order approximation in some high-order perturbative (or all-order "exact") correction. To make this clearer, consider the ground state energy of the Dirac Hamiltonian:
$$
\begin{aligned}
{\cal{E}}
&=
m\sqrt{1-(Z\alpha)^2} \\
&=
m-\frac{1}{2}m(Z\alpha)^2-\frac{1}{8}m(Z\alpha)^4-\frac{1}{16}m(Z\alpha)^6+{\cal{O}}(m(Z\alpha)^8) \ .
\end{aligned}
$$
If instead of the exact energy, we just wanted to find the energy to say $(Z\alpha)^6$ accuracy, then we should construct the appropriate effective operators to that order in the Foldy-Wouthuysen framework:
$$
{\cal{H}}_\text{eff}=m+{\cal{H}}^{(2)}+{\cal{H}}^{(4)}+{\cal{H}}^{(6)} \ ,
$$
where
$$
{\cal{H}}^{(2)}=\frac{1}{2m}\vec{p}^2+U \ ,
$$
$$
{\cal{H}}^{(4)}=-\frac{1}{8m^3}\vec{p}{}^4+\frac{1}{8m^2}\nabla^2U+\frac{1}{2m^2}\vec{S}\cdot(\nabla U\times\vec{p}) \ ,
$$
$$
{\cal{H}}^{(6)}=\frac{1}{16m^5}\vec{p}{}^6-\frac{3}{64m^4}\{\vec{p}{}^2,\nabla^2U+4\vec{S}\cdot(\nabla U\times\vec{p})\}+\frac{5}{128m^4}[\vec{p}{}^2,[\vec{p}{}^2,U]]+\frac{(\nabla U)^2}{8m^3} \ .
$$
Here, ${\cal{H}}^{(2)}$ is the usual nonrelativistic Hamiltonian, and ${\cal{H}}^{(4)}$ contains the mass-velocity, Darwin, and spin-orbit terms; ${\cal{H}}^{(6)}$ is much more complicated with no simple interpretation.
The $(Z\alpha)^2$ and $(Z\alpha)^4$ contributions to the energy are simply given by first-order perturbation theory with the non-relativistic ground state:
$$
E^{(2)}=\langle\phi_0|{\cal{H}}^{(2)}|\phi_0\rangle=-\frac{1}{2}m(Z\alpha)^2 \ ,
$$
$$
E^{(4)}=\langle\phi_0|{\cal{H}}^{(4)}|\phi_0\rangle=-\frac{1}{8}m(Z\alpha)^4 \ .
$$
However, to get the complete $(Z\alpha)^6$ contribution, a second-order correction must be calculated with ${\cal{H}}^{(4)}$ as well, since ${\cal{H}}^{(4)}\sim(Z\alpha)^4$ and $\varepsilon_0-\varepsilon_k\sim(Z\alpha)^2$:
$$
E^{(6)}=\langle\phi_0|{\cal{H}}^{(6)}|\phi_0\rangle+\sum_{k\neq0}\frac{\langle\phi_0|{\cal{H}}^{(4)}|\phi_k\rangle\langle\phi_k|{\cal{H}}^{(4)}|\phi_0\rangle}{\varepsilon_0-\varepsilon_k} \ .
$$
We saw on the example of the Darwin interaction that the second term is divergent; it turns out that the first term is divergent too, but when added together in an appropriate regularization (like $U_\lambda(\vec{r})= U(\vec{r})(1-\exp(-\lambda r))$ with $\lambda\rightarrow\infty$ in the end), their sum ends up finite and gives the correct $E^{(6)}=-m(Z\alpha)^6/16$ value (Pachucki, Phys. Rev. A 56 297 (1997)).
Similar divergences and cancellations are expected in higher $Z\alpha$ orders as well. When you tried to solve the Schrodinger equation for the Darwin term (without any higher-order operators) you used an effective interaction beyond its range of applicability, and no meaningful, consistent result can be expected from this.
For high $Z$ (like $Z=92$ in your example) the convergence of the $Z\alpha$ expansion is very slow, and there is no really better choice than to use the full solution of the Dirac equation. Of course, in that solution you cannot separate the effect of mass-velocity, Darwin and spin-orbit terms anymore, but this is natural, as such a separation is only possible in leading order: in higher orders, these effects combine in very complicated ways as already shown by ${\cal{H}}^{(6)}$.
(About the possible failure of the Dirac equation when $Z\alpha>1$: currently we do not even know what is beyond $Z=118$, so no one has any idea what to expect from $Z>137$. But note that solving the Dirac equation with not a point-like but a finite nucleus pushes the value of critical $Z$ quite higher.)
Addendum 1: analytic derivation of the first-order wave function
Here I derive the first-order wave function for the hydrogenic ground state in a Dirac delta perturbation to complement the result of Schwartz (Eq. (15) in his paper; he only gave the result without derivation). The Hamiltonian in atomic units reads
$$
H=-\frac{1}{2}\nabla^2-\frac{1}{r}+\lambda\delta(\vec{r})\equiv H^{(0)}+W \ ,
$$
and the zeroth-order wave function and zeroth- and first-order energy is
$$
\phi^{(0)}(\vec{r})=\frac{1}{\sqrt{\pi}}\exp(-r) \ \ , \ \ E^{(0)}=-\frac{1}{2} \ \ , \ \ E^{(1)}=\langle\phi^{(0)}|W|\phi^{(0)}\rangle=\frac{\lambda}{\pi} \ .
$$
If
$$
H^{(0)}\psi^{(1)}+W\psi^{(0)}=E^{(0)}\psi^{(1)}+E^{(1)}\psi^{(0)}
$$
is solved for the first-order wave function $\psi^{(1)}$, then the second-order energy can be readily computed as
$$
E^{(2)}=\langle\phi^{(0)}|W|\phi^{(1)}\rangle-E^{(1)}\langle\phi^{(0)}|\phi^{(1)}\rangle \ .
\tag{$\star\star\star$}
$$
Note that I am not using intermediate normalization, so in general $\langle\phi^{(0)}|\phi^{(1)}\rangle\neq0$.
The equation to solve is
$$
\left[-\frac{1}{2}\nabla^2-\frac{1}{r}+\frac{1}{2}\right]\phi^{(1)}(\vec{r})+\frac{1}{\sqrt{4\pi}}\frac{\lambda}{\pi}\left[2\pi\delta(\vec{r})-2\right]\exp(-r)=0 \ .
$$
Since $\phi^{(1)}$ must have $S$ symmetry, we may substitute the Ansatz
$$
\phi^{(1)}(\vec{r})=\frac{1}{\sqrt{4\pi}}\frac{\lambda}{\pi}\chi(r)\exp(-r) \ .
$$
After carrying out the differentiations, we get
$$
\nabla^2\chi(r)-2\chi'(r)=4\pi\delta(\vec{r})-4 \ ,
$$
where $'$ means differentation with respect to $r$, and actually
$$
\nabla^2\chi(r)=\chi''(r)+\frac{2}{r}\chi'(r) \ .
$$
Substituting $\chi(r)=\zeta(r)-1/r$ and using
$$
\nabla^2\frac{1}{r}=-4\pi\delta(\vec{r}) \ \ \Rightarrow \ \
\nabla^2\chi(r)=\zeta''(r)+\frac{2}{r}\zeta'(r)+4\pi\delta(\vec{r})
$$
leads to
$$
\zeta''(r)+\frac{2}{r}\zeta'(r)-2\zeta'(r)-\frac{2}{r^2}+4=0 \ .
$$
One more substitution $\Omega(r)=\zeta'(r)$ yields the first-order differential equation
$$
\Omega'(r)+\frac{2}{r}\Omega(r)-2\Omega(r)-\frac{2}{r^2}+4=0 \ ,
$$
which is solved by
$$
\Omega(r)=c'\frac{\exp(2r)}{r^2}+\frac{2}{r}+2 \ .
$$
We must set $c'=0$, otherwise the solution is not normalizable. Then
$$
\zeta(r)=2r+2\ln(r)+c \ ,
$$
and
$$
\begin{aligned}
\phi^{(1)}(\vec{r})
&=\frac{1}{\sqrt{4\pi}}\frac{\lambda}{\pi}\left[2r+2\ln(r)-\frac{1}{r}+c\right]\exp(-r) \\
&=\frac{1}{\sqrt{4\pi}}\frac{\lambda}{\pi}\left[2r+2\ln(r)-\frac{1}{r}\right]\exp(-r)+\tilde{c}\phi^{(0)}(\vec{r})
\ .
\end{aligned}
$$
There is an arbitrary integration constant in the solution that gives a component to $\phi^{(1)}$ which is proportional to $\phi^{(0)}$. You could fix this constant by e.g. requiring the intermediate normalization $\langle\phi^{(0)}|\phi^{(1)}\rangle=0$; however, by inspecting $(\star\star\star)$, it is easy to see that $E^{(2)}$ is independent of this constant, so I just set $c=0$:
$$
\phi^{(1)}(\vec{r})=\frac{1}{\sqrt{4\pi}}\frac{\lambda}{\pi}\left[2r+2\ln(r)-\frac{1}{r}\right]\exp(-r) \ .
$$
This corresponds to Eq. (15) of Schwartz. Now you can see that $\psi(\vec{r})\sim1/r$ for small $r$ as claimed earlier, and therefore $E^{(2)}$ is divergent when you try to compute it with $(\star\star\star)$.
Actually, the first-order wave function also determines the third-order energy as a consequence of the $2n+1$ theorem of Wigner. If we choose $\phi^{(1)}$ to satisfy the intermediate normalization (by setting $c=2\ln(2)+2\gamma-5$, where $\gamma$ is the Euler-Mascheroni constant), then
$$
E^{(3)}=\langle\phi^{(1)}|W|\phi^{(1)}\rangle-E^{(1)}\langle\phi^{(1)}|\phi^{(1)}\rangle \ ,
$$
and it is immediately seen, that this, too, is divergent.
Addendum 2: numerical results
As mentioned previously, Atkinson and Crater tried to make the "Coulomb+delta" eigenvalue problem well-defined by regularizing the delta function. They treated it as the limit of a constant potential barrier, and found its effect to be vanishing when the limit is taken, thereby recovering the hydrogenic ground state. Here, I numerically confirm this result in two different ways.
Approach 1: variational optimization
Since $H$ is bounded from below for $\lambda>0$, we have by the variational theorem
$$
{\cal{E}}_\text{trial}=\langle\phi_\text{trial}|H|\phi_\text{trial}\rangle\geq E
$$
for any trial function, where $E$ is the exact ground state energy. We can then try converging to $E$ by minimizing the expectation value. The trial function is parametrized as the ${\cal{N}}$-term linear combination
$$
\phi_\text{trial}(\vec{r})=\sum_{\mu=1}^{\cal{N}}c_\mu\chi_\mu(\vec{r}) \ ,
$$
where the $S$ symmetry basis functions are simple normalized Gaussians:
$$
\chi_\mu(\vec{r})=\left(\frac{2\alpha_\mu}{\pi}\right)^{3/4}\exp(-\alpha_\mu r^2)
$$
with the adjustable parameter $\alpha_\mu>0$. The Schrodinger equation then becomes a generalized matrix eigenvalue equation:
$$
\boldsymbol{Hc}=E\boldsymbol{Sc} \ ,
$$
where $H_{\mu\nu}=\langle\chi_\mu|H|\chi_\nu\rangle$, $S_{\mu\nu}=\langle\chi_\mu|\chi_\nu\rangle$ can be calculated analytically:
$$
S_{\mu\nu}=\left(\frac{4\alpha_\mu\alpha_\nu}{(\alpha_\mu+\alpha_\nu)^2}\right)^{3/4} \ ,
$$
$$
H_{\mu\nu}=\left[3\frac{\alpha_\mu\alpha_\nu}{\alpha_\mu+\alpha_\nu}-2\sqrt{\frac{\alpha_\mu+\alpha_\nu}{\pi}}+\lambda\left(\frac{\alpha_\mu+\alpha_\nu}{\pi}\right)^{3/2}\right]S_{\mu\nu} \ ,
$$
the three terms corresponding to the kinetic energy, Coulomb attraction and delta interaction.
Due to my job, I have access to a very flexible variational solver, where basically everything was already implemented; I only had to build the very simple matrix elements of the delta interaction. In the calculation, the values of $c_\mu$ are determined by repeatedly solving the matrix eigenvalue equation, while the Gaussian parameters $\alpha_\mu$ are refined in a stochastic optimization cycle. I performed this variational calculation for more and more basis functions, and the ground state energy is seen to converge to the pure hydrogenic limit.
\begin{array}{r|ccccc}
{\cal{N}} & E \, (\lambda=0.0) & E \, (\lambda=0.01) & E \, (\lambda=0.1) & E \, (\lambda=1.0) & E \, (\lambda=10.0) \\
\hline
1 & -0.424413182 & -0.423651782 & -0.417062450 & -0.368511994 & -0.226365895 \\
5 & -0.499809832 & -0.498997446 & -0.499277868 & -0.499277712 & -0.499277627 \\
10 & -0.499999161 & -0.499997094 & -0.499992511 & -0.499997327 & -0.499997229 \\
20 & -0.499999986 & -0.499999655 & -0.499999241 & -0.499999424 & -0.499999071
\end{array}
For $\lambda<0$, the calculation of course quickly fails, since $H$ is not bounded from below.
Just to be sure, I checked the method for another perturbing potential, $W=\lambda r$. For small $\lambda$, the numerical variational results were in very good agreement with the known analytical PT corrections for this linear potential:
$$
{\cal{E}}=-\frac{1}{2}+\frac{3}{2}\lambda-\frac{3}{2}\lambda^2+\frac{27}{4}\lambda^3+{\cal{O}(\lambda^4)} \ .
$$
For example, for $\lambda=0.01$, I found $E_\text{var}=-0.485144$ and $E_\text{PT}=-0.485143$.
Approach 2: Dirac delta as a limit of Gaussians
When writing $\psi(\vec{r})=\chi(r)/(\sqrt{4\pi}r)$, the radial Schrodinger equation reads
$$
\left[-\frac{1}{2}\frac{\partial^2}{\partial r^2}-\frac{1}{r}+\lambda\frac{\delta(r)}{4\pi r^2}\right]\chi(r)=E\chi(r) \ .
$$
I regularized the radial delta function as
$$
\delta_\xi(r)=\frac{1}{\sqrt{\pi}\xi}\exp(-r^2/\xi^2) \ ,
$$
and calculated the lowest eigenvalue with the NDEigensystem routine of Mathematica for increasingly smaller values of $\xi$. The table shows once again that the ground state energy converges to that of a simple hydrogen as $\xi\rightarrow0^+$.
\begin{array}{c|ccc}
\xi & E \, (\lambda=0.1) & E \, (\lambda=1.0) & E \, (\lambda=10.0) \\
\hline
10^{-1} & -0.488128 & -0.447406 & -0.397968 \\
10^{-2} & -0.493161 & -0.484291 & -0.476065 \\
10^{-3} & -0.498289 & -0.497295 & -0.496416 \\
10^{-4} & -0.499724 & -0.499631 & -0.499550 \\
10^{-5} & -0.499965 & -0.499957 & -0.499950
\end{array}
It seems the delta function is "lost" within the Coulomb singularity. Note that this is not a numerical artifact, if we placed the Dirac delta somewhere else (by $\delta(r)\rightarrow\delta(r-r_0)$), the ground state energy would actually change. For $r_0=1$ and $\lambda=10.0$, I found that the energy converged to $E\approx-0.265$.
Addendum 3: final remarks
So, we have an eigenvalue problem that produces severe divergences when attacked with Green's function techniques or perturbation theory, but immediately collapses to a trivial, "delta funcion-less" problem when treated variationally, or with a regularization; results from this latter, trivial problem, however, are in contradiction with the observation that the finite first-order energy correction gives correct values. Actually, just after doing the numerical work, I found another paper (Velenik, Zivkovic, de Jeu, Murrell, Mol. Phys. 18 5 693 (1970)) which explicitly shows that variational calculations with the "Coulomb+delta" potential just collapse to the hydrogenic ground state and that beyond-leading-order perturbation theory becomes useless (see also Blinder, Phys. Rev. A 18 853 (1978) and Blinder, Theor. Chim. Acta 53 159 (1979)).
This eigenvalue problem is ill-posed, and just shows that one should not use low-order effective operators beyond first-order PT. Delta function interactions appear in many areas of QM (fine- and hyperfine structure of atoms, the annihilation contribution to the hyperfine splitting of positronium, vacuum polarization and other effects in QED, etc.) as low-energy approximations to some more fundamental process. There is no reason to believe that these interactions could be consistently used beyond their range of applicability. The fact that the second-order energy correction with a delta function perturbation is divergent has been known at least since the 1950's. In the context of higher-order hyperfine effects, people tried to smooth out the delta potential, so that the nucleus is not treated anymore as a point-like magnetic dipole, while for the Darwin term (and many other similar terms in bound state QED), the divergences are cancelled by higher-order effective operators in the Foldy-Wouthuysen formalism. There is nothing more to gain from low-order approximate interactions.
Addendum 4: self-adjoint extensions
It just came to my attention that while the original eigenvalue problem of $H$ with the usual boundary conditions is ill-posed, there exists a modification of the problem via the self-adjoint extension of $H$, which is well defined, and could, in principle, provide well-behaving eigenvalues (see e.g. here and here). This is a very interesting mathematical topic on its own right, but from a physical point of view, I doubt the eigenvalues would carry useful information beyond their first-order contributions. The sources of trouble in the original, ill-defined problem are the missing finer physical effects. Also, I would be very surprised if the eigenvalues of the modified, well-defined problem turned out to be different from the "delta-less" eigenvalues found with the variational or regularization approach.
Addendum 5: addressing the comments
1. Normalizability of the wave functon with continuum components
There are quite a few questions/doubts in the comments about the role of continuum states in perturbation theory, the first being about the normalizability of a wave function that has both discrete and continuum state contributions in some basis.
Let $H^{(0)}$ have both discrete and continuum eigenstates ($\{|\phi^{(0)}_k\rangle\}$ and $\{|\phi^{(0)}(\alpha)\rangle\}$ for some continuous $\alpha$ over some domain $\Omega$). These states are orthonormal in the following sense:
$$
\langle\phi^{(0)}_k|\phi^{(0)}_l\rangle=\delta_{kl} \ \ , \ \ \langle\phi^{(0)}(\alpha)|\phi^{(0)}(\alpha)\rangle=\delta(\alpha-\alpha') \ \ , \ \
\langle\phi^{(0)}_k|\phi^{(0)}(\alpha)\rangle=0 \ ,
$$
and complete, meaning that an arbitrary state of the function space can be expanded in them as
$$
|\Psi\rangle=\sum_kc_k|\phi^{(0)}_k\rangle+\int_\Omega\mathrm{d}\alpha \, c(\alpha)|\phi^{(0)}(\alpha)\rangle \ ,
$$
with components $c_k=\langle\phi^{(0)}_k|\Psi\rangle$, $c(\alpha)=\langle\phi^{(0)}(\alpha)|\Psi\rangle$.
Nothing prevents this state from being normalizable; for example, $\langle\Psi|\Psi\rangle=1$ if
$$
\sum_k|c_k|^2+\int_\Omega\mathrm{d}\alpha \, |c(\alpha)|^2=1 \ ,
$$
as follows from the orthogonality relations. For example, $|\Psi\rangle$ could be a bound eigenstate of some other Hamiltonian $H=H^{(0)}+W$. In this basis, the identity operator has the resolution
$$
I=\sum_k|\phi^{(0)}_k\rangle\langle\phi^{(0)}_k|+\int_\Omega\mathrm{d}\alpha \, |\phi^{(0)}(\alpha)\rangle\langle\phi^{(0)}(\alpha)| \ ,
$$
since this maps every state onto itself: $I|\Psi\rangle=|\Psi\rangle$. Note that continuum states are normalized to delta function, but this is not a problem since they always appear behind an integral over $\alpha$.
The first-order wave function and the second-order energy corrections formally read
$$
|\phi^{(1)}_n\rangle=\frac{P_n^\perp}{E^{(0)}_n-H^{(0)}}W|\phi^{(0)}_n\rangle \ ,
$$
$$
E^{(2)}_n=\langle\phi^{(0)}_n|W|\phi^{(1)}_n\rangle=\left\langle\phi^{(0)}_n\left|W\frac{P_n^\perp}{E^{(0)}_n-H^{(0)}}W\right|\phi^{(0)}_n\right\rangle \ ,
$$
where $P_n^\perp=I-|\phi^{(0)}_n\rangle\langle\phi^{(0)}_n|$ is a projector onto the subspace orthogonal to $|\phi^{(0)}_n\rangle$ (we assumed no degeneracies), and we used intermediate normalization (so that $P_n^\perp|\phi^{(1)}_n\rangle=|\phi^{(1)}_n\rangle$). Substituting the above resolution of identity gives
$$
|\phi^{(1)}_n\rangle=\sum_{k\neq n}\frac{\langle \phi_k^{(0)}|W|\phi_n^{(0)}\rangle}{E^{(0)}_n-E^{(0)}_k}|\phi_k^{(0)}\rangle+\int_\Omega\mathrm{d}\alpha \, \frac{\langle \phi^{(0)}(\alpha)|W|\phi_n^{(0)}\rangle}{E^{(0)}_n-E^{(0)}(\alpha)}|\phi(\alpha)^{(0)}\rangle \ ,
$$
$$
E^{(2)}_n=\sum_{k\neq n}\frac{|\langle \phi_k^{(0)}|W|\phi_n^{(0)}\rangle|^2}{E^{(0)}_n-E^{(0)}_k}+\int_\Omega\mathrm{d}\alpha \, \frac{|\langle \phi^{(0)}(\alpha)|W|\phi_n^{(0)}\rangle|^2}{E^{(0)}_n-E^{(0)}(\alpha)} \ .
$$
These continuum contributions tend to be discussed in textbooks (see e.g. Landau & Lifshitz III or Chapters 31 and 33 of Schiff).
Additionally, think of how the basis of plane waves is continuous in the infinite box limit, but you can still construct localized, normalizable wave packets (e.g. Gaussians) with them.
It is easy to check that the analytical first-order wave function I wrote down in Addendum 1 for the delta perturbation has a non-zero overlap with the hydrogenic continuum states, therefore it has continuum components. Just pick some value for $\varepsilon$, and compute the integral $\langle\phi^{(0)}_{\varepsilon00}|\phi^{(1)}\rangle$ numerically.
2. Box normalization of scattering states and the continuum limit
You seem to claim that continuum contributions can be neglected based on the following argument: you first put the system in a box, so that scattering states (when normalized as $\langle\psi_\alpha|\psi_{\alpha'}\rangle=\delta_{\alpha\alpha'}$) scale as $|\psi_{\alpha}\rangle\sim1/\sqrt{L}$, and then $|\psi_\alpha\rangle\langle\psi_\alpha|\sim1/L$ vanishes in the limit $L\rightarrow\infty$. The flaw in this argument is that this factor of $1/L$ is absorbed in the summation over $\alpha$ as it turns into an integration in the continuum limit.
To see this, consider first plane waves in 1D. In case of a finite box, we have
$$
\phi_n(x)=\frac{1}{\sqrt{L}}\exp(ik_nx) \ \ , \ \ k_n=\frac{2\pi n}{L} \ ,
$$
with orthonormality and completeness relations
$$
\int_{-\frac{L}{2}}^{+\frac{L}{2}}\mathrm{d}x \, \phi^*_n(x)\phi_m(x)=\delta_{mn} \ \ , \ \
\sum_{n=-\infty}^{+\infty}\phi_n(x)\phi^*_n(x')=\delta(x-x') \ .
$$
In the infinite space limit, we should have
$$
\phi(x;k)=\frac{1}{\sqrt{2\pi}}\exp(ikx) \ ,
$$
with
$$
\int_{-\infty}^{+\infty}\mathrm{d}x \, \phi^*(x;k)\phi(x;k')=\delta(k-k') \ \ , \ \
\int_{-\infty}^{+\infty}\mathrm{d}k \, \phi(x;k)\phi^*(x';k)=\delta(x-x') \ .
$$
To actually get this second relation (and not identically zero!), we must absorb the $1/L$ in the summation over $n$ as we pass to the continuum limit:
$$
\sum_{n=-\infty}^{+\infty}\phi_n(x)\phi^*_n(x')=
\underbrace{\sum_{n=-\infty}^{+\infty}\frac{2\pi}{L}}_{\rightarrow\int_{-\infty}^{+\infty}\mathrm{d}k}\underbrace{\frac{1}{2\pi}\exp(ik_n(x-x'))}_{\rightarrow\phi(x;k)\phi^*(x';k)}
\rightarrow\delta(x-x') \ .
$$
This is exactly what happens in the PT formulae too. In a box, the "original" continuum states have the property
$$
\langle\phi^{(0)}_\alpha|\phi^{(0)}_{\alpha'}\rangle=\delta_{\alpha\alpha'}\frac{L}{\kappa} \ ,
$$
where $\kappa>0$ is some uninteresting numerical factor. Of course, in the finite box treatment, you can re-normalize these states as
$$
|\psi^{(0)}_\alpha\rangle=\sqrt{\frac{\kappa}{L}}|\phi^{(0)}_{\alpha}\rangle \ \ , \ \ \langle\psi^{(0)}_\alpha|\psi^{(0)}_{\alpha'}\rangle=\delta_{\alpha\alpha'} \ ,
$$
and then you can treat everything as a sum; e.g.
$$
I=\sum_k|\phi^{(0)}_k\rangle\langle\phi^{(0)}_k|+\sum_\alpha |\psi^{(0)}_\alpha\rangle\langle\psi^{(0)}_\alpha| \ .
$$
But if you want to take the $L\rightarrow\infty$ limit, then that factor of $\kappa/L$ will be a part of the integration measure, and nothing tends to zero:
$$
\sum_\alpha\frac{|\langle\psi^{(0)}_\alpha|W|\phi^{(0)}_n\rangle|^2}{E_n^{(0)}-E_\alpha^{(0)}}=
\underbrace{\sum_\alpha\frac{\kappa}{L}}_{\rightarrow\int_\Omega\mathrm{d}\alpha}
\underbrace{\frac{|\langle\phi^{(0)}_\alpha|W|\phi^{(0)}_n\rangle|^2}{E_n^{(0)}-E_\alpha^{(0)}}}_{\rightarrow\frac{|\langle\phi^{(0)}(\alpha)|W|\phi^{(0)}_n\rangle|^2}{E_n^{(0)}-E^{(0)}(\alpha)}}
\rightarrow
\int_\Omega\mathrm{d}\alpha \, \frac{|\langle \phi^{(0)}(\alpha)|W|\phi_n^{(0)}\rangle|^2}{E^{(0)}_n-E^{(0)}(\alpha)} \ .
$$
