The density gradient equations to be solved are of the following form:

\Lambda_e = -b_n\frac{\nabla^2 \sqrt{n}}{\sqrt{n}}

\Lambda_h = -b_p\frac{\nabla^2 \sqrt{p}}{\sqrt{p}}

where b_n, b_p are coefficients relating the gradient in \sqrt{n} with a quantum correction.

b_n = \frac{\gamma_n\hbar^2}{6 m_n}

b_p = \frac{\gamma_p\hbar^2}{6 m_p}

Using:

\frac{\nabla^2 \sqrt{n}}{\sqrt{n}} = \frac{1}{2} \left\{\nabla^2 \log n + \frac{1}{2} \left(\nabla \log n\right)^2\right\}

and

\Phi_e = \log n

\Phi_h = \log p

The equations become:

\Lambda_e + \frac{b_n}{2} \left( \nabla \cdot \nabla \Phi_e + \frac{1}{2} \left( \nabla \Phi_e \right)^2 \right) = 0

\Lambda_h + \frac{b_p}{2} \left( \nabla \cdot \nabla \Phi_h + \frac{1}{2} \left( \nabla \Phi_h \right)^2 \right) = 0

Considering the discretized equation along the edge between nodes i and j:

\Lambda_e \cdot \textrm{NV}_{i} + \sum_j \frac{b_{n_{i,j}}}{2}\left( \textrm{SA}_{i,j} \cdot \nabla_{i,j} \Phi_{e_{i,j}} - \frac{1}{2} \textrm{EV}_{i,j} \left( \nabla_{i,j} \Phi_e \right)^2 \right) - \textrm{INT}_{i,k} \frac{b_{n_{ox}}}{x_n} = 0

 \Lambda_h \cdot \textrm{NV}_{i} + \sum_j \frac{b_{p_{i,j}}}{2}\left( \textrm{SA}_{i,j} \cdot \nabla_{i,j} \Phi_{h_{i,j}} - \frac{1}{2} \textrm{EV}_{i,j} \left( \nabla_{i,j} \Phi_h \right)^2 \right) - \textrm{INT}_{i,k} \frac{b_{p_{ox}}}{x_p} = 0

The gradient is then:

\nabla_{i,j} \Phi_x = \frac{\Phi_{x_j} - \Phi_{x_i}}{L_{i_j}}

The symbols have the following meaning:

Note that the sign on the \textrm{EV} term depends on how the edge volume assembly is performed.

Wettstein :cite:`WettsteinVLSI2002` solves the equations in both the insulator in semiconductor materials. For now, we consider the approach of other authors and assume that the carriers diminish quickly in the insulator :cite:`GarciaAsenov2011`. Solving the equations in the insulator would be important for resonant tunneling between adjacent semiconductor regions :cite:`Hohr:sispad:2002`.

The calculations in the oxide are :cite:`jin2004`:

b_{n_{ox}} = \frac{\gamma_{n_{ox}} \hbar^2}{m_{ox}}

x_n = \frac{\hbar}{\sqrt{2 m_{ox} \Phi_B}}

At ohmic contacts the following boundary conditions are required to meet the equilibrium boundary condition assumption:

\begin{eqnarray}
\Lambda_e =& 0\\
\Lambda_h =& 0
\end{eqnarray}

The values of \Phi_e and \Phi_h are from the calculation of the electron and hole density from their respective Fermi levels.

\begin{eqnarray}
\Phi_e =& \frac{k T \log N_C + E_{F_n} - E_{C} - \Lambda_e}{k T}\\
\Phi_h =& \frac{k T \log N_V - E_{F_p} + E_{V} - \Lambda_h}{k T}
\end{eqnarray}

The convention chosen in this description is the \Lambda_e and \Lambda_h act to reduce the electron and hole concentration.

\begin{eqnarray}
n \propto \exp \left( -\Lambda_e \right)\\
p \propto \exp \left( -\Lambda_h \right)
\end{eqnarray}

The calculated values of \Lambda_e and \Lambda_h modify the driving force for current, as in other quantum correction models.

In :cite:`Wettstein:sispad:2002`, the authors point out that n_i for recombination must be scaled, to prevent large recombination near the interface. This is since:

R \propto n p - n_i^2

Therefore:

n_i^2 \propto \exp \left( \frac{-\Lambda_e - \Lambda_h}{k T} \right)

In prototyping the DG equations, 2 points of numerical overflow where discovered. When calculating the edge volume contribution, a heuristic like:

\left(\nabla_{i,j} \Phi_x\right)^2 = \left(b_x 10^8\right) \left(10^{-4}\left(\frac{\Phi_{x_j} - \Phi_{x_i}}{L_{i_j}}\right)\right)^2

was used. In addition, it is important to make sure that updates in \Lambda_e and \Lambda_h are not too large to cause an overflow in the calculation of n and p.

\begin{eqnarray}
n \propto \textrm{limexp}\left(\Phi_e\right) \\
p \propto \textrm{limexp}\left(\Phi_h\right)
\end{eqnarray}

where limiting the \exp function is necessary to prevent overflow.

Since the classical and density gradient solutions are different, it is necessary to ramp the parameters \gamma_n, and \gamma_p to improve convergence. A ramping strategy should be considered where the step change in \gamma_n and \gamma_p may be adjusted when a simulation fails to converge.

To ensure accurate simulation results, it may be necessary to apply mesh refinements away from the interface, where the maximum in carrier concentrations occur.

For a MOS device, it makes sense to only solve the DG equation for the carrier in the channel of the device.

\begin{eqnarray}
\frac{\nabla^2 \sqrt{n}}{\sqrt{n}} =& u \nabla \cdot \vec{v} \\
S_n =& \sqrt{n}\\
u   =& \frac{1}{S_n}\\
\nabla u =& -\frac{\nabla S_n}{S_n^2}\\
\vec{v} =& \nabla S_n\\
\nabla \cdot \vec{v} =& \nabla \cdot \nabla S_n\\
\nabla \cdot \left(u \vec{v}\right) =& u\nabla \cdot \vec{v} + \vec{v} \cdot \nabla u\\
u\nabla \cdot \vec{v} =& \nabla \cdot \left(u \vec{v}\right) - \vec{v} \cdot \nabla u\\
u \vec{v} =& \nabla \log S_n = \frac{1}{2} \nabla \log n\\
\vec{v} \cdot \nabla u =& -\nabla S_n \cdot \frac{\nabla S_n}{S_n^2} = -\left(\nabla \log S_n\right)^2 = -\frac{1}{4}\left(\nabla \log n\right)^2\\
\frac{\nabla^2 \sqrt{n}}{\sqrt{n}} =& \frac{1}{2} \nabla^2 \log n + \frac{1}{4} \left(\nabla \log n\right)^2
\end{eqnarray}

We start with deriving the equations in :cite:`WettsteinVLSI2002`.

\begin{eqnarray}
\Lambda_e = -b_n\frac{\nabla^2 \sqrt{n}}{\sqrt{n}}\\
\Lambda_h = -b_p\frac{\nabla^2 \sqrt{p}}{\sqrt{p}}
\end{eqnarray}

where

\begin{eqnarray}
b_n = \frac{\gamma_n\hbar^2}{6 m_n}\\
b_p = \frac{\gamma_p\hbar^2}{6 m_p}
\end{eqnarray}

The effect is such that:

\begin{eqnarray}
n = N_{C} \exp\left(\frac{E_F - E_C - \Lambda_e }{k T}\right)\\
p = N_{V} \exp\left(\frac{E_V - E_F - \Lambda_h }{k T}\right)
\end{eqnarray}

and the intrinsic carrier density is now:

n_i^2 \propto \exp\left(\frac{-\Lambda_e - \Lambda_h}{k T} \right)

This is especially important for recombination.

For convenience we define:

\begin{eqnarray}
n = \exp\left(\frac{\Phi_e }{k T}\right)\\
p = \exp\left(\frac{\Phi_h }{k T}\right)
\end{eqnarray}

where

\begin{eqnarray}
\Phi_e = E_F - E_C - \Phi_C - \Lambda_e \\
\Phi_h = E_V - E_F + \Phi_V - \Lambda_h \\
\end{eqnarray}

where

\begin{eqnarray}
\Phi_C =& -k T \log(N_C)\\
\Phi_V =& k T \log(N_V)
\end{eqnarray}

For current conduction, the effect is that:

\begin{eqnarray}
J_n =& n \mu \nabla \left(\Phi_e\right)  + q D_n \nabla n \\
J_p =& p \mu \nabla \left(\Phi_h\right)  - q D_p \nabla p
\end{eqnarray}

where

\begin{eqnarray}
\Phi_C =& -k T \log(N_C)\\
\Phi_V =& k T \log(N_V)
\end{eqnarray}

In the derivation which follows, they exploit the following relation:

\frac{\nabla^2 \sqrt{n}}{\sqrt{n}} = \frac{1}{2} \left\{\nabla^2 \log n + \frac{1}{2} \left(\nabla \log n\right)^2\right\}

For a volume integration:

\begin{eqnarray}
%\int \Lambda_e \partial v = -\frac{b_n}{2} \int \left\{\nabla^2 \log n + \frac{1}{2} \left(\nabla \log n\right)^2\right\} \partial v\\
%\int \Lambda_e \partial v = -\frac{b_n}{2} \int  \left\{\nabla \cdot \nabla \log n + \frac{1}{2} \left(\nabla \log n\right)^2\right\} \partial v \\
\int \Lambda_e \partial v = - \frac{b_n}{2} \left\{ \int \nabla \log n \cdot \partial s + \frac{1}{2} \int \left(\nabla \log n\right)^2 \partial v \right\}\\
\int \Lambda_h \partial v = - \frac{b_p}{2} \left\{ \int \nabla \log p \cdot \partial s + \frac{1}{2} \int \left(\nabla \log p\right)^2 \partial v \right\}
\end{eqnarray}

When assembled onto node i with respect to nodes j.

\begin{eqnarray}
\Lambda_{e,i} \Omega_i = \sum_j \frac{b_n \sigma_{i,j}}{2 l_{i,j}} \left\{ \left(\frac{\Phi_{e,j} -\Phi_{e,i}}{k T} \right) - \frac{1}{4} \left(\frac{\Phi_{e,j} -\Phi_{e,i}}{k T} \right)^2 \right\}\\
\Lambda_{h,i} \Omega_i = \sum_j \frac{b_p \sigma_{i,j}}{2 l_{i,j}} \left\{ \left(\frac{\Phi_{h,j} -\Phi_{h,i}}{k T} \right) - \frac{1}{4} \left(\frac{\Phi_{h,j} -\Phi_{h,i}}{k T} \right)^2 \right\}
\end{eqnarray}

Following :cite:`WettsteinDissertation` they discretize:

\begin{eqnarray}
%\int \Lambda_e \partial v = -b_n \int \frac{\nabla^2 \sqrt{n}}{\sqrt{n}} \partial v\\
%\int \Lambda_e \partial v = -b_n \int \frac{\nabla \cdot \nabla \sqrt{n}}{\sqrt{n}} \partial v\\
\int \Lambda_e \partial v = -b_n \int \frac{\nabla \sqrt{n}}{\sqrt{n}} \cdot \partial s\\
\int \Lambda_h \partial v = -b_p \int \frac{\nabla \sqrt{p}}{\sqrt{p}} \cdot \partial s
\end{eqnarray}

When assembled onto node i with respect to nodes j.

\begin{eqnarray}
\Lambda_i \Omega_i = \sum_j \frac{b_n \sigma_{i,j}}{l_{i,j}} \left(\frac{\sqrt{n_i} - \sqrt{n_j}}{\sqrt{n_i}} \right)\\
\Lambda_i \Omega_i = \sum_j \frac{b_n \sigma_{i,j}}{l_{i,j}} \left(1 - \frac{\sqrt{n_j}}{\sqrt{n_i}} \right)
\end{eqnarray}

Which then leads to:

\begin{eqnarray}
\Lambda_{e,i} \Omega_i = \sum_j \frac{b_n \sigma_{i,j}}{l_{i,j}} \left\{ 1 - \exp\left(\frac{\Phi_{e,i}}{2 k T} - \frac{\Phi_{e,j}}{2 k T} \right) \right\}\\
\Lambda_{h,i} \Omega_i = \sum_j \frac{b_n \sigma_{i,j}}{l_{i,j}} \left\{ 1 - \exp\left(\frac{\Phi_{h,i}}{2 k T} - \frac{\Phi_{h,j}}{2 k T} \right) \right\}
\end{eqnarray}