Semiconductor Physics Equations

README.rst

1D Semiconductor Physics Equations

Drift-diffusion Model

Equations are here (4):

http://certik.github.com/theoretical-physics/book/src/elmag/elmag.html#equations

This is called a drift-diffusion model and can model any semiconductor. For example pn-junctions:

http://en.wikipedia.org/wiki/P-n_junction

transistors:

http://en.wikipedia.org/wiki/Transistor

or thin films of silicon, very simplified model is here:

http://certik.github.com/theoretical-physics/book/src/elmag/elmag.html#example-1

The last two equations (and read the assumptions above that had to be made). The 3 nonlinear PDEs is a way more general and precise model, that also gives the charges at the boundaries etc.

1D Case

For the 1D pn-junction, the equations are here:

http://certik.github.com/theoretical-physics/book/src/elmag/elmag.html#example-3

it's a set of 6 nonlinear first-order ODEs. The boundary conditions are not exactly clear to me yet, but here is what I should get:

http://www.nextnano.de/nextnano3/tutorial/1Dtutorial_pn_junction.htm

so I tried to solve it using the Euler method (for ODEs), see the attached python script. But unfortunately, p(x) explodes, but it should go to zero, as x->320nm. So either the equations are wrong, or I am setting the model constants wrong.

The constants (in the code) are:

Constants (SI units)

Basic:

q: 1.6022e-19
k_b: 1.3806503e-23
eps_0: 8.854187817e-12
T: 298.0

GaAs:

eps_r: 12.93
mu_p: 0.045
mu_n: 0.14
D_p: 0.00115556862583
D_n: 0.00359510239146
N_A: 5e+23
N_D: 2e+24

semiconductors.py

#from scipy.integrate import odeint
from numpy import arange, pi, arctan, array, zeros
atan = arctan
from pylab import plot, savefig, legend

def euler(derivs, y_0, t_list):
    """
    It has the same interface as odeint, but uses a simple Euler method, so
    that it's robust.
    """
    R = zeros((len(t_list), len(y_0)))
    y = array(y_0)
    R[0, :] = y
    for i, t in enumerate(t_list[:-1]):
        delta_t = t_list[i+1]-t_list[i]
        yp = array(derivs(y, t))
        y += yp*delta_t
        R[i+1, :] = y
    return R


# everything in SI units:
nm = 1e-9
cm = 1e-2

q = 1.6022e-19
k_b = 1.3806503e-23
T = 298.
eps_0 = 8.854187817e-12
R = 0.

N_A = 0.5*1e18 * cm**(-3)
N_D = 2.0*1e18 * cm**(-3)
eps_r = 12.93
eps = eps_r*eps_0
mu_n = 1400 * cm**2
mu_p = 450 * cm**2
D_n = mu_n * k_b * T/q
D_p = mu_p * k_b * T/q
print "1D pn-junction simulation"
print
print "Constants (SI units)"
print "--------------------"
print "Basic:"
print "q:", q
print "k_b:", k_b
print "eps_0:", eps_0
print "T:", T
print
print "GaAs:"
print "eps_r:", eps_r
print "mu_p:", mu_p
print "mu_n:", mu_n
print "D_p:", D_p
print "D_n:", D_n
print "N_A:", N_A
print "N_D:", N_D

def C(t):
    return (N_D+N_A)*(atan((t-150*nm)/(10*nm))/pi +0.5) - N_A

def deriv(y, t):
    yp = [0, 0, 0, 0, 0, 0]
    yp[5] = q/eps * (y[2] - y[0] + C(t))
    yp[0] = y[1]
    yp[1] = R/D_n+mu_n/D_n * (y[1]*y[5]+y[0]*yp[5])
    yp[2] = y[3]
    yp[3] = R/D_p-mu_p/D_p * (y[3]*y[5]+y[2]*yp[5])
    yp[4] = y[5]
    return yp

t = arange(0, 200*nm, 300*nm/2000)
_eps = N_A/1.e3
# The n(0) condition is a problem, as we don't now it, and if we set it to 0,
# it will always be 0:
y = euler(deriv, [-1e-20, 0., N_A, 0., 0., 0.], t)
n = y[:, 0]
n_p = y[:, 1]
p = y[:, 2]
p_p = y[:, 3]
phi = y[:, 4]
phi_p = y[:, 5]
plot(t, n, "r-", lw=2, label="n")
plot(t, p, "b-", lw=2, label="p")
#plot(t, p_p, "bx", label="p_p")
plot(t, phi, "k-", label="phi")
plot(t, C(t), "g-", lw=2, label="C")
legend()
savefig("sol.png")

sol.png