D.81 Dirac fine struc­ture Hamil­ton­ian

This note de­rives the fine struc­ture Hamil­ton­ian of the hy­dro­gen atom. This Hamil­ton­ian fixes up the main rel­a­tivis­tic er­rors in the clas­si­cal so­lu­tion of chap­ter 4.3. The de­riva­tion is based on the rel­a­tivis­tic Dirac equa­tion from chap­ter 12.12 and uses non­triv­ial lin­ear al­ge­bra.

Ac­cord­ing to the Dirac equa­tion, the rel­a­tivis­tic Hamil­ton­ian and wave func­tion take the form

\begin{displaymath}
H_D = m_{\rm e}c^2 \left(\begin{array}{cr} 1&0\ 0&-1 \end{...
...n{array}{c} \vec\psi^{ p}\ \vec\psi^{ n} \end{array}\right)
\end{displaymath}

where $m_{\rm e}$ is the mass of the elec­tron when at rest, $c$ the speed of light, and the $\sigma_i$ are the 2 $\times$ 2 Pauli spin ma­tri­ces of chap­ter 12.10. Sim­i­larly the ones and ze­ros in the shown ma­tri­ces are 2 $\times$ 2 unit and zero ma­tri­ces. The wave func­tion is a four-di­men­sion­al vec­tor whose com­po­nents de­pend on spa­tial po­si­tion. It can be sub­di­vided into the two-di­men­sion­al vec­tors $\vec\psi^{ p}$ and $\vec\psi^{ n}$. The two com­po­nents of $\vec\psi^{ p}$ cor­re­spond to the spin up and spin down com­po­nents of the nor­mal clas­si­cal elec­tron wave func­tion; as noted in chap­ter 5.5.1, this can be thought of as a vec­tor if you want. The two com­po­nents of the other vec­tor $\vec\psi^{ n}$ are very small for the so­lu­tions of in­ter­est. These com­po­nents would be dom­i­nant for states that would have neg­a­tive rest mass. They are as­so­ci­ated with the anti-par­ti­cle of the elec­tron, the positron.

The Dirac equa­tion is solv­able in closed form, but that so­lu­tion is not some­thing you want to con­tem­plate if you can avoid it. And there is re­ally no need for it, since the Dirac equa­tion is not ex­act any­way. To the ac­cu­racy it has, it can eas­ily be solved us­ing per­tur­ba­tion the­ory in es­sen­tially the same way as in de­riva­tion {D.79}. In this case, the small pa­ra­me­ter is 1/$c$: if the speed of light is in­fi­nite, the non­rel­a­tivis­tic so­lu­tion is ex­act. And if you ball­park a typ­i­cal ve­loc­ity for the elec­tron in a hy­dro­gen atom, it is only about one per­cent or so of the speed of light.

So, fol­low­ing de­riva­tion {D.79}, take the Hamil­ton­ian apart into suc­ces­sive pow­ers of 1$\raisebox{.5pt}{$/$}$$c$ as $H_D$ $\vphantom0\raisebox{1.5pt}{$=$}$ $H_{D,0}+H_{D,1}+H_{D,2}$ with

\begin{displaymath}
H_{D,0} = \left(\begin{array}{cc} m_{\rm e}c^2&0\ 0&-m_{\r...
...H_{D,2} = \left(\begin{array}{cc} V&0\ 0&V \end{array}\right)
\end{displaymath}

and sim­i­larly for the wave func­tion vec­tor:

\begin{displaymath}
\vec\psi_D =
\left(\begin{array}{c}\vec\psi^{ p}_0\ \vec...
...vec\psi^{ p}_4\ \vec\psi^{ n}_4\end{array}\right) +
\ldots
\end{displaymath}

and the en­ergy:

\begin{displaymath}
E_D = E_{D,0} + E_{D,1}+ E_{D,2} + E_{D,3} + E_{D,4} + \ldots
\end{displaymath}

Sub­sti­tu­tion into the Hamil­ton­ian eigen­value prob­lem $H_D\vec\psi_D$ $\vphantom0\raisebox{1.5pt}{$=$}$ $E_D\vec\psi_D$ and then col­lect­ing equal pow­ers of 1$\raisebox{.5pt}{$/$}$$c$ to­gether pro­duces again a sys­tem of suc­ces­sive equa­tions, just like in de­riva­tion {D.79}:

\begin{eqnarray*}
c^2:&&
\left[
\left(\begin{array}{cc} m_{\rm e}c^2&0\ 0&-m...
..._0\ \vec\psi^{ n}_0\end{array}\right)
= 0
\hspace{.45truein}
\end{eqnarray*}


\begin{eqnarray*}
c^1:&&
\left[
\left(\begin{array}{cc} m_{\rm e}c^2&0\ 0&-m...
...n{array}{c}\vec\psi^{ p}_0\ \vec\psi^{ n}_0\end{array}\right)
\end{eqnarray*}


\begin{eqnarray*}
c^0:&&
\left[
\left(\begin{array}{cc} m_{\rm e}c^2&0\ 0&-m...
...n{array}{c}\vec\psi^{ p}_0\ \vec\psi^{ n}_0\end{array}\right)
\end{eqnarray*}


\begin{eqnarray*}
c^{-1}:&&
\left[
\left(\begin{array}{cc} m_{\rm e}c^2&0\ 0...
...n{array}{c}\vec\psi^{ p}_0\ \vec\psi^{ n}_0\end{array}\right)
\end{eqnarray*}


\begin{eqnarray*}
c^{-2}:&&
\left[
\left(\begin{array}{cc} m_{\rm e}c^2&0\ 0...
...n{array}{c}\vec\psi^{ p}_0\ \vec\psi^{ n}_0\end{array}\right)
\end{eqnarray*}


\begin{eqnarray*}
c^{-3}:&& \smash{\cdots}\hspace{4.13truein}
\end{eqnarray*}

The first, or­der $c^2$, eigen­value prob­lem has en­ergy eigen­val­ues $\pm{m}_ec^2$, in other words, plus or mi­nus the rest mass en­ergy of the elec­tron. The so­lu­tion of in­ter­est is the phys­i­cal one with a pos­i­tive rest mass, so the de­sired so­lu­tion is

\begin{displaymath}
E_{D,0}=m_{\rm e}c^2
\qquad
\vec\psi^{ p}_0 = \mbox{still arbitrary}
\qquad
\vec\psi^{ n}_0 = 0
\end{displaymath}

Plug that into the or­der $c^1$ equa­tion to give, for top and bot­tom sub­vec­tors

\begin{displaymath}
0 = E_{D,1} \vec\psi^{ p}_0
\qquad
- 2 m_{\rm e}c^2 \vec\psi^{ n}_1 = -\sum_i c{\widehat p}_i\sigma_i \vec\psi^{ p}_0
\end{displaymath}

It fol­lows from the first of those that the first or­der en­ergy change must be zero be­cause $\vec\psi^{ p}_0$ can­not be zero; oth­er­wise there would be noth­ing left. The sec­ond equa­tion gives the lead­ing or­der val­ues of the sec­ondary com­po­nents, so in to­tal

\begin{displaymath}
E_{D,1} = 0
\qquad
\vec\psi^{ p}_1 = \mbox{still arbitra...
...j \frac{1}{2m_{\rm e}c}{\widehat p}_j\sigma_j \vec\psi^{ p}_0
\end{displaymath}

where the sum­ma­tion in­dex $i$ was re­named to $j$ to avoid am­bi­gu­ity later.

Plug all that in the or­der $c^0$ equa­tion to give

\begin{displaymath}
0 = - \frac{1}{2m_{\rm e}}\sum_i\sum_j{\widehat p}_i{\wideh...
...j \frac{1}{2m_{\rm e}c}{\widehat p}_j\sigma_j \vec\psi^{ p}_1
\end{displaymath}

The first of these two equa­tions is the non­rel­a­tivis­tic Hamil­ton­ian eigen­value prob­lem of chap­ter 4.3. To see that, note that in the dou­ble sum the terms with $j$ $\raisebox{.2pt}{$\ne$}$ $i$ pair­wise can­cel since for the Pauli ma­tri­ces, $\sigma_i\sigma_j+\sigma_j\sigma_i$ $\vphantom0\raisebox{1.5pt}{$=$}$ 0 when $j$ $\raisebox{.2pt}{$\ne$}$ $i$. For the re­main­ing terms in which $j$ $\vphantom0\raisebox{1.5pt}{$=$}$ $i$, the rel­e­vant prop­erty of the Pauli ma­tri­ces is that $\sigma_i\sigma_i$ is one (or the 2 $\times$ 2 unit ma­trix, re­ally,) giv­ing

\begin{displaymath}
\frac{1}{2m_{\rm e}}\sum_i\sum_j{\widehat p}_i{\widehat p}_...
...
= \frac{1}{2m_{\rm e}}\sum_i{\widehat p}_i^2 + V
\equiv H_0
\end{displaymath}

where $H_0$ is the non­rel­a­tivis­tic hy­dro­gen Hamil­ton­ian of chap­ter 4.3.

So the first part of the or­der $c^0$ equa­tion takes the form

\begin{displaymath}
H_0\vec\psi^{ p}_0=E_{D,2}\vec\psi^{ p}_0
\end{displaymath}

The en­ergy $E_{D,2}$ will there­fore have to be a Bohr en­ergy level $E_n$ and each com­po­nent of $\vec\psi^{ p}_0$ will have to be a non­rel­a­tivis­tic en­ergy eigen­func­tion with that en­ergy:

\begin{displaymath}
E_{D,2} = E_n
\qquad
\vec\psi^{ p}_0 =
\sum_l\sum_m c_{...
...{nlm}{\uparrow}+ \sum_l\sum_m c_{lm{-}} \psi_{nlm}{\downarrow}
\end{displaymath}

The sum mul­ti­ply­ing ${\uparrow}$ is the first com­po­nent of vec­tor $\vec\psi^{ p}_0$ and the sum mul­ti­ply­ing ${\downarrow}$ the sec­ond. The non­rel­a­tivis­tic analy­sis in chap­ter 4.3 was in­deed cor­rect as long as the speed of light is so large com­pared to the rel­e­vant ve­loc­i­ties that 1$\raisebox{.5pt}{$/$}$$c$ can be ig­nored.

To find out the er­ror in it, the rel­a­tivis­tic ex­pan­sion must be taken to higher or­der. To or­der $c^{-1}$, you get for the top vec­tor

\begin{displaymath}
0 = - (H_0 - E_n) \vec\psi^{ p}_1 + E_{D,3} \vec\psi^{ p}_0
\end{displaymath}

Now if $\vec\psi^{ p}_1$ is writ­ten as a sum of the eigen­func­tions of $H_0$, in­clud­ing $\vec\psi^{ p}_0$, the first term will pro­duce zero times $\vec\psi^{ p}_0$ since $(H_0-E_n)\vec\psi^{ p}_0$ $\vphantom0\raisebox{1.5pt}{$=$}$ 0. That means that $E_{D,3}$ must be zero. The ex­pan­sion must be taken one step fur­ther to iden­tify the rel­a­tivis­tic en­ergy change. The bot­tom vec­tor gives

\begin{displaymath}
\vec\psi^{ n}_3 = \sum_j \frac{1}{2m_{\rm e}c}{\widehat p}...
...j \frac{1}{2m_{\rm e}c}{\widehat p}_j\sigma_j \vec\psi^{ p}_0
\end{displaymath}

To or­der $c^{-2}$, you get for the top vec­tor

\begin{displaymath}
0 = - (H_0 - E_n) \vec\psi^{ p}_2
- \sum_i\sum_j
{\wideh...
...ehat p}_j\sigma_j \vec\psi^{ p}_0
+ E_{D,4} \vec\psi^{ p}_0
\end{displaymath}

and that de­ter­mines the ap­prox­i­mate rel­a­tivis­tic en­ergy cor­rec­tion.

Now re­call from de­riva­tion {D.79} that if you do a non­rel­a­tivis­tic ex­pan­sion of an eigen­value prob­lem $(H_0+H_1)\psi$ $\vphantom0\raisebox{1.5pt}{$=$}$ $E\psi$, the equa­tions to solve are (D.55) and (D.56);

\begin{displaymath}
(H_0-E_{{\vec n},0})\psi_{{\vec n},0} = 0
\qquad
(H_0-E_{...
...})\psi_{{\vec n},1}
= - (H_1-E_{{\vec n},1})\psi_{{\vec n},0}
\end{displaymath}

The first equa­tion was sat­is­fied by the so­lu­tion for $\vec\psi^{ p}_0$ ob­tained above. How­ever, the sec­ond equa­tion presents a prob­lem. Com­par­i­son with the fi­nal Dirac re­sult sug­gests that the fine struc­ture Hamil­ton­ian cor­rec­tion $H_1$ should be iden­ti­fied as

\begin{displaymath}
H_1 \stackrel{\mbox{?}}{=}
\sum_i\sum_j {\widehat p}_i\sigma_i \frac{V-E_n}{4m_{\rm e}^2c^2}{\widehat p}_j\sigma_j
\end{displaymath}

but that is not right, since $E_n$ is not a phys­i­cal op­er­a­tor, but an en­ergy eigen­value for the se­lected eigen­func­tion. So map­ping the Dirac ex­pan­sion straight­for­wardly onto a clas­si­cal one has run into a snag.

It is maybe not that sur­pris­ing that a two-di­men­sion­al wave func­tion can­not cor­rectly rep­re­sent a truly four-di­men­sion­al one. But clearly, what­ever is se­lected for the fine struc­ture Hamil­ton­ian $H_1$ must at least get the en­ergy eigen­val­ues right. To see how this can be done, the op­er­a­tor ob­tained from the Dirac equa­tion will have to be sim­pli­fied. Now for any given $i$, the sum over $j$ in­cludes a term $j$ $\vphantom0\raisebox{1.5pt}{$=$}$ $i$, a term $j$ $\vphantom0\raisebox{1.5pt}{$=$}$ ${\overline{\imath}}$, where ${\overline{\imath}}$ is the num­ber fol­low­ing $i$ in the cyclic se­quence $\ldots123123\ldots$, and it in­volves a term $j$ $\vphantom0\raisebox{1.5pt}{$=$}$ ${\overline{\overline{\imath}}}$ where ${\overline{\overline{\imath}}}$ pre­cedes $i$ in the se­quence. So the Dirac op­er­a­tor falls apart into three pieces:

\begin{displaymath}
H_1 \stackrel{\mbox{?}}{=}
\sum_i {\widehat p}_i\sigma_i \...
...ne{\overline{\imath}}}}\sigma_{{\overline{\overline{\imath}}}}
\end{displaymath}

or us­ing the prop­er­ties of the Pauli ma­tri­ces that $\sigma_i\sigma_i$ $\vphantom0\raisebox{1.5pt}{$=$}$ 1, $\sigma_i\sigma_{{\overline{\imath}}}$ $\vphantom0\raisebox{1.5pt}{$=$}$ ${\rm i}\sigma_{{\overline{\overline{\imath}}}}$, and $\sigma_i\sigma_{{\overline{\overline{\imath}}}}$ $\vphantom0\raisebox{1.5pt}{$=$}$ $-{\rm i}\sigma_{{\overline{\overline{\imath}}}}$ for any $i$,
\begin{displaymath}
H_1 \stackrel{\mbox{?}}{=}
\sum_i {\widehat p}_i \frac{V-E...
...{{\overline{\overline{\imath}}}}\sigma_{{\overline{\imath}}} %
\end{displaymath} (D.59)

The ap­proach will now be to show first that the fi­nal two terms are the spin-or­bit in­ter­ac­tion in the fine struc­ture Hamil­ton­ian. Af­ter that, the much more tricky first term will be dis­cussed. Reno­tate the in­dices in the last two terms as fol­lows:

\begin{displaymath}
H_{1,\mbox{\scriptsize spin-orbit}} =
{\rm i}\sum_i {\wide...
...}{4m_{\rm e}^2c^2}{\widehat p}_{{\overline{\imath}}}\sigma_{i}
\end{displaymath}

Since the rel­a­tive or­der of the sub­scripts in the cy­cle was main­tained in the reno­ta­tion, the sums still con­tain the ex­act same three terms, just in a dif­fer­ent or­der. Take out the com­mon fac­tors;

\begin{displaymath}
H_{1,\mbox{\scriptsize spin-orbit}} =
\frac{{\rm i}}{4m_{\...
...}}}}(V-E_n){\widehat p}_{{\overline{\imath}}}\right]\sigma_{i}
\end{displaymath}

Now ac­cord­ing to the gen­er­al­ized canon­i­cal com­mu­ta­tor of chap­ter 4.5.4:

\begin{displaymath}
{\widehat p}_i (V-E_n) = (V-E_n) {\widehat p}_i - {\rm i}\hbar\frac{\partial (V-E_n)}{\partial r_i}
\end{displaymath}

where $E_n$ is a con­stant that pro­duces a zero de­riv­a­tive. So ${\widehat p}_{{\overline{\imath}}}$, re­spec­tively ${\widehat p}_{{\overline{\overline{\imath}}}}$ can be taken to the other side of $V-E_n$ as long as the ap­pro­pri­ate de­riv­a­tives of $V$ are added. If that is done, $(V-E_n){\widehat p}_{{\overline{\imath}}}{\widehat p}_{{\overline{\overline{\imath}}}}$ and $-(V-E_n){\widehat p}_{{\overline{\overline{\imath}}}}{\widehat p}_{{\overline{\imath}}}$ can­cel since lin­ear mo­men­tum op­er­a­tors com­mute. What is left are just the added de­riv­a­tive terms:

\begin{displaymath}
H_{1,\mbox{\scriptsize spin-orbit}} =
\frac{\hbar}{4m_{\rm...
...math}}}}}{\widehat p}_{{\overline{\imath}}}
\right]\sigma_{i}
\end{displaymath}

Note that the er­rant eigen­value $E_n$ mer­ci­fully dropped out. Now the hy­dro­gen po­ten­tial $V$ only de­pends on the dis­tance $r$ from the ori­gin, as 1/$r$, so

\begin{displaymath}
\frac{\partial V}{\partial r_i} = - \frac{V}{r^2}r_i
\end{displaymath}

and plug­ging that into the op­er­a­tor, you get

\begin{displaymath}
H_{1,\mbox{\scriptsize spin-orbit}} =
- \frac{\hbar V}{4m_...
...{\imath}}}}{\widehat p}_{{\overline{\imath}}}\right]\sigma_{i}
\end{displaymath}

The term be­tween the square brack­ets can be rec­og­nized as the $i$th com­po­nent of the an­gu­lar mo­men­tum op­er­a­tor; also the Pauli spin ma­trix $\sigma_i$ is de­fined as ${\widehat S}_i$$\raisebox{.5pt}{$/$}$${\textstyle\frac{1}{2}}\hbar$, so

\begin{displaymath}
H_{1,\mbox{\scriptsize spin-orbit}}
= - \frac{V}{2m_{\rm e}^2c^2r^2} \sum_i \L _i{\widehat S}_i
\end{displaymath}

Get rid of $c^2$ us­ing $\vert E_1\vert$ $\vphantom0\raisebox{1.5pt}{$=$}$ ${\textstyle\frac{1}{2}}\alpha^2{m_{\rm e}}c^2$, of $V$ us­ing $V$ $\vphantom0\raisebox{1.5pt}{$=$}$ $-2\vert E_1\vert a_0$$\raisebox{.5pt}{$/$}$$r$, and $m_{\rm e}$ us­ing $\vert E_1\vert$ $\vphantom0\raisebox{1.5pt}{$=$}$ $\hbar^2$$\raisebox{.5pt}{$/$}$$2{m_{\rm e}}a_0^2$ to get the spin-or­bit in­ter­ac­tion as claimed in the sec­tion on fine struc­ture.

That leaves the term

\begin{displaymath}
\sum_i {\widehat p}_i \frac{V-E_n}{4m_{\rm e}^2c^2}{\widehat p}_i
\end{displaymath}

in (D.59). Since $V$ $\vphantom0\raisebox{1.5pt}{$=$}$ $H_0-{\widehat p}^2/2m_{\rm e}$, it can be writ­ten as

\begin{displaymath}
\sum_i {\widehat p}_i \frac{H_0-E_n}{4m_{\rm e}^2c^2}{\wide...
...i
- \frac{\left({\widehat p}^{ 2}\right)^2}{8m_{\rm e}^3c^2}
\end{displaymath}

The fi­nal term is the claimed Ein­stein cor­rec­tion in the fine struc­ture Hamil­ton­ian, us­ing $\vert E_1\vert$ $\vphantom0\raisebox{1.5pt}{$=$}$ ${\textstyle\frac{1}{2}}\alpha^2{m_{\rm e}}c^2$ to get rid of $c^2$.

The first term,

\begin{displaymath}
H_{1,{\rm Darwin}} \stackrel{\mbox{?}}{=}
\sum_i {\widehat p}_i \frac{H_0-E_n}{4m_{\rm e}^2c^2}{\widehat p}_i
\end{displaymath}

is the sole re­main­ing prob­lem. It can­not be trans­formed into a de­cent phys­i­cal op­er­a­tor. The ob­jec­tive is just to get the en­ergy cor­rec­tion right. And to achieve that re­quires only that the Hamil­ton­ian per­tur­ba­tion co­ef­fi­cients are eval­u­ated cor­rectly at the $E_n$ en­ergy level. Specif­i­cally, what is needed is that

\begin{displaymath}
H_{\underline{\vec n}{\vec n},1,{\rm Darwin}} \equiv
\lang...
...t{\widehat p}_i(H_0-E_n){\widehat p}_i\psi_{{\vec n},0}\rangle
\end{displaymath}

for any ar­bi­trary pair of un­per­turbed hy­dro­gen en­ergy eigen­func­tions $\psi_{\underline{\vec n},0}$ and $\psi_{{\vec n},0}$ with en­ergy $E_n$. To see what that means, the lead­ing Her­mit­ian op­er­a­tor ${\widehat p}_i$ can be taken to the other side of the in­ner prod­uct, and in half of that re­sult, $H_0-E_n$ will also be taken to the other side:

\begin{displaymath}
H_{\underline{\vec n}{\vec n},1,{\rm Darwin}} = \frac{1}{8m...
...\vec n},0}\vert{\widehat p}_i\psi_{{\vec n},0}\rangle
\right)
\end{displaymath}

Now if you sim­ply swap the or­der of the fac­tors in $(H_0-E_n){\widehat p}_i$ in this ex­pres­sion, you get zero, be­cause both eigen­func­tions have en­ergy $E_n$. How­ever, swap­ping the or­der of $(H_0-E_n){\widehat p}_i$ brings in the gen­er­al­ized canon­i­cal com­mu­ta­tor $[V,{\widehat p}_i]$ that equals ${\rm i}\hbar\partial{V}$$\raisebox{.5pt}{$/$}$$\partial{r}_i$. There­fore, writ­ing out the re­main­ing in­ner prod­uct you get

\begin{displaymath}
H_{\underline{\vec n}{\vec n},1,{\rm Darwin}} = \frac{-\hba...
...}^*\psi_{{\vec n},0}}{\partial r_i}
{ \rm d}^3{\skew0\vec r}
\end{displaymath}

Now, the po­ten­tial $V$ be­comes in­fi­nite at $r$ $\vphantom0\raisebox{1.5pt}{$=$}$ 0, and that makes math­e­mat­i­cal ma­nip­u­la­tion dif­fi­cult. There­fore, as­sume for now that the nu­clear charge $e$ is not a point charge, but spread out over a very small re­gion around the ori­gin. In that case, the in­ner prod­uct can be rewrit­ten as

\begin{displaymath}
H_{\underline{\vec n}{\vec n},1,{\rm Darwin}} = \frac{-\hba...
...ec n},0}^*\psi_{{\vec n},0}
\right] { \rm d}^3{\skew0\vec r}
\end{displaymath}

and the first term in­te­grates away since $\psi_{\underline{\vec n},0}^*\psi_{{\vec n},0}$ van­ishes at in­fin­ity. In the fi­nal term, use the fact that the de­riv­a­tives of the po­ten­tial en­ergy $V$ give $e$ times the elec­tric field of the nu­cleus, and there­fore the sec­ond or­der de­riv­a­tives give $e$ times the di­ver­gence of the elec­tric field. Maxwell’s first equa­tion (13.5) says that that is $e$$\raisebox{.5pt}{$/$}$$\epsilon_0$ times the nu­clear charge den­sity. Now if the re­gion of nu­clear charge is al­lowed to con­tract back to a point, the charge den­sity must still in­te­grate to the net pro­ton charge $e$, so the charge den­sity be­comes $e\delta^3({\skew0\vec r})$ where $\delta^3({\skew0\vec r})$ is the three-di­men­sion­al delta func­tion. There­fore the Dar­win term pro­duces Hamil­ton­ian per­tur­ba­tion co­ef­fi­cients as if its Hamil­ton­ian is

\begin{displaymath}
H_{1,{\rm Darwin}} = \frac{\hbar^2e^2}{8m_{\rm e}^2c^2\epsilon_0} \delta^3({\skew0\vec r})
\end{displaymath}

Get rid of $c^2$ us­ing $\vert E_1\vert$ $\vphantom0\raisebox{1.5pt}{$=$}$ $\frac12\alpha^2{m_{\rm e}}c^2$, of $e^2$$\raisebox{.5pt}{$/$}$$\epsilon_0$ us­ing $e^2$$\raisebox{.5pt}{$/$}$$4\pi\epsilon_0$ $\vphantom0\raisebox{1.5pt}{$=$}$ $2\vert E_1\vert a_0$, and $m_{\rm e}$ us­ing $\vert E_1\vert$ $\vphantom0\raisebox{1.5pt}{$=$}$ $\hbar^2$$\raisebox{.5pt}{$/$}$$2{m_{\rm e}}a_0^2$ to get the Dar­win term as claimed in the sec­tion on fine struc­ture. It will give the right en­ergy cor­rec­tion for the non­rel­a­tivis­tic so­lu­tion. But you may rightly won­der what to make of the im­plied wave func­tion.