SEM Dopant Contrast due to Patch Fields in Semiconductors

Re-implementation of MatSci Part III research project, originally done as a MuPad program, now re-done for fun using Python tools.

The project did a simulation of electrons that are emitted when an electron beam hits a p-n junction in a scanning electron microscope.

This file is basically a translation of my original MuPad source code, but using features of IPython:

  • embedded LaTeX, using equations borrowed from my final report
  • generating inline some of the plots I had to generate in other ways for my final report.
  • doing everything in Python, instead of a mixture of MuPad and Perl scripts.

It isn't complete, but most of the work is done. For reference and explanation you'll need to see my original 2002 report.

In [1]:
import math
import numpy as np
from scipy import integrate
import sympy
In [2]:
# Physical constants

q_e = -1.60e-19 # charge of electron
m_e = 9.11e-31 # mass of electron
In [3]:
# Simulation constants. SI units unless otherwise stated
working_distance = 4.5e-3 # Electron microscope geometry
num_events = 400
num_energy_bins = 20
num_angle_bins = 10
t_finish = 2e-9
y0 = 1e-12 # arbitrary
V_BI_base = 0.7 # Corresponds to band gap for silicon


# parameters that we will vary, but have some defaults here for sake of graphs
work_function_base = 4.52 # electron volts
V_extract_base = 20.0
excitation_point_base = 0.0

Probability density for angle/energy of electrons:

Standard results for back-scattered electrons in this configuration:

$$ N(E) \propto \frac{E}{(E+\phi)^{4}} $$$$ N(\alpha) \propto cos(\alpha) $$
In [4]:
def make_prob_distributions(work_function):
    prob_angle = lambda a: np.cos(a)
    prob_energy = lambda en: en/(en + work_function) ** 4
    return prob_angle, prob_energy
    
angle_min = - math.pi / 2
angle_range = math.pi
angle_max = angle_min + angle_range
energy_min = 0
# energy_range and energy_max to be calculated.

Normalise:

$$ n_1 = \int_0^\infty p_{energy}(e) {d}e $$$$ n_2 = \int_{-\pi/2}^{\pi/2} p_{angle}(\alpha) {d}\alpha $$
In [5]:
def make_normalised_probability_distribution(prob_angle, prob_energy):
    n1, _ = integrate.quad(prob_energy, 0, 10^9)
    n2, _ = integrate.quad(prob_angle, -np.pi/2, np.pi/2)

    return lambda angle, en: prob_energy(en) * prob_angle(angle) / (n1 * n2)

Check that things look the way they should:

In [6]:
prob_angle_tmp, prob_energy_tmp = make_prob_distributions(work_function_base)
prob_tmp = make_normalised_probability_distribution(prob_angle_tmp, prob_energy_tmp)


%matplotlib inline
from matplotlib import pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

angle_axis = np.linspace(-np.pi/2, np.pi/2, 50)
energy_axis = np.linspace(0, 20, 100)
X, Y = np.meshgrid(angle_axis, energy_axis)
Z = prob_tmp(X, Y)
fig = plt.figure(figsize=(12, 9))
ax = plt.subplot(111, projection='3d')
ax.view_init(30, 30)
ax.plot_surface(X, Y, Z, cstride=2, rstride=5, cmap=plt.cm.coolwarm)
plt.title("Probability of electron")
ax.set_xlabel("Angle")
ax.set_ylabel("Energy (eV)")
None

Split into buckets

We want to adjust the energy range so that we get a self-similar distribution of energies whatever the work function. We scale it so that we have at most 1 electron in the box that is at $e = energyMax, \alpha=0$.

We require zero electrons in the box that would be at $e=energyMax$, $\alpha=0$. This requires:

$$ prob(0, energyMax) . numEvents . \frac{energyMax}{numEnergyBins} . \frac{angleRange}{numAngleBins} < 1/2 $$

since < 0.5 would be rounded down to 0.

Using a symbolic integration on $n_1$ above, it can be shown that the equation we need to solve is:

$$ 6 (energyMax . workFunction)^2 . numEvents . angleMax - numEnergyBins . numAngleBins .(energyMax + workFunction)^4 = 0 $$

There are 4 solutions to this, 2 negative and 2 positive, and we need the largest positive one.

In [7]:
energyMax, angleRange, numEnergyBins, numAngleBins, workFunction, numEvents = sympy.symbols("energyMax angleRange numEnergyBins numAngleBins workFunction numEvents")
In [8]:
energyMaxSolutions = sympy.solve(6 * (energyMax * workFunction)**2 * numEvents * angleRange - numEnergyBins * numAngleBins * (energyMax + workFunction)**4, energyMax)
In [9]:
def get_energy_max(work_function):
    return float(max([s.evalf(subs=dict(workFunction=work_function,
                                  angleRange=angle_range,
                                  numAngleBins=num_angle_bins,
                                  numEnergyBins=num_energy_bins,
                                  numEvents=num_events))
                      for s in energyMaxSolutions]))

Build probability distribution:

In [10]:
def make_energy_bins(work_function):
    energy_max = get_energy_max(work_function)
    prob_angle, prob_energy = make_prob_distributions(work_function)
    prob = make_normalised_probability_distribution(prob_angle, prob_energy)
    
    d_e = (energy_max - energy_min)/ num_energy_bins
    d_a = angle_range / num_angle_bins
    bins = np.empty([num_energy_bins, num_angle_bins])

    for i in range(0, num_energy_bins):
        for j in range(0, num_angle_bins):
            # Really should do double integration over prob(),
            # but a good enough numerical approximation is to
            # average height at the 4 corners and multiple by base area
            v = (prob(angle_min + j * d_a,       energy_min + i * d_e) +
                 prob(angle_min + (j + 1) * d_a, energy_min + i * d_e) +
                 prob(angle_min + j * d_a,       energy_min + (i + 1) * d_e) +
                 prob(angle_min + (j + 1) * d_a, energy_min + (i + 1) * d_e)) / 4
            # Scale by num_events        
            bins[i, j] = np.round(float(v * d_a * d_e * num_events))
    return bins, (angle_min, angle_max), (energy_min, energy_max)
In [11]:
# Plot it:

bins_tmp, _, _ = make_energy_bins(work_function_base)

def hist_3d(matrix):
    fig = plt.figure(figsize=(12, 9))
    ax = fig.add_subplot(111, projection='3d')
    x_data, y_data = np.meshgrid(np.arange(matrix.shape[1]),
                                 np.arange(matrix.shape[0]) )
    x_data = x_data.flatten()
    y_data = y_data.flatten()
    z_data = matrix.flatten()
    ax.bar3d(x_data, y_data, np.zeros(len(z_data)),
             1, 1, z_data)
    return ax
ax = hist_3d(bins_tmp)
ax.view_init(30, 30)

Build a list of angles and energies of electrons according to the distribution above. Each electron is given an energy randomly selected from within the bit it is in.

In [12]:
def generate_electrons(work_function):
    np.random.seed(1234) # same each time
    r = np.random.random_sample
    bins, (angle_min, angle_max), (energy_min, energy_max) = make_energy_bins(work_function)
    d_a = (angle_max - angle_min) / num_angle_bins
    d_e = (energy_min - energy_max) / num_energy_bins
    return [(float(angle_min + (j + r()) * d_a),
             float(energy_min + (i + r()) * d_e))
             for i in range(0, num_energy_bins) for j in range(0, num_angle_bins) for k in range(0, int(bins[i, j]))]

Electric potential

The form of the potential has been calculated theoretically:

$$ V(x, y) = \frac{V_{BI}}{\pi} arctan(\frac{x}{y}), y > 0 $$

The extractor bias $V_{extract}$ and thus the extractor field $E_{extract}$ also needs to be added to this:

$$ V(x,y) = \frac{V_{BI}}{\pi} \arctan \left(\frac{x}{y}\right) - E_{extract}y $$
In [13]:
# Symbolic manipulation:
x, y, V_BI, V_extract, workingDistance = sympy.symbols("x y V_BI V_extract workingDistance")
E_extract = -1 * (V_extract / workingDistance)
V = V_BI / sympy.pi * sympy.atan(x/y) - E_extract * y
E_x = -1 * sympy.diff(V, x)
E_y = -1 * sympy.diff(V, y)
E_z = 0

from sympy.interactive import printing
printing.init_printing(use_latex=True)
E_x, E_y
/home/luke/tmpstore/anaconda/lib/python2.7/site-packages/IPython/core/formatters.py:239: FormatterWarning: Exception in image/png formatter: 
\begin{pmatrix}- \frac{V_{BI}}{\pi y \left(\frac{x^{2}}{y^{2}} + 1\right)}, & \frac{V_{BI} x}{\pi y^{2} \left(\frac{x^{2}}{y^{2}} + 1\right)} - V_{extract} / workingDistance\end{pmatrix}
^
Unknown symbol: \begin (at char 0), (line:1, col:1)
  FormatterWarning,
Out[13]:
$$\begin{pmatrix}- \frac{V_{BI}}{\pi y \left(\frac{x^{2}}{y^{2}} + 1\right)}, & \frac{V_{BI} x}{\pi y^{2} \left(\frac{x^{2}}{y^{2}} + 1\right)} - \frac{V_{extract}}{workingDistance}\end{pmatrix}$$
In [14]:
# As Python expression:
printing.init_printing(pretty_print=False, use_latex=False)
E_x, E_y
Out[14]:
(-V_BI/(pi*y*(x**2/y**2 + 1)), V_BI*x/(pi*y**2*(x**2/y**2 + 1)) - V_extract/workingDistance)
In [15]:
def make_field(V_BI, V_extract):
    # Method 1: could use sympy.evalf:
    # def E_field(x, y):
    #    subs = dict(workingDistance=working_distance, V_BI=V_BI, V_extract=V_extract, x=x, y=y)
    #    return (E_x.evalf(subs=subs),
    #            E_y.evalf(subs=subs))
    
    # For speed, though, we copy-paste output from above and fix where needed:
    def E_field(x, y):
        return (-V_BI/(math.pi*y*(x**2/y**2 + 1)),
                 V_BI*x/(math.pi*y**2*(x**2/y**2 + 1)) - V_extract/working_distance)
    
    return E_field
In [16]:
#  Plot it:
   
E_field_tmp = make_field(V_BI_base, V_extract_base)
E_field_x = lambda x, y: E_field_tmp(x, y)[0]
E_field_y = lambda x, y: E_field_tmp(x, y)[1]
In [17]:
fig = plt.figure(figsize=(15,8))
ax = fig.add_subplot(111, projection='3d')
x = np.linspace(-0.0001, 0.0001, num=50)
y = np.linspace(0.0000001, 0.0002, num=50)
X, Y = np.meshgrid(x, y)
Z = E_field_x(X, Y).clip(-30000, 0)
plt.title("$E_x$")
ax.plot_surface(X, Y, Z, cstride=2, rstride=2, cmap=plt.cm.Reds)
ax.set_xlabel("x (m)")
ax.set_ylabel("y (m)")
None
In [18]:
fig = plt.figure(figsize=(15,10))
ax = fig.add_subplot(111, projection='3d')
x = np.linspace(-4e-10, 4e-10, num=50)
y = np.linspace(1e-11, 5e-10, num=50)
X, Y = np.meshgrid(x, y)
Z = E_field_y(X, Y)
Z = Z.clip(-1.5e9, 1.5e9) # clip to avoid singularity
plt.title("$E_y$")
ax.plot_surface(X, Y, Z, cstride=2, rstride=2, cmap=plt.cm.Greens)
ax.set_xlabel("x (m)")
ax.set_ylabel("y (m)")
None

Equations of motion

Equation of motion for electron:

$$ m_e \frac{{d}^2\mathbf{r}}{{d}t^2} = q_e\mathbf{E}(\mathbf{r}) $$

where

$$ \mathbf{r} = \left( \begin{array}{l}x\\y\end{array} \right) $$

We can convert this second order o.d.e in 2 dependent variables to a first order o.d.e. in 4 equations:

$$ \frac{d}{dt}\left( \begin{array}{l} x\\x'\\y\\y'\\ \end{array} \right) = \left( \begin{array}{l} x'\\\frac{\displaystyle q_e}{\displaystyle m_e}E_x(x,y)\\y'\\\frac{\displaystyle q_e}{\displaystyle m_e}E_y(x,y) \end{array} \right)\\ $$
In [19]:
def make_ode(V_BI, V_extract):
    E_field = make_field(V_BI, V_extract)
    def ode_func(t, vec):
        x, y, vx, vy = vec
        E_x, E_y = E_field(x, y)
        return np.array([vx,
                         vy,
                         q_e/m_e * E_x,
                         q_e/m_e * E_y])
    return ode_func

Initial conditions - initial position and velocity

In [20]:
ev_to_joule = lambda energy_ev: abs(q_e) * energy_ev
v_x_0 = lambda angle, energy: np.sqrt(2 * np.abs(energy) / m_e) * np.sin(angle)
v_y_0 = lambda angle, energy: np.sqrt(2 * np.abs(energy) / m_e) * np.cos(angle)
In [21]:
def make_starting_point(electron_angle, electron_energy, excitation_point):
    energy = ev_to_joule(electron_energy)
    return np.array([excitation_point, # The point where electron beam is hitting the silicon, with p-n junction at origin
                     y0,
                     v_x_0(electron_angle, energy),
                     v_y_0(electron_angle, energy)
                    ])

Detection routine

In [22]:
# Electron micropscope geometry:
DETECTION_HEIGHT = working_distance
DETECTION_WIDTH = 2.5e-3 

t_finish_last = 2e-9

# Use a search to try to home in on the point when
# electron will be at the detector height.
def detect(starting_point, ode_func, return_trajectory=False):
    global t_finish_last
    t_min = 0
    num_guesses = 0
    t_guess = t_finish_last
    t_max = None
    t_last_min = None
    t_give_up = 0.01
    detected = False
    while True:
        num_guesses += 1
        
        def f(t_final):
            r = integrate.ode(ode_func)
            r.set_integrator('vode')
            r.set_initial_value(starting_point, 0)
            return r.integrate(t_final)
        
        x, y, vx, vy = f(t_guess)
        if y > DETECTION_HEIGHT:
            if y/DETECTION_HEIGHT > 1.1:
                # Need to go back
                t_max = t_guess
                t_guess = (t_min + t_max)/2
            else:
                # Right height
                
                if abs(x) < DETECTION_WIDTH:
                    # Correct width - caught
                    detected = True
                    #print "FOUND after {0}".format(num_guesses)
                    
                else:
                    detected = False
                    #print "NOTFOUND after {0}".format(num_guesses)
                
                # Either way, we're done
                t_finish_last = t_guess
                
                break
        else:
            # Need to go higher
            if t_guess > t_give_up:
                #print "ABANDONED after {0}".format(num_guesses)
                break
            if t_guess > t_min:
                t_min = t_guess
            if t_max is None: # We have no idea what to use for t_max
                t_guess = t_guess * 2
            else:
                t_guess = (t_min + t_max) /2
    if return_trajectory:
        return detected, np.array([f(i) for i in np.arange(0, t_guess, t_guess/20)])

    return detected
In [23]:
V_extract_test = 20.0
electrons_test = generate_electrons(work_function_base)
e_angle, e_energy = electrons_test[0]
ode_func_test = make_ode(V_BI_base, V_extract_test)
ic = make_starting_point(e_angle, e_energy, 0)
In [24]:
# t_finish_last = 1e-12
# %timeit detect(ic, ode_func_test)
In [148]:
# Sanity check: plot some trajectories

V_extract_test = 20.0
electrons_test = generate_electrons(work_function_base)
ode_func_test = make_ode(V_BI_base, V_extract_test)

fig = plt.figure(figsize=(10, 12))

for i, (e_angle, e_energy) in enumerate(electrons_test[::4]):
    ic = make_starting_point(e_angle, e_energy, 0)
    detected, path = detect(ic, ode_func_test, return_trajectory=True)
    xs = path[:, 0].clip(-0.004, 0.004)
    ys = path[:, 1].clip(0, 0.005)
    plt.plot(xs, ys, color="red" if detected else "blue")