# 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))
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))
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))
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")