#!/usr/bin/env python3 # -*- coding: utf-8 -*- """ hetero_pn_band_finite.py Finite-thickness abrupt hetero-pn junction under the depletion approximation. Coordinate ---------- x = 0 : n-layer outer contact/surface x = d_n : n/p heterointerface x = d_n + d_p : p-layer outer contact/surface Bias sign --------- V > 0 : forward bias V < 0 : reverse bias The electrostatic junction drop is Vj = Vbi - V Energy reference ---------------- EV_n, EV_p, EF_n, EF_p before contact must be specified on one common energy reference. Usually EV_n = 0 eV is convenient. Then Vbi [V] = EF_n - EF_p (energies in eV) Finite-thickness treatment -------------------------- This program treats four regimes: 1) n partial, p partial 2) n fully depleted, p partial 3) n partial, p fully depleted 4) n fully depleted, p fully depleted When both sides are only partially depleted, the ordinary depletion-charge neutrality condition applies: q ND Wn = q NA Wp If one side becomes fully depleted, additional bias is supported partly by charge on the external contact of the fully depleted layer. Therefore, after the onset of full depletion, q ND d_n = q NA Wp (or the p-side counterpart) is NOT generally maintained. It is true only at the onset where the electric field at that external surface is zero. Assumptions ----------- - abrupt heterointerface - uniform fully ionized dopants - depletion approximation - no fixed sheet charge at the heterointerface - ideal external contacts that can supply whatever surface charge is needed to satisfy the imposed terminal voltage - electrostatic potential is continuous at the heterointerface - displacement D = eps E is continuous at the heterointerface - no transport/recombination calculation; under finite bias the plotted quasi-Fermi level is only a visual guide """ from dataclasses import dataclass import argparse import math import numpy as np import matplotlib.pyplot as plt EPS0 = 8.8541878128e-12 Q = 1.602176634e-19 @dataclass class Layer: name: str doping_cm3: float EV_eV: float Eg_eV: float EF_eV: float thickness_nm: float eps_r: float @property def doping_m3(self): return self.doping_cm3 * 1e6 @property def thickness_m(self): return self.thickness_nm * 1e-9 @property def eps(self): return self.eps_r * EPS0 @property def EC_eV(self): return self.EV_eV + self.Eg_eV @dataclass class JunctionResult: regime: str Vbi: float V: float Vj: float Wn_m: float Wp_m: float D_interface_C_m2: float D_left_C_m2: float D_right_C_m2: float Qn_ion_C_m2: float Qp_ion_C_m2: float Vn: float Vp: float @property def sigma_left_contact(self): # surface free charge on the left metal, with outward normal into semiconductor return self.D_left_C_m2 @property def sigma_right_contact(self): # sign convention reported as semiconductor-side D at the right boundary; # metal surface charge has the opposite normal convention. return -self.D_right_C_m2 def _standard_partial(n, p, Vj): ND, NA = n.doping_m3, p.doping_m3 denom = 1/(n.eps*ND) + 1/(p.eps*NA) sigma = math.sqrt(2*Q*Vj/denom) Wn = sigma/(Q*ND) Wp = sigma/(Q*NA) return Wn, Wp, sigma def _solve_n_full_p_partial(n, p, Vj): """ n layer fully depleted; p side partial. Boundary condition E(Wp)=0 in neutral p bulk. Let W = Wp. Vj = q/eps_n * (NA*W*d_n - ND*d_n^2/2) + q*NA*W^2/(2 eps_p) """ ND, NA = n.doping_m3, p.doping_m3 dn = n.thickness_m a = Q*NA/(2*p.eps) b = Q*NA*dn/n.eps c = -Q*ND*dn*dn/(2*n.eps) - Vj disc = b*b - 4*a*c if disc < 0: raise RuntimeError("No real p-depletion-width solution.") roots = [(-b + math.sqrt(disc))/(2*a), (-b - math.sqrt(disc))/(2*a)] roots = [w for w in roots if w >= 0] if not roots: raise RuntimeError("No positive p-depletion-width solution.") Wp = max(roots) D0 = Q*NA*Wp Dleft = D0 - Q*ND*dn Dright = 0.0 Vn = D0*dn/n.eps - Q*ND*dn*dn/(2*n.eps) Vp = Q*NA*Wp*Wp/(2*p.eps) return dn, Wp, D0, Dleft, Dright, Vn, Vp def _solve_n_partial_p_full(n, p, Vj): """ n side partial; p layer fully depleted. Boundary condition E(-Wn)=0 in neutral n bulk. """ ND, NA = n.doping_m3, p.doping_m3 dp = p.thickness_m a = Q*ND/(2*n.eps) b = Q*ND*dp/p.eps c = -Q*NA*dp*dp/(2*p.eps) - Vj disc = b*b - 4*a*c if disc < 0: raise RuntimeError("No real n-depletion-width solution.") roots = [(-b + math.sqrt(disc))/(2*a), (-b - math.sqrt(disc))/(2*a)] roots = [w for w in roots if w >= 0] if not roots: raise RuntimeError("No positive n-depletion-width solution.") Wn = max(roots) D0 = Q*ND*Wn Dleft = 0.0 Dright = D0 - Q*NA*dp Vn = Q*ND*Wn*Wn/(2*n.eps) Vp = D0*dp/p.eps - Q*NA*dp*dp/(2*p.eps) return Wn, dp, D0, Dleft, Dright, Vn, Vp def _solve_both_full(n, p, Vj): """ Both finite layers fully depleted. D0 = displacement at heterointerface. n side: D(x) = D0 + q ND (x-d_n), 0 <= x <= d_n p side, y=x-d_n: D(y) = D0 - q NA y, 0 <= y <= d_p Vj = integral(E dx) over the whole semiconductor. """ ND, NA = n.doping_m3, p.doping_m3 dn, dp = n.thickness_m, p.thickness_m cap_geom = dn/n.eps + dp/p.eps fixed = Q/2 * (ND*dn*dn/n.eps + NA*dp*dp/p.eps) D0 = (Vj + fixed) / cap_geom Dleft = D0 - Q*ND*dn Dright = D0 - Q*NA*dp Vn = D0*dn/n.eps - Q*ND*dn*dn/(2*n.eps) Vp = D0*dp/p.eps - Q*NA*dp*dp/(2*p.eps) return dn, dp, D0, Dleft, Dright, Vn, Vp def solve_depletion(n: Layer, p: Layer, V: float): ND, NA = n.doping_m3, p.doping_m3 dn, dp = n.thickness_m, p.thickness_m if min(ND, NA, dn, dp, n.eps_r, p.eps_r) <= 0: raise ValueError("Dopings, thicknesses and dielectric constants must be positive.") Vbi = n.EF_eV - p.EF_eV if Vbi <= 0: raise ValueError( f"Computed Vbi={Vbi:.6g} V <= 0. Check the common energy reference." ) Vj = Vbi - V if Vj < 0: raise ValueError( f"Vbi - V = {Vj:.6g} V < 0. " "Strong forward-bias transport is outside this depletion model." ) if Vj == 0: return JunctionResult( "flat-band", Vbi, V, Vj, 0, 0, 0, 0, 0, 0, 0, 0, 0 ) Wn0, Wp0, sigma0 = _standard_partial(n, p, Vj) tol = 1e-12 if Wn0 <= dn*(1+tol) and Wp0 <= dp*(1+tol): Wn = min(Wn0, dn) Wp = min(Wp0, dp) D0 = Q*ND*Wn Vn = Q*ND*Wn*Wn/(2*n.eps) Vp = Q*NA*Wp*Wp/(2*p.eps) return JunctionResult( "partial / partial", Vbi, V, Vj, Wn, Wp, D0, 0.0, 0.0, Q*ND*Wn, -Q*NA*Wp, Vn, Vp ) # The side with smaller maximum ionized sheet charge reaches full depletion first. Qnmax = Q*ND*dn Qpmax = Q*NA*dp if Qnmax <= Qpmax: vals = _solve_n_full_p_partial(n, p, Vj) Wn, Wp, D0, Dl, Dr, Vn, Vp = vals if Wp <= dp*(1+tol): Wp = min(Wp, dp) return JunctionResult( "n full / p partial", Vbi, V, Vj, Wn, Wp, D0, Dl, Dr, Q*ND*dn, -Q*NA*Wp, Vn, Vp ) vals = _solve_both_full(n, p, Vj) else: vals = _solve_n_partial_p_full(n, p, Vj) Wn, Wp, D0, Dl, Dr, Vn, Vp = vals if Wn <= dn*(1+tol): Wn = min(Wn, dn) return JunctionResult( "n partial / p full", Vbi, V, Vj, Wn, Wp, D0, Dl, Dr, Q*ND*Wn, -Q*NA*dp, Vn, Vp ) vals = _solve_both_full(n, p, Vj) Wn, Wp, D0, Dl, Dr, Vn, Vp = vals return JunctionResult( "n full / p full", Vbi, V, Vj, Wn, Wp, D0, Dl, Dr, Q*ND*dn, -Q*NA*dp, Vn, Vp ) def make_profiles(n: Layer, p: Layer, r: JunctionResult, npts=3000): dn, dp = n.thickness_m, p.thickness_m xi = dn x = np.linspace(0, dn+dp, npts) D = np.zeros_like(x) # n side nmask = x <= xi xn = dn - r.Wn_m if r.Wn_m < dn: neutral_n = x < xn dep_n = (x >= xn) & (x <= xi) D[neutral_n] = 0.0 D[dep_n] = Q*n.doping_m3*(x[dep_n] - xn) else: D[nmask] = r.D_interface_C_m2 + Q*n.doping_m3*(x[nmask] - dn) # p side y = x - dn xp = r.Wp_m pmask = x >= xi if r.Wp_m < dp: dep_p = (y >= 0) & (y <= xp) neutral_p = y > xp D[dep_p] = r.D_interface_C_m2 - Q*p.doping_m3*y[dep_p] D[neutral_p] = 0.0 else: D[pmask] = r.D_interface_C_m2 - Q*p.doping_m3*y[pmask] E = np.where(nmask, D/n.eps, D/p.eps) # phi(0)=0 and E=-dphi/dx phi = np.zeros_like(x) dx = np.diff(x) phi[1:] = -np.cumsum(0.5*(E[:-1]+E[1:])*dx) EV0 = np.where(nmask, n.EV_eV, p.EV_eV) EC0 = np.where(nmask, n.EC_eV, p.EC_eV) EV = EV0 - phi EC = EC0 - phi # Quasi-Fermi guide only. At V=0 it is exactly flat. EFn = n.EF_eV EFp = EFn - r.V EFguide = np.linspace(EFn, EFp, len(x)) if abs(r.V) < 1e-15: EFguide[:] = EFn return dict( x_m=x, D_C_m2=D, E_V_m=E, phi_V=phi, EV_eV=EV, EC_eV=EC, EFguide_eV=EFguide, interface_m=xi, xn_m=max(0.0, xn), xp_m=min(dn+dp, dn+xp) ) def print_result(n, p, r): print("\n=== Finite-thickness hetero-pn depletion result ===") print(f"Regime = {r.regime}") print(f"Applied bias V = {r.V: .6g} V (positive = forward)") print(f"Built-in voltage Vbi = {r.Vbi: .6g} V") print(f"Junction drop Vbi - V = {r.Vj: .6g} V") print() print(f"n depletion width Wn = {r.Wn_m*1e9: .6g} nm") print(f"p depletion width Wp = {r.Wp_m*1e9: .6g} nm") print(f"n thickness dn = {n.thickness_nm: .6g} nm") print(f"p thickness dp = {p.thickness_nm: .6g} nm") print() print(f"n-side ionized charge = {r.Qn_ion_C_m2: .6g} C/m^2") print(f"p-side ionized charge = {r.Qp_ion_C_m2: .6g} C/m^2") print(f"D at interface = {r.D_interface_C_m2: .6g} C/m^2") print(f"D at n outer surface = {r.D_left_C_m2: .6g} C/m^2") print(f"D at p outer surface = {r.D_right_C_m2: .6g} C/m^2") print() print(f"Potential drop in n layer = {r.Vn: .6g} V") print(f"Potential drop in p layer = {r.Vp: .6g} V") print(f"check Vn + Vp = {r.Vn+r.Vp: .6g} V") if "full" in r.regime: print("\nFinite-thickness note:") print(" Once a layer is fully depleted, extra bias is supported by") print(" external-contact charge as well as ionized dopant charge.") print(" Hence depletion charges of n and p layers need not remain equal.") def plot_band_diagram(n, p, r, outfile=None, show=True, npts=3000, neg_phi=False): prof = make_profiles(n, p, r, npts) xnm = prof["x_m"]*1e9 fig, ax1 = plt.subplots(figsize=(10.8, 6.6)) ax1.plot(xnm, prof["EC_eV"], label=r"$E_C(x)$") ax1.plot(xnm, prof["EV_eV"], label=r"$E_V(x)$") ax1.plot(xnm, prof["EFguide_eV"], "--", label=r"$E_F$ / quasi-Fermi guide") ax1.axvline(n.thickness_nm, linestyle=":", linewidth=1.2, label="n/p interface") if r.Wn_m < n.thickness_m: ax1.axvline((n.thickness_m-r.Wn_m)*1e9, linestyle="--", linewidth=0.9) if r.Wp_m < p.thickness_m: ax1.axvline((n.thickness_m+r.Wp_m)*1e9, linestyle="--", linewidth=0.9) ax1.set_xlabel("Depth from n-layer surface (nm)") ax1.set_ylabel("Energy (eV)") ax1.grid(True, alpha=0.25) ax2 = ax1.twinx() if neg_phi: ax2.plot(xnm, -prof["phi_V"], "-.", label=r"$-\phi(x)$") ax2.set_ylabel(r"Negative electrostatic potential $-\phi$ (V)") else: ax2.plot(xnm, prof["phi_V"], "-.", label=r"$\phi(x)$") ax2.set_ylabel(r"Electrostatic potential $\phi$ (V)") deltaE1 = ax1.get_ylim() deltaE2 = ax2.get_ylim() range1 = deltaE1[1] - deltaE1[0] range2 = deltaE2[1] - deltaE2[0] center2 = (deltaE2[0] + deltaE2[1]) / 2 ax2.set_ylim( center2 - range1 / 2, center2 + range1 / 2 ) ax1.set_title( f"Hetero p-n: {r.regime}, V={r.V:g} V, Vbi={r.Vbi:.4g} V, " f"Wn={r.Wn_m*1e9:.3g} nm, Wp={r.Wp_m*1e9:.3g} nm" ) fig.tight_layout() if outfile: fig.savefig(outfile, dpi=180, bbox_inches="tight") print(f"\nSaved plot: {outfile}") if show: plt.show() else: plt.close(fig) def parser(): a = argparse.ArgumentParser() a.add_argument("--ND", type=float, required=True) a.add_argument("--n-EV", type=float, default=0.0) a.add_argument("--n-Eg", type=float, required=True) a.add_argument("--n-EF", type=float, required=True) a.add_argument("--n-d", type=float, required=True, help="nm") a.add_argument("--n-eps", type=float, required=True) a.add_argument("--NA", type=float, required=True) a.add_argument("--p-EV", type=float, required=True) a.add_argument("--p-Eg", type=float, required=True) a.add_argument("--p-EF", type=float, required=True) a.add_argument("--p-d", type=float, required=True, help="nm") a.add_argument("--p-eps", type=float, required=True) a.add_argument("-V", "--bias", type=float, default=0.0) a.add_argument("-o", "--output", default="hetero_pn_band_finite.png") a.add_argument("--no-show", action="store_true") a.add_argument("--npts", type=int, default=3000) a.add_argument("--neg-phi", type=int, default=1, help="Plot -phi(x) instead of electrostatic potential phi(x)" ) return a def main(): args = parser().parse_args() n = Layer("n", args.ND, args.n_EV, args.n_Eg, args.n_EF, args.n_d, args.n_eps) p = Layer("p", args.NA, args.p_EV, args.p_Eg, args.p_EF, args.p_d, args.p_eps) r = solve_depletion(n, p, args.bias) print_result(n, p, r) plot_band_diagram(n, p, r, args.output, not args.no_show, args.npts, args.neg_phi) if __name__ == "__main__": main()