Binding energy from Dirac equation as a function of Z

Here we shall plot the electron binding energy predicted by the Dirac equation for the hydrogen atom:

\[E_{nj}=mc^2\left\{\left[1+\left(\frac{\alpha Z}{n-(j+1/2)+\sqrt{(j+1/2)^2-(\alpha Z)^2}}\right)^2\right]^{\frac{-1}{2}}-1\right\}\]

import matplotlib as mpl

params = {
    "font.family": "serif",
    "text.usetex": True,
    "font.serif": 'Palatino',
    "figure.figsize": [5.5,4],
}
mpl.rcParams.update(params)

import numpy as np
from matplotlib import pyplot as plt

# Setting a nice plot style
plt.style.use('ggplot')

# Define a function to calculate E(Z) given j and n

def E(Z,j,n):
    
    # Defining some constants 
    
    c = 3*(10**8) # Speed of light
    m = 9.1*(10**(-31)) # Mass of the electron
    alpha = 137**-1 # Fine structure constant 
    
    # Defining our binding energy function, E(Z) 
    # derived from the Dirac equation.
    
    denom = n - (j + .5) + np.sqrt((j + 0.5)**2 - (alpha*Z)**2)
    frac = alpha*Z/denom
    return m*(c**2)*((1+frac**2)**(-0.5) - 1)  

# Now defining a function to plot E(Z)

def binding_energy_plot(j,n):
    
    # Create an even spaced array of numbers in the range 
    # 0 - 137, in steps of 0.1.
    Z = np.arange(0.0, 137.0, 0.1)  
    
    
    # Plot the function
    plt.plot(Z, E(Z,j,n))
    
    # Label the x and y axies
    plt.xlabel("Z")
    plt.ylabel("Binding Energy (Joules)")
    
binding_energy_plot(0.5,1)     
plt.tight_layout()
plt.savefig('dirac_binding_energy.png')

We can see that the binding energy goes asymptotically to \(-\infty\) as \(Z\rightarrow 137\) (that is, as \(\alpha Z\rightarrow 1\)). In this regime, the coupling is so strong that we can have spontaneous production of electron-positron pairs from the vacuum, and we must turn to quantum field theory instead.