1
$\begingroup$

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. My solution

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"][0]
kBev = physical_constants["Boltzmann constant in eV/K"][0]

# 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[0], 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))
        quad = np.zeros_like(self.C)
        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 = fig.add_subplot(111)
        #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][0]
                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][0]
            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 = fig.add_subplot(111)
#ax.plot(ps.xks, ps.C)
#ax.set_xlabel("Position [$\mu m$]")
#ax.set_ylabel("Doping concentration [idk]")
#plt.show()
ps.solve(relax=1.0)

Thank you in advance for your help!

$\endgroup$
2
  • $\begingroup$ So there is a bug in your code somewhere. Are you asking us to debug the code for you? $\endgroup$ Sep 24 '20 at 16:20
  • $\begingroup$ @WolfgangBangerth I'd be grateful for a little code review in terms of physics and numerics. I don't necessarily think there is a bug in the code, I've been debugging it for multiple hours now (accumulated over days). Mostly, I'm unsure about how they obtained a positive electronegativity and what I'm doing wrong in terms of the non-linearity. $\endgroup$
    – aaron
    Sep 24 '20 at 17:02

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.