# Solution of non-linear Poisson equation does not match reference

I'm trying to solve the non-linear Poisson equation as a first step to solve the drift-diffusion equations for semiconductors. For reference, I'm using a preprint from the Weierstrass Institut (which I will refer to as [W]) that can be found here. All the model parameters and equations can be found there. I'll make references to specific parts of this paper.

I started implementing the algorithm in 2.4 in [W]. The discretization is done using FVM and the solution is then calculated using the Newton method. I ran into two issues:

1. The local electroneutrality in my case is the negative of what they get in [W]. What could be the cause of that? Do I have to switch the sign?

2. The non-linearity (_assemble_nonlinear) does not contribute correctly to the solution. Comparing my solution to the one in [W], Figure 2, it is clear that something went wrong. I'm suspecting, that the problem lies in the non-linearity since without it, the solution is a homogenous solution to the Poisson equation. Here is my code with annotations refering to [W]:

import numpy as np
import sys
from typing import Callable
from scipy.constants import e as q, epsilon_0 as eps0
from scipy.constants import physical_constants
from scipy.linalg import solve_banded
from scipy.sparse import diags
import matplotlib.pyplot as plt
import matplotlib as mpl

np.set_printoptions(suppress=True,linewidth=np.nan,threshold=sys.maxsize)

kB = physical_constants["Boltzmann constant"]
kBev = physical_constants["Boltzmann constant in eV/K"]

# Setup for plotting
plt.style.use("ggplot")
plt.rcParams["text.usetex"] = True
plt.rcParams["font.family"] = "serif"
cmap = mpl.cm.get_cmap("summer")
plt.rcParams["axes.prop_cycle"] = mpl.cycler(
color=[cmap(v) for v in np.linspace(0, 0.7, 10)]
)

class PoissonSolver:
def __init__(self,
grid : np.ndarray,
C : Callable,
T : float,
Ev : float,
Ec : float,
Nv : float,
Nc : float,
epsr : float
):
"""
grid : np.ndarray
Grid for the calculation
C : Callable
Doping profile
T : float
Temperature in K
Ev : float
Valence band energy in eV
Ec : float
Conduction band energy in eV
Nv : float
Nc : float
epsr : float
Relative permettivity
"""
self.N = len(grid)  # Number of grid points
self.xks = grid
self.xkk1s = np.pad(1 / 2 * (self.xks[1:] + self.xks[:-1]), (1, 1), constant_values=(self.xks, self.xks[-1]))
self.hs = self.xks[1:] - self.xks[:-1]
self.omegas = self.xkk1s[1:] - self.xkk1s[:-1]
self.C = C(grid) # Doping profile
self.T = T  # Temperature in K
self.U_T =  kB * self.T / q  # Thermal voltage
self.Ev = q * Ev  # Valence Band energy in J
self.Ec = q * Ec  # Conduction Band energy in J
self.Nv = Nv
self.Nc = Nc
self.epss = epsr * eps0  # Permettivity
self.psis = []

def _calc_psi0(self):
# Intrisic carrier density squared from Eq. (9)
Nisq = self.Nc * self.Nv * np.exp(-(self.Ec - self.Ev) / (kB * T))
nvals = np.where(self.C <= 0)
pvals = np.where(self.C > 0)
# Numerically stable quadratic formula used in Eq. (8)
quad[nvals] = -2 * np.sqrt(Nisq) / (self.C[nvals] - np.sqrt(self.C[nvals] ** 2 + 4 * Nisq))
quad[pvals] = (self.C[pvals] + np.sqrt(self.C[pvals] ** 2 + 4 * Nisq)) / (2 * np.sqrt(Nisq))
# Actual Eq. (8)
psi0 = - (self.Ec - self.Ev) / (2 * q) - 1 / 2 * self.U_T * np.log(self.Nc / self.Nv) + self.U_T * np.log(quad)

# Plot potential
#fig = plt.figure()
#ax.plot(self.xks, psi0)
#ax.set_xlabel("Position [$$\mu m$$]")
#ax.set_ylabel("Initial Potential [V]")
#plt.show()

return psi0.reshape((self.N, 1))

def _assemble_matrix(self):
M = np.zeros((3, self.N))  # Only store diagonals
K = np.zeros((self.N, self.N))  # Store entire matrix
for i in range(self.N):
# M is Jacobian of Eq. (12)
# K is the matrix for the linear part of Eq. (12)
if i == 0:
M[1, i] = 1
M[2, i] = - self.epss / self.hs[i]
elif i == self.N - 1:
M[0, i] = - self.epss / self.hs[i - 1]
M[1, i] = 1
else:
# Derivatives of the mobile carrier concentration
En = q * self.psis[-1][i]
dpdpsi = -Nv / self.U_T * np.exp((self.Ev - En) / (kB * T))
dndpsi =  Nc / self.U_T * np.exp((En - self.Ec) / (kB * T))

M[1, i] = + self.epss * (1 / self.hs[i] + 1 / self.hs[i - 1]) + q * self.omegas[i] * (dpdpsi - dndpsi)
if i != 1:
M[0, i] = - self.epss / self.hs[i]
if i != self.N - 2:
M[2, i] = - self.epss / self.hs[i - 1]

K[i, i - 1] = - self.epss / self.hs[i - 1]
K[i, i] = + self.epss * (1 / self.hs[i] + 1 / self.hs[i - 1])
K[i, i + 1] = - self.epss / self.hs[i]

return M, K

def _assemble_nonlinear(self):
# Non-linear part of Eq.(12)
S = np.zeros((self.N, 1))
for i in range(1, self.N - 1):
# Mobile carrier concentration
En = q * self.psis[-1][i]
p = Nv * np.exp((self.Ev - En) / (kB * T))
n = Nc * np.exp((En - self.Ec) / (kB * T))

S[i] = -q * (self.C[i] + p - n) * self.omegas[i]
plt.show()

return S

def solve(self, relax=1.0):
self.psis = [self._calc_psi0()]
for n in range(9):
# Newton iterations
dFdpsi, K = self._assemble_matrix()
S = self._assemble_nonlinear()
F = np.dot(K, self.psis[-1]) + S
dpsi = solve_banded((1, 1), dFdpsi, -relax * F)
psi = self.psis[-1] + dpsi
self.psis.append(psi)

self.psis = np.array(self.psis)
for psi in self.psis:
plt.plot(self.xks, psi)
plt.show()

# Same parameters as in Fig. 2
T = 300
Ev = 0
Ec = 1.424
Nv = 9e18
Nc = 4.7e17
epsr = 12.9
NA = 1e17
ND = 1e16
grid = np.linspace(0, 10e-6, 51)

def C(x):
return -ND * np.heaviside(-(x - 1e-6), 1) + NA * np.heaviside((x - 9e-6), 1)

ps = PoissonSolver(grid, C=C, T=T, Ev=Ev, Ec=Ec, Nv=Nv, Nc=Nc, epsr=epsr)
#fig = plt.figure()
#ax.set_xlabel("Position [$$\mu m$$]")