Self-consistently solve coupled Poisson-Drift Diffusion eqns. for semiconductor device analysis

Question

I have a system of coupled non-linear ODEs.

$$J = \mu e \left( n(x) E(x) + \frac{K T}{e} \frac{dn(x)}{dx}\right)$$ $$\frac{dE(x)}{dx} = \frac{4 \pi e}{\epsilon}[N_D(x) -n(x)]$$

The first equation is drift-diffusion, the second is Gauss' Law. I am interested in solving this system self-consistently for the carrier profile, $n(x)$, and the electric field, $E(x)$, using finite differences. The current, $J$, the mobility, $\mu$, elementary charge, $e$, are known constants. The doping profile $N_D(x)$ is known, and I have boundary conditions for the carrier concentration $n(0) = N_D(0)$ and $n(L) = N_D(L)$. I am unsure how to proceed with solving these equations.

I am following the analysis of Baranger and Wilkins in the attached review. Appendix B, specifically the section dealing with guessing an initial density.

"Equation (Bl) is solved self-consistently with the Poisson equation [Eq. (2)] for n (x) and E(x) using a finite difference method. As boundary conditions we fix the current JDD at the desired value and set the density equal to N+ at the edges of the structure."

EDIT: I have taken the advise of @hypernova and implemented his discretization scheme. In order to implement, I have non-dimensionalized the equations to prevent problems with numerical implementation.

Non-dimensionalizing the parameters: $$N_D^* = \frac{N_D(x)}{N_c}$$ $$n^* = \frac{n(x)}{N_c}$$ $$x^* = \frac{x}{L}$$ $$E^* = \frac{E(x) e L}{K T}$$

Applying the scaling: $$\theta = \frac{J L}{\mu N_c K T}$$ $$\gamma = \frac{K T \epsilon}{4 \pi e^2 L^2 N_c}$$

We obtain the non-dimensional equations: $$\frac{dn^*}{dx^*} = \theta - n^* E^*$$ $$\frac{dE^*}{dx^*} = \frac{1}{\gamma} [N_D^* - n^*]$$

We make the initial guess $n^{*,(0)} = N_D^*(0)$.

Dicretize the problem as per the answer below:

$$\frac{E^{*,(k,*)}_{j+1} - E^{*,(k,*)}_{j}}{h^*} = \frac{1}{2 \gamma} [N_{D,j+1}^* + N_{D,j}^* - n_{j+1}^{*,(k)} -n_{j}^{*,(k)}]$$ for $j = 0,1,...,N-1$.

We make the correction:

$$C^{(k)} = \frac{N_D^*{(1)} - N_D^*{(0)} - A^{(k)}}{B^{(k)}}$$ $$A^{(k)} = h^* \{\frac{1}{2} [ 2 \theta - n_0^*{(k)} E_0^{*,(k,*)} - n_N^*{(k)} E_N^{*,(k,*)}] + \sum_{j=1}^{N-1} (\theta - n_j^{*,(k)} E_j^{*,(k,*)}) \}$$ $$B^{(k)} = h^* \{ \frac{1}{2} [-n_0^{*,(k)}-n_N^{*,(k)}] + \sum_{j=1}^{N-1} (-n_j^{*,(k)}) \}$$

$$E_j^{*,(k)} = E_j^{*,(k,*)} + C^{(k)}$$

And finally update the density: $$\frac{n^{*,(k+1)}_{j+1} - n^{*,(k+1)}_{j}}{h^*} = \frac{1}{2} [2 \theta - n^{*,(k)}_{j+1} E^{*,(k)}_{j+1} - n^{*,(k)}_{j} E^{*,(k)}_{j}]$$ for $j = 0,1,...,N-1$. With condition $n^{*,(k+1)}_{0} = N_D^*(0)$.

I believe that this discretized scheme is correct. However, I have had trouble implementing. Namely, after only two iterations, the density, correction factors, and E-field quantities explode. Furthermore, another curious aspect is that the densities become strongly negative in the center of the domain, and the right boundary does not exactly match the imposed condition.

Some parameter values for implementation: $$J = 10^3 [A/cm^2]$$ $$K = 1.380*10^-23 [J/K]$$ $$T = 300 [K]$$ $$\mu = 7500 [cm^2 V^-1 s^-1]$$ $$L = 4.4*10^-6 [m]$$ $$e = 1.602*10^-19 [C]$$ $$\epsilon = 8.85*10^-12 [F/m]$$ $$N_C = 10^{15} [cm^-3]$$ $$N_{0,L} = 10^{18} [cm^-3]$$

The doping profile, $N_D$ is a step function:

$ N_D(x) \left\{ \begin{array}{ll} N_{0,L} & 0 \leq x\leq 2 \mu m \\ N_C & 2 < x < 2.4 \mu m \\ N_{0,L} & 2.4 \mu m\leq x \leq L \\ \end{array} \right. $

Answer 1

Restatement of equations

Define \begin{align} \alpha&=-\frac{4\pi e}{\epsilon},\\ \beta&=-\frac{1}{KT},\\ f&=-\alpha N_D,\\ g&=-\frac{\beta J}{\mu e}, \end{align} and the system reads equivalently \begin{align} \frac{{\rm d}E}{{\rm d}x}&=\alpha n+f,&&x\in\left(0,L\right),\\ \frac{{\rm d}n}{{\rm d}x}&=\beta nE+g,&&x\in\left(0,L\right),\\ n&=N_D(0),&&x=0,\\ n&=N_D(L),&&x=L. \end{align}

Iterative scheme: Summary

Since the system is nonlinear, one must use either a nonlinear solver or an iterative scheme. Based on my experience, a Picard-based iteration (see here and here for more information) would suffice to provide accurate results with efficiency.

The idea goes as follows.

Initially, guess an expression for $n(x)$. This guess could be completely arbitrary, and would not have impact on the computational result. Typically, a linear function that satisfies the boundary conditions appears to be an ideal candidate, i.e., $$ n^{(0)}(x)=N_D(0)+\frac{N_D(L)-N_D(0)}{L}x. $$

Thanks to the initial guess, \begin{align} \frac{{\rm d}E}{{\rm d}x}&=\alpha n^{(0)}+f,&&x\in\left(0,L\right) \end{align} requires no more than an indefinite integration. Artificially fake $E(0)=0$, and the integration would give you $E^{(0,*)}$. Of course, this $E^{(0,*)}$ is not satisfactory at all, because $E(0)=0$ may very well be a wrong implementation. We need to fix this bug and correct $E^{(0,*)}$ to some $E^{(0)}$, making it consistent with the boundary conditions for $n$.

Note that $$ N_D(L)-N_D(0)=n(L)-n(0)=\int_0^L\frac{{\rm d}n}{{\rm d}x}{\rm d}x=\int_0^L\left(\beta nE+g\right){\rm d}x. $$ Therefore, it is natural to require our expected $E^{(0)}$ to satisfy $$ N_D(L)-N_D(0)=\int_0^L\left(\beta n^{(0)}E^{(0)}+g\right){\rm d}x. $$ With this in mind, simply put $$ E^{(0)}=E^{(0,*)}+C, $$ where $$ C=\frac{N_D(L)-N_D(0)-\int_0^L\left(n^{(0)}E^{(0,*)}+g\right){\rm d}x}{\int_0^L\beta n^{(0)}{\rm d}x}. $$ You may check that this $E^{(0)}$ does meet our need.

Now as you have $n^{(0)}$ and $E^{(0)}$, \begin{align} \frac{{\rm d}n}{{\rm d}x}&=\beta n^{(0)}E^{(0)}{{\rm d}x}+g,&&x\in\left(0,L\right),\\ n&=N_D(0),&&x=0 \end{align} is another indefinite integration, whose solution reads $n^{(1)}$. Besides, thanks to our correction of $E^{(0,*)}$ to $E^{(0)}$, the other boundary condition $n(L)=N_D(L)$ will automatically be satisfied. Thus this $n^{(1)}$ does not need any further correction, and can be made use of directly.

So you have $n^{(1)}$ by hand, you may solve the first system (together with correction, of course) and get $E^{(1)}$. And as you have $n^{(1)}$ and $E^{(1)}$ by hand, you may solve the second system and get $n^{(2)}$... You keep running this iterative process, until, say, $$ \max_{x\in\left(0,L\right)}\left\{\left|n^{(k)}-n^{(k-1)}\right|+\left|E^{(k)}-E^{(k-1)}\right|\right\} $$ is smaller than some tolerance (e.g., $10^{-8}$ or as per your demand for accuracy) for some $k$.

Iterative scheme: Discretization

Let $0=x_0<x_1<x_2<\cdots<x_N=L$ be the equi-spaced grid points for $\left[0,L\right]$, with $x_j=jh$ and $h=L/N$. Let $n_j^{(k)}$ be an approximation of $n^{(k)}(x_j)$, and $E_j^{(k)}$ be an approximation of $E^{(k)}(x_j)$.

Then the first system \begin{align} \frac{{\rm d}E^{(k,*)}}{{\rm d}x}&=\alpha n^{(k)}+f,&&x\in\left(0,L\right),\\ E^{(k,*)}&=0,&&x=0 \end{align} observes the following discretization \begin{align} \frac{E_{j+1}^{(k,*)}-E_j^{(k,*)}}{h}&=\frac{1}{2}\left[\left(\alpha n_{j+1}^{(k)}+f_{j+1}\right)+\left(\alpha n_j^{(k)}+f_j\right)\right],&&j=0,1,\cdots,N-1,\\ E_0^{(k,*)}&=0. \end{align}

The correction goes that \begin{align} E_j^{(k)}&=E_j^{(k,*)}+C^{(k)},&&j=0,1,\cdots,N, \end{align} where $$ C^{(k)}=\frac{N_D(L)-N_D(0)-A^{(k)}}{B^{(k)}}, $$ with $$ A^{(k)}=\left\{\frac{1}{2}\left[\left(\beta n_0^{(k)}E_0^{(k,*)}+g_0\right)+\left(\beta n_N^{(k)}E_N^{(k,*)}+g_N\right)\right]+\sum_{j=1}^{N-1}\left(\beta n_j^{(k)}E_j^{(k,*)}+g_j\right)\right\}h $$ and $$ B^{(k)}=\left[\frac{1}{2}\left(\beta n_0^{(k)}+\beta n_N^{(k)}\right)+\sum_{j=1}^{N-1}\beta n_j^{(k)}\right]h. $$

The second system \begin{align} \frac{{\rm d}n^{(k+1)}}{{\rm d}x}&=\beta n^{(k)}E^{(k)}+g,&&x\in\left(0,L\right),\\ n&=N_D(0),&&x=0 \end{align} observes the following discretization \begin{align} \frac{n_{j+1}^{(k+1)}-n_j^{(k+1)}}{h}&=\frac{1}{2}\left[\left(\beta n_{j+1}^{(k)}E_{j+1}^{(k)}+g_{j+1}\right)+\left(\beta n_j^{(k)}E_j^{(k)}+g_j\right)\right],&&j=0,1,\cdots,N-1,\\ n^{(k+1)}_0&=N_D(0). \end{align}

The relative error goes that $$ \varepsilon^{(k+1)}=\max_j\left\{\left|n_j^{(k+1)}-n_j^{(k)}\right|+\left|E_j^{(k+1)}-E_j^{(k)}\right|\right\}. $$

That is it! Hope this could be helpful for you.