Charge densities

Classical

Things are much more simpler.

When any kind of Schrödinger equation is not solved, the charge carrier densities are estimated from the position-dependent conduction and valence band edges $E_{\text{c}}^i(\mathbf{x})$ and $E_{\text{v}}^i(\mathbf{x})$, quasi-Fermi levels, and the electrostatic potential $\varphi (\mathbf{x})$ in the context of Thomas-Fermi approximation.

These classical charge carrier densities are calculated as

(2.2.3)$$n(\mathbf{x})=\sum _{i\in \text{CB}}{N}_{\text{c}}^i(T){\mathcal{F}}_{1/2}([-E_{\text{c}}^i(\mathbf{x})+e\varphi (\mathbf{x})+E_{\text{F},n}(\mathbf{x})]/kT)$$

(2.2.4)$$p(\mathbf{x})=\sum _{i\in \text{VB}}{N}_{\text{v}}^i(T){\mathcal{F}}_{1/2}([E_{\text{v}}^i(\mathbf{x})-e\varphi (\mathbf{x})-E_{\text{F},p}(\mathbf{x})]/kT).$$

Here $N_{\text{v}}^i(T)$ and $N_{\text{v}}^i(T)$ are the equivalent density of states at the conduction and valence band edges, which are given by

(2.2.5)$$N_l^i(T)=g_l^i(\frac{m_{\text{dos,}\lambda }^{i}kT}{2\pi \hslash }{)}^{2/3}(l,\lambda )=(\text{v,h}),\text{or}(\text{c,e}).$$

Here $m_{\text{dos,}\lambda }^i$ is the density-of-mass for $d=3$ defined in (2.2.15).

This calculation of carrier densities is much faster than the quantm mechanical calculation, but the quantum effect such as energy quantization, carrier leackage into the barrier, etc. cannot be taken into account.

Also in this case, the carrier densities can be written as $n(\mathbf{x},\varphi )$ and $p(\mathbf{x},\varphi )$, which enters into the non-linear Poisson equation introduced next.

Moreover, when the current equation is included in the calculation scheme, seeing the carrier densities as $n(\mathbf{x},\varphi ,E_{\text{F},n})$ and $p(\mathbf{x},\varphi ,E_{\text{F},p})$ makes it easy to understand what the self-consistent calculation is actually doing.

Quantum mechanical

Multi-band model ($\mathbf{k}\cdot \mathbf{p}$ model)

Once the $\mu$-th component envelope function of the $j$-th eigenstate of electron ($l=\text{c}$) or hole ($l=\text{v}$) in the $i$-th band is obtained as $(F_{\mu }{)}_{l,j}^i(\mathbf{x})$ from the multi-band Schrödinger equation, the probability distribution of this $j$-th eigenstate reads

(2.2.6)$$p_{l,j}^i(\mathbf{x})=\sum _{\mu }|(F_{\mu }{)}_{l,j}^i(\mathbf{x}){|}^2.$$

where we are assuming 3D structure so far.

Then the quantum mechanical carrier densities for 3D structure are defined from these probability densities, energy eigenvalues $E_{\text{c},j}$ and $E_{\text{v},j}$, position-dependent quasi-Fermi levels $E_{\text{F},n}(\mathbf{x})$ and $E_{\text{F},p}(\mathbf{x})$ as

(2.2.7)$$n(\mathbf{x})=\sum _{i\in \text{CB}}{g}_{\text{c}}^i\sum _j{p}_{\text{c},j}^i(\mathbf{x})f([E_{\text{c},j}^i-E_{\text{F},n}(\mathbf{x})]/kT)$$

(2.2.8)$$p(\mathbf{x})=\sum _{i\in \text{VB}}{g}_{\text{v}}^i\sum _j{p}_{\text{v},j}^i(\mathbf{x})f([-E_{\text{v},j}^i+E_{\text{F},n}(\mathbf{x})]/kT)$$

where $f(E)$ is the Fermi-Dirac distribution at temperature $T$, $g_{\text{c}}^i$ and $g_{\text{v}}^i$ represent the possible spin and valley degeneracies.

When the simulation is over 1D structure, the wave function can be separated into the plane wave specified with the lattice wave vector ${\mathbf{k}}_{\parallel }$ in the lateral 2D direction and the quantized wave function in the growth direction, which has the ${\mathbf{k}}_{\parallel }$-dependency. Then the charge carrier densitiy is obtained by the following integral over ${\mathbf{k}}_{\parallel }$:

(2.2.9)$$n(x)=\sum _{i\in \text{CB}}{g}_{\text{c}}^i\sum _j\frac{1}{(2\pi {)}^2}{\int }_{{\mathrm{\Omega }}_{BZ}}{d}^2{\mathbf{k}}_{\parallel }{p}_{\text{c},j}^i(x,{\mathbf{k}}_{\parallel })f([E_{\text{c},j}^i({\mathbf{k}}_{\parallel })-E_{\text{F},n}(x)]/kT)$$

(2.2.10)$$p(x)=\sum _{i\in \text{VB}}{g}_{\text{v}}^i\sum _j\frac{1}{(2\pi {)}^2}{\int }_{{\mathrm{\Omega }}_{BZ}}{d}^2{\mathbf{k}}_{\parallel }{p}_{\text{v},j}^i(x,{\mathbf{k}}_{\parallel })f([-E_{\text{v},j}^i({\mathbf{k}}_{\parallel })+E_{\text{F},n}(x)]/kT)$$

Here the integration is over the two-dimensional Brillouin zone ${\mathrm{\Omega }}_{BZ}$.

Similarly, the charge carrier densities for 2D structure is calculated by the integral over the 1-dimensional Brillouin zone as

(2.2.11)$$n(\mathbf{x})=\sum _{i\in \text{CB}}{g}_{\text{c}}^i\sum _j\frac{1}{2\pi }{\int }_{{\mathrm{\Omega }}_{BZ}}d\text{k}{p}_{\text{c},j}^i(\mathbf{x},\text{k})f([E_{\text{c},j}^i(k)-E_{\text{F},n}(\mathbf{x})]/kT)$$

(2.2.12)$$p(\mathbf{x})=\sum _{i\in \text{VB}}{g}_{\text{v}}^i\sum _j\frac{1}{2\pi }{\int }_{{\mathrm{\Omega }}_{BZ}}d\text{k}{p}_{\text{v},j}^i(\mathbf{x},\text{k})f([-E_{\text{v},j}^i(k)+E_{\text{F},p}(\mathbf{x})]/kT)$$

Single-band model

Things are simpler.

When the single-band Schrödinger equation is set to be solved, the envelope function of the $j$-th eigenstate has only one component $F_{l,j}^i(\mathbf{x})$. Also, the k-integration in (2.2.9) to (2.2.12) can be done analytically due to the parabolic dispersion according to the effective mass tensor ${\underset{―}{m}}_e^{\ast i}$ and ${\underset{―}{m}}_h^{\ast i}$.

Thanks to this simpicity the quantum mechanical charge carrier densities for $d$-dimensional simulation can be written up by the following expression:

(2.2.13)$$n(\mathbf{x})=\sum _{i\in \text{CB}}{g}_{\text{c}}^i(\frac{m_{\text{dos,e}}kT}{2\pi {\hslash }^2}{)}^{(3-d)/2}\sum _j{p}_{\text{c},j}^i(\mathbf{x}){\mathcal{F}}_{(1-d)/2}([E_{\text{c},j}^i-E_{\text{F},n}(\mathbf{x})]/kT)$$

(2.2.14)$$p(\mathbf{x})=\sum _{i\in \text{VB}}{g}_{\text{v}}^i(\frac{m_{\text{dos,h}}kT}{2\pi {\hslash }^2}{)}^{(3-d)/2}\sum _j{p}_{\text{v},j}^i(\mathbf{x}){\mathcal{F}}_{(1-d)/2}([-E_{\text{v},j}^i+E_{\text{F},p}(\mathbf{x})]/kT)$$

TODO: The sign in the fermi-dirac integral might be opposite. check the source code.

Here ${\mathcal{F}}_n(E)$ denotes the Fermi-Dirac integral of order $n$ and $m_{\text{dos,}\lambda }^i$ is so-called density-of-states mass defined as

(2.2.15)$$m_{\text{dos,}\lambda }^i=(\text{det}{\overline{m}}_{\lambda }^{\ast i})\lambda =\text{e,h}$$

where ${\overline{m}}_{\lambda }^{\ast i}$ describes the $2\times 2$ or $1\times 1$ submatrix of the effective mass tensor ${\underset{―}{m}}_{\lambda }^{\ast i}$ in the direction of ${\mathbf{k}}_{\parallel }$.

In any cases, the carrier densities are dependent on the electrostatic potential $\varphi (\mathbf{x})$ through the wave function, which is obtained from the $\varphi$-dependent Hamiltonian $H(\varphi )$. Thus we can also write them as $n(\mathbf{x},\varphi )$ and $p(\mathbf{x},\varphi )$, which enters into the non-linear Poisson equation introduced later.

Moreover, when the current equation is included in the calculation scheme, seeing the carrier densities as $n(\mathbf{x},\varphi ,E_{\text{F},n})$ and $p(\mathbf{x},\varphi ,E_{\text{F},p})$ makes it easy to understand what the self-consistent calculation is actually doing.

---

Last update: 04/12/2024