# Jacobian of the electron and hole drift diffusion equations with respect to potential in semiconductors

These are the discretized drift diffusion equations as taken from the book "Analysis and Simulation of Semiconductor Devices". The electron continuity equation is:

$$((\Delta_{j-1}^y+\Delta_j^y)/2\Delta_i^x).D_{i+1/2,j}^n.B((\varphi_{i+1,j}-\varphi_{i,j})/V_T).n_{i+1,j} + ((\Delta_{j-1}^y+\Delta_j^y)/2\Delta_{i-1}^x).D_{i-1/2,j}^n.B((\varphi_{i-1,j}-\varphi_{i,j})/V_T).n_{i-1,j} - [(((\Delta_{j-1}^y+\Delta_j^y)/2\Delta_i^x).D_{i+1/2,j}^n.B((\varphi_{i,j}-\varphi_{i+1,j})/V_T)) + (((\Delta_{j-1}^y+\Delta_j^y)/2\Delta_{i-1}).D_{i-1/2,j}^n.B((\varphi_{i,j}-\varphi_{i-1,j})/V_T)) + (((\Delta_{i-1}^x+\Delta_i^x)/2\Delta_j^y).D_{i,j+1/2}^n.B((\varphi_{i,j}-\varphi_{i,j+1})/V_T)) + (((\Delta_{i-1}^x+\Delta_i^x)/2\Delta_{j-1}^y).D_{i,j-1/2}^n.B((\varphi_{i,j}-\varphi_{i,j-1})/V_T))].n_{i,j} + ((\Delta_{i-1}^x+\Delta_i^x)/2\Delta_j^y).D_{i,j+1/2}^n.B((\varphi_{i,j+1}-\varphi_{i,j})/V_T).n_{i,j+1} + ((\Delta_{i-1}^x+\Delta_i^x)/2\Delta_{j-1}^y).D_{i,j-1/2}^n.B((\varphi_{i,j-1}-\varphi_{i,j})/V_T).n_{i,j-1} = R_{i,j}^{EFFECTIVE}.((\Delta_{i-1}^x+\Delta_i^x)/2).((\Delta_{j-1}^y+\Delta_j^y)/2)$$

The hole continuity equation is:

$$((\Delta_{j-1}^y+\Delta_j^y)/2\Delta_i^x).D_{i+1/2,j}^p.B((\varphi_{i,j}-\varphi_{i+1,j})/V_T).p_{i+1,j} + ((\Delta_{j-1}^y+\Delta_j^y)/2\Delta_{i-1}^x).D_{i-1/2,j}^p.B((\varphi_{i,j}-\varphi_{i-1,j})/V_T).p_{i-1,j} - [(((\Delta_{j-1}^y+\Delta_j^y)/2\Delta_i^x).D_{i+1/2,j}^p.B((\varphi_{i+1,j}-\varphi_{i,j})/V_T)) + (((\Delta_{j-1}^y+\Delta_j^y)/2\Delta_{i-1}^x).D_{i-1/2,j}^p.B((\varphi_{i-1,j}-\varphi_{i,j})/V_T)) + (((\Delta_{i-1}^x+\Delta_i^x)/2\Delta_j^y).D_{i,j+1/2}^p.B((\varphi_{i,j+1}-\varphi_{i,j})/V_T)) + (((\Delta_{i-1}^x+\Delta_i^x)/2\Delta_{j-1}^y).D_{i,j-1/2}^p.B((\varphi_{i,j-1}-\varphi_{i,j})/V_T))].p_{i,j} + \\ ((\Delta_{i-1}^x+\Delta_i^x)/2\Delta_j^y).D_{i,j+1/2}^p.B((\varphi_{i,j}-\varphi_{i,j+1})/V_T).p_{i,j+1} + ((\Delta_{i-1}^x+\Delta_i^x)/2\Delta_{j-1}^y).D_{i,j-1/2}^p.B((\varphi_{i,j}-\varphi_{i,j-1})/V_T).p_{i,j-1} = R_{i,j}^{EFFECTIVE}.((\Delta_{i-1}^x+\Delta_i^x)/2).((\Delta_{j-1}^y+\Delta_j^y)/2\ ,$$

where $\Delta_i$ and $\Delta_j$ are the mesh spacing for finite difference grid along X and Y directions. $D_n$,$D_p$ are electron and hole diffusion coefficients.$B(x) = x/(e^x-1)$ is Bernoulli's function and $V_T$ is the thermal voltage.

Also, $$R^{EFFECTIVE} = R^{SRH} + R^{AU} + R^{OPT} + R^{II}$$ where $R^{EFFECTIVE}$ is the effective generation/recombination term. $R^{SRH}$, $R^{AU}$, $R^{OPT}$ and $R^{II}$ are SRH, Auger, Radiative and Impact Ionization Recombination terms.

Now, I am applying the coupled Newton's iterative scheme to solve the above equations where I need to differentiate the equations. When differentiating with respect to potential '$\varphi$', I noticed that $\varphi$ appears in three places namely,

1. High field Electron/Hole Mobility, $\mu_{n,p}^{high}=\mu_{n,p}^{low}⁄ [1 + (\mu_{n,p}^{LICN}.E_{n,p}/v_{n,p}^{sat})^{\beta_{n,p}}]^{1/\beta_{n,p}}$ and $D_{n,p} = \mu_{n,p}.V_T$;

2. Bernoulli's function, $B(x)$ where 'x' here means difference between potentials at neighboring grid points and

3. Impact Ionization $R^{II} = -\alpha_n.(|J_{i,j}^n|/q) - \alpha_p.(|J_{i,j}^p|/q)$, $\alpha_n = a_n.e^{-b_n/|E_{i,j}|}$, $\alpha_p = a_p.e^{-b_p/|E_{i,j}|}$ where $J^{n,p}$ and $E$ represent electron/hole currents and electric field respectively. Current terms are made of potential terms as in, $J_{i+1/2,j}^{nx}=D_{i+1/2,j}^n.[B(\frac{\varphi_{i,j}-\varphi_{i+1,j}}{V_T}).n_{i,j} - B(\frac{\varphi_{i+1,j}-\varphi_{i,j}}{V_T}).n_{i+1,j}]/\Delta_i^x \rightarrow (31)$, $J_{i,j+1/2}^{ny}=D_{i,j+1/2}^n.[B(\frac{\varphi_{i,j}-\varphi_{i,j+1}}{V_T}).n_{i,j} - B(\frac{\varphi_{i,j+1}-\varphi_{i,j}}{V_T}).n_{i,j+1}]/\Delta_j^y \rightarrow (32)$, $J_{i+1/2,j}^{px}=D_{i+1/2,j}^p.[B(\frac{\varphi_{i,j}-\varphi_{i+1,j}}{V_T}).p_{i+1,j} - B(\frac{\varphi_{i+1,j}-\varphi_{i,j}}{V_T}).p_{i,j}]/\Delta_i^x \rightarrow (33)$, $J_{i,j+1/2}^{py}=D_{i,j+1/2}^p.[B(\frac{\varphi_{i,j}-\varphi_{i+1,j}}{V_T}).p_{i,j+1} - B(\frac{\varphi_{i+1,j}-\varphi_{i,j}}{V_T}).p_{i,j}]/\Delta_j^y \rightarrow (34)$

So, basically I need to differentiate all of these to find the Jacobian with respect to $\varphi$. But, the problem arises when the difference between potentials become zero. This means that the entire electron/hole equation becomes a constant when the potential difference is zero. Hence, the Jacobian is also zero.

This would mean that the final Jacobian matrix would have entire columns and rows made out of zeros thus making the matrix singular and it is not possible to solve a singular matrix.

Is this the correct approach to solving using the Newton's method. What should I do when the potential difference is zero?

## References

1. Selberherr, Siegfried. Analysis and simulation of semiconductor devices. Springer Science & Business Media, 2012.
• First, it is very hard to read such equations after discretization. Maybe you can add the original equations as well. Second, you should calculate the derivative of $\phi$ with respect to the PDE before discretization. This way you might see already the subtleties causing your Jacobian to fail and you will have no problems with discretization afterwards. If the Jacobian still fails, you need to think hard if it still describes proper physics within the assumptions of the model or how the model needs to behave in such a situation. – Bort Aug 14 '17 at 14:38
• Have you tried a Computer Algebra System (CAS) to compute the Jacobian? You can compute it and then automatically generate the code fromit. – nicoguaro Aug 14 '17 at 17:48
• I think you're just not computing the Jacobian correctly. You should take the derivative and then evaluate at zero potential. It sounds like you are evaluating and then taking the derivative, which doesn't make sense. – Dave Kielpinski Aug 14 '17 at 22:43
• @Dave Kielpinski - Taking the derivative is not possible for zero potential as because the Bernoulli's function is not valid at zero potential. At x=0, B(x) is taken to be 1 which makes all the equations constant with respect to potential thus making the differentiation 0. – P. Biswas Aug 17 '17 at 10:19
• Also, I am using Matlab to do all the compuations. – P. Biswas Aug 17 '17 at 10:20