RasmahRasmah

Semiconductor Drift-Diffusion

Overview

A semiconductor at thermal equilibrium sets up an electrostatic potential that makes the electron and hole currents vanish. The carrier densities follow Boltzmann statistics (Selberherr 1984),

\[n = n_i e^{\phi/V_T}, \qquad p = n_i e^{-\phi/V_T}, \qquad V_T = \frac{kT}{q}\]

and the potential solves the nonlinear Poisson equation

\[-\nabla\cdot(\varepsilon \nabla\phi) = q\,(p - n + N_D^+ - N_A^-)\]

where n_i is the intrinsic carrier density, V_T the thermal voltage, and N_D⁺ − N_A⁻ the net ionized doping. Rasmah solves this zero-bias Poisson–Boltzmann problem with Newton.

using Rasmah

A manufactured solution

The strongest test is to manufacture an exact solution: choose φ(x,y,z) = (1−x²)(1−y²)(1−z²) on the cube [−1,1]³ and pick the source f that makes it solve the equation exactly. The Newton solve must recover it.

ε, q, ni, VT = 1.0, 1.0, 1.0, 1.0
doping(x) = 0.0

# φ = (1−x²)(1−y²)(1−z²) is zero on the boundary and has ∇²φ = −2·sum of pairwise products
∇²φ(x) = -2 * ((1 - x[2]^2) * (1 - x[3]^2) + (1 - x[1]^2) * (1 - x[3]^2) + (1 - x[1]^2) * (1 - x[2]^2))
φ_exact(x) = (1 - x[1]^2) * (1 - x[2]^2) * (1 - x[3]^2)
f(x) = ε * ∇²φ(x) / q + ni * exp(-φ_exact(x) / VT) - ni * exp(φ_exact(x) / VT)

mesh = tetrahedralize(brep(box(2.0, 2.0, 2.0)), 5)

n = size(mesh.nodes, 2)
fixed = Int[]; g = Float64[]
for i in 1:n
    x, y, z = mesh.nodes[1, i], mesh.nodes[2, i], mesh.nodes[3, i]
    if abs(x) > 1 - 1e-6 || abs(y) > 1 - 1e-6 || abs(z) > 1 - 1e-6
        push!(fixed, i); push!(g, 0.0)
    end
end

φ = solve_semiconductor(mesh, ε, q, ni, VT, doping; fixed = fixed, g = g, f = f)
216-element Vector{Float64}:
 0.0
 0.0
 0.0
 0.0
 0.0
 0.0
 0.0
 0.0
 0.0
 0.0
 ⋮
 0.0
 0.0
 0.0
 0.0
 0.0
 0.0
 0.0
 0.0
 0.0

Verification — the potential is recovered

The interior potential must match φ_exact. At the cube centre it is 1, and the linear-element solve recovers it to within discretisation error:

center = argmin(sum(abs2, mesh.nodes[:, i]) for i in 1:n)
println("φ at centre = ", round(φ[center]; digits = 4), "  (analytic = 1.0)")
φ at centre = 0.8637  (analytic = 1.0)

Next steps