Use of simulation to explain concepts — PN junction to device scaling

import numpy as np
import matplotlib.pyplot as plt
import opennemo as nemo
from opennemo.visualization import style

style.use(notebook=True)
from scipy.linalg import solve_banded

props = nemo.structures.properties
T = 300.0
VT = nemo.core.units.KB * T                       # thermal voltage (V)
Q = nemo.core.units.Q                             # elementary charge (C)
EPS0 = nemo.core.units.EPS0                       # vacuum permittivity (F/m)
NC, NV = props.value("Si", "effective_density_conduction"), props.value("Si", "effective_density_valence")
EG, EPS = props.value("Si", "band_gap"), props.value("Si", "eps_static")
NI = np.sqrt(NC * NV) * np.exp(-EG / (2 * VT))   # intrinsic density (cm⁻³)
print(f"Si at 300 K: N_C = {NC:.3g}, N_V = {NV:.3g} cm⁻³, E_g = {EG:.4f} eV, n_i = {NI:.3e} cm⁻³")

def caughey_thomas(n, mu_min, mu_max, n_ref, alpha):
    return mu_min + (mu_max - mu_min) / (1 + (n / n_ref) ** alpha)


def bernoulli(z):
    z = np.asarray(z, float)
    small = np.abs(z) < 1e-6
    return np.where(small, 1 - z / 2, z / np.expm1(np.where(small, 1.0, z)))


def tridiag_solve(lo, diag, up, rhs, left, right):
    # banded solve with Dirichlet values at both ends
    a = np.zeros((3, len(diag)))
    a[0, 1:], a[1], a[2, :-1] = up, diag, lo
    a[1, 0], a[0, 1], rhs[0] = 1.0, 0.0, left
    a[1, -1], a[2, -2], rhs[-1] = 1.0, 0.0, right
    return solve_banded((1, 1), a, rhs)


class Diode:
    def __init__(self, na, nd, lp=6000.0, ln=6000.0, dx=5.0, tau=1e-7):
        self.x = np.arange(0.0, lp + ln + dx / 2, dx)           # nm
        self.lp, self.na, self.nd, self.tau = lp, na, nd, tau
        p_side = self.x < lp
        doping = np.where(p_side, na, nd)
        self.mun = caughey_thomas(doping, 68.5, 1414.0, 9.2e16, 0.711)   # cm²/Vs
        self.mup = caughey_thomas(doping, 44.9, 470.5, 2.23e17, 0.719)
        self.poisson = nemo.Poisson(nemo.Grid(self.x), EPS, np.where(p_side, -na, nd), T)
        self.ends = np.zeros(len(self.x), bool)
        self.ends[[0, -1]] = True
        self.psi_p, self.psi_n = -VT * np.log(na / NI), VT * np.log(nd / NI)
        self.psi = np.where(p_side, self.psi_p, self.psi_n)
        self.phin, self.phip = np.zeros(len(self.x)), np.zeros(len(self.x))

    def carriers(self):
        return NI * np.exp((self.psi - self.phin) / VT), NI * np.exp((self.phip - self.psi) / VT)

    def continuity(self, v):
        x, d = self.x * 1e-7, np.diff(self.psi) / VT               # cm; potential step in units of VT
        h = np.diff(x)
        vol = np.zeros(len(x))
        vol[1:] += h / 2
        vol[:-1] += h / 2
        n, p = self.carriers()
        den = self.tau * (n + NI) + self.tau * (p + NI)
        dn = VT * 0.5 * (self.mun[1:] + self.mun[:-1]) / h
        dp = VT * 0.5 * (self.mup[1:] + self.mup[:-1]) / h
        bp, bm = bernoulli(d), bernoulli(-d)
        # electrons: F_{i+1/2} = D/h [n_{i+1} B(d) - n_i B(-d)],  F_{i+1/2} - F_{i-1/2} = R vol
        diag = np.zeros(len(x))
        diag[:-1] -= dn * bm
        diag[1:] -= dn * bp
        n_new = tridiag_solve(dn * bm, diag - p / den * vol, dn * bp, -NI**2 / den * vol,
                              NI * np.exp((self.psi[0] - v) / VT), NI * np.exp(self.psi[-1] / VT))
        # holes: G_{i+1/2} = D/h [p_i B(d) - p_{i+1} B(-d)],  G_{i+1/2} - G_{i-1/2} = -R vol
        diag = np.zeros(len(x))
        diag[:-1] += dp * bp
        diag[1:] += dp * bm
        p_new = tridiag_solve(-dp * bp, diag + n_new / den * vol, -dp * bm, NI**2 / den * vol,
                              NI * np.exp((v - self.psi[0]) / VT), NI * np.exp(-self.psi[-1] / VT))
        self.phin = self.psi - VT * np.log(np.maximum(n_new, 1e-300) / NI)
        self.phip = self.psi + VT * np.log(np.maximum(p_new, 1e-300) / NI)
        jn = Q * dn * (n_new[1:] * bp - n_new[:-1] * bm)
        jp = Q * dp * (p_new[:-1] * bp - p_new[1:] * bm)
        return jn + jp                                             # A/cm² between the nodes

    def solve(self, v, tol=1e-9, max_iter=300):
        self.poisson.set_dirichlet(self.ends, np.where(self.x < self.x.mean(), self.psi_p + v, self.psi_n))
        self.phin[0] = self.phip[0] = v
        for it in range(max_iter):
            n, p = self.carriers()
            new = self.poisson.solve(self.psi, electrons=n, holes=p).potential
            change, self.psi = np.abs(new - self.psi).max(), new
            j = self.continuity(v)
            if change < tol and it > 1:
                break
        return j.mean(), np.ptp(j) / abs(j.mean()), it + 1

    def iv(self, volts):
        return np.array([self.solve(v)[0] for v in volts])

for na, nd in ((1e15, 1e15), (1e17, 1e17)):
    d = Diode(na, nd)
    d.solve(0.0)
    vbi = VT * np.log(na * nd / NI**2)
    w = np.sqrt(2 * EPS * EPS0 * vbi / (Q * 1e6) * (na + nd) / (na * nd)) * 1e9    # nm
    field = np.abs(np.gradient(d.psi, d.x * 1e-7))                                   # V/cm
    n, p = d.carriers()
    rho = p - n + np.where(d.x < d.lp, -na, nd)                    # space charge / q (cm⁻³)
    w_sim = np.trapezoid(np.abs(rho), d.x) / na                    # nm: the depleted width it amounts to
    print(f"N_A = N_D = {na:.0e}: V_bi {np.ptp(d.psi):.4f} V (formula {vbi:.4f});  W = ∫|ρ|dx/qN "
          f"{w_sim:.0f} nm (formula {w:.0f});  E_max {field.max():.3e} V/cm "
          f"(formula {2 * vbi / (w * 1e-7):.3e})")

volts = np.round(np.arange(0.0, 0.901, 0.05), 3)
base = Diode(1e15, 1e15)
j_base = base.iv(volts)

na = nd = 1e15
dn, dp = VT * base.mun[0], VT * base.mup[-1]
vbi = VT * np.log(na * nd / NI**2)
w = np.sqrt(2 * EPS * EPS0 * vbi / (Q * 1e6) * 2 / na) * 1e2                        # cm
wp = wn = 6e-4 - w / 2
j0 = Q * NI**2 * (dn / (na * wp) + dp / (nd * wn))
print(f"D_n = {dn:.1f}, D_p = {dp:.1f} cm²/s, diffusion length sqrt(D_n tau) = "
      f"{np.sqrt(dn * 1e-7) * 1e4:.0f} µm, J_0 = {j0:.3e} A/cm²")
print(" V (V)   J simulation (A/cm²)   J Shockley (A/cm²)   ratio")
for v, j in zip(volts, j_base):
    if v in (0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9):
        shockley = j0 * np.expm1(v / VT)
        print(f"{v:5.2f}   {j:18.3e}   {shockley:18.3e}   {j / shockley:6.2f}")
ideality = np.gradient(volts, np.log(np.maximum(j_base, 1e-30))) / VT
print("local ideality factor n = (q/kT) dV/dlnJ at 0.2, 0.5, 0.8 V:",
      np.round(ideality[[4, 10, 16]], 2))

def bands(d, v):
    d.solve(v)
    ec = -d.psi + EG / 2 + VT / 2 * np.log(NV / NC)       # conduction band edge (eV), E_i = -q psi
    return ec, ec - EG, -d.phin, -d.phip


def junction_voltage(d, v):
    # quasi-Fermi level splitting at the junction: n p = n_i² exp(q V_j / kT) there
    d.solve(v)
    k = np.argmin(np.abs(d.x - d.lp))
    return d.phip[k] - d.phin[k]


fig, axes = plt.subplots(1, 2, figsize=(11, 4))
cols = style.sweep(4)
for ax, v in zip(axes, (0.5, 0.9)):
    ec, ev, efn, efp = bands(base, v)
    um = base.x / 1e3
    ax.plot(um, ec, color=style.OPEN_NEMO, label="$E_c$")
    ax.plot(um, ev, color=style.REFERENCE, label="$E_v$")
    ax.plot(um, efn, color=cols[1], ls="--", label="$E_{Fn}$")
    ax.plot(um, efp, color=cols[2], ls=":", label="$E_{Fp}$")
    ax.set(xlabel="position (µm)", ylabel="energy (eV)", title=f"$N_A = N_D = 10^{{15}}$ cm⁻³, V = {v} V")
    ax.legend(fontsize=8)
fig.tight_layout()
for v in (0.5, 0.6, 0.7, 0.8, 0.9):
    vj = junction_voltage(base, v)
    print(f"V = {v}: junction voltage {vj:.3f} V, dropped in the neutral regions {v - vj:.3f} V")

designs = {"10¹⁵ cm⁻³, 6 + 6 µm": j_base,
           "10¹⁷ cm⁻³, 6 + 6 µm": Diode(1e17, 1e17).iv(volts),
           "10¹⁵ cm⁻³, 1 + 1 µm": Diode(1e15, 1e15, 1000.0, 1000.0, dx=2.0).iv(volts),
           "10¹⁷ cm⁻³, 1 + 1 µm": Diode(1e17, 1e17, 1000.0, 1000.0, dx=2.0).iv(volts)}
print(" V (V)  " + "  ".join(f"{k:>19s}" for k in designs))
for i, v in enumerate(volts):
    if v in (0.3, 0.5, 0.6, 0.7, 0.8, 0.9):
        print(f"{v:5.2f}  " + "  ".join(f"{j[i]:19.3e}" for j in designs.values()))
fig, ax = plt.subplots(figsize=(7, 4.2))
for c, (label, j) in zip(style.sweep(len(designs)), designs.items()):
    ax.semilogy(volts[1:], j[1:], color=c, label=label)
ax.set(xlabel="forward bias V (V)", ylabel="current density J (A/cm²)")
ax.legend(fontsize=8)
fig.tight_layout()
r = designs["10¹⁷ cm⁻³, 6 + 6 µm"] / designs["10¹⁵ cm⁻³, 6 + 6 µm"]
print(f"\ndoping x100 changes J by x{r[volts == 0.5][0]:.3f} at 0.5 V and x{r[volts == 0.9][0]:.2f} at 0.9 V")