Diskretisierung

Bei der numerischen Lösung von Differentialgleichungen wird der Träger der gesuchten Lösung zerlegt in finite Differenzen [Sel84]. Unser Modell ist eindimensional, daß heißt, das Bauteil hat nur zwei Oberflächen: den linken Rand $L$ und den rechten Rand $R$. Die Diskretisierung des kompakten Trägers ist in Abbildung 1.4 dargestellt. Der Abstand zweier benachbarter Stützstellen ist für das gesamte Bauteil nicht äqui-distant (inhomogenes Gitter).

Abbildung 1.4: Diskretisierung des kompakten Trägers der Poisson-Gleichung.

Da im Gegensatz zu einer analytischen Lösung, die auf einem kontinuierlichen Träger lebt, die numerische Lösung auf einer endlichen Anzahl von Stützstellen gegeben ist, müssen die Werkzeuge der Infinitesimalrechnung entsprechend modifiziert werden. Bei der Ableitung einer Funktion $f$ wird aus dem Differentialquotient der Differenzenquotient

$$\left. \partial_z f (z) \right\vert _{z_m} = \left. \frac{df}{dz}... ... z} \right\vert _{z_m} = \frac{ f(z_{m+1}) - f(z_{m-1})}{h_m + h_{m+1}} \qquad,$$ (1.59)

mit dem Index der Stützstellen m und dem Abstand zweier Stützstellen $h_m := z_{m+1} - z_m$. Bei der Integration wird anstelle der Grenzwertbetrachtung einer Riemann-Zwischensumme die endliche Summation über Treppenfunktionen durchgeführt

$$\int_0^z f ( z ) \ dz = \sum_m f(z_m) \left( \frac{h_{m-1} + h_m}{2} \right) \qquad.$$ (1.60)

Zum Lösen der Poisson-Gleichung (1.10) wird diese über jeweils eine Stützstelle $z_m$ stückweise integriert [Sel84], der Übersichtlichkeit halber gelte die Konvention $z_{m + 1/2} = z_m + (h_m / 2)$,

$$\epsilon_0 \partial_z \lbrack \epsilon( z ) \partial_z \phi ( z ) \rbrack = - \rho ( \phi ( z ) )$$ (1.61)

$$\Rightarrow \int_{z_{m-1/2}}^{z_{m+1/2}} \epsilon_0 \partial_z \lbrack \epsilon( z ) \partial_z \phi ( z ) \rbrack \ dz = - \int_{z_{m-1/2}}^{z_{m+1/2}} \rho ( \phi ( z ) ) \ dz$$ (1.62)

$$\overset{(\ref{DefNumInt})}{\Rightarr... ...z_{m-1/2}} = - \rho ( \phi (z_m) ) \left( \frac{{h_{m-1}} + {h_{m}}}{2} \right)$$

$$\overset{(\ref{DefNumAbl})}{\Rightarrow} \epsilon_0 \epsilon (z_{m+1/2}) \frac{\phi(z_{m+1})-\phi(z_m)}{h_m} - \epsilon_0 \epsilon (z_{m-1/2}) \frac{\phi(z_{m})-\phi(z_{m-1})}{h_{m-1}} =$$

$$- \rho ( \phi (z_m) ) \left( \frac{{h_{m-1}} + {h_{m}}}{2} \right)$$ (1.63)

$$\Leftrightarrow \epsilon_0 \epsilon (z_{m-1/2}) \frac{\phi(z_{m-1... ...ac{ \epsilon (z_{m+1/2})}{h_{m}} + \frac{\epsilon (z_{m-1/2})}{h_{m-1}} \right]$$

$$+ \epsilon_0 \epsilon (z_{m+1/2}) \frac{\phi(z_{m+1})}{h_{m}} = - \rho ( \phi (z_m) ) \left( \frac{{h_{m-1}} + {h_{m}}}{2} \right) \qquad.$$ (1.64)

Dadurch läßt sich der diskretisierte, eindimensionale Laplace-Operator in Matrixform darstellen

$$\epsilon_0 \partial_z \lbrack \epsilon( z ) \partial_z \dots \rbr... ...uad \qquad\qquad \qquad \qquad \qquad \qquad \qquad \qquad \qquad \qquad \qquad$$

$$\epsilon_0 \underbrace{\left( \begin{array}{ccccccc} - \frac{... ... \ddots & \ddots \end{array}\right)}_{=:\underline{\mathbf{A}}}$$

Zusammen mit dem Vektoren $\mathbf{\Phi}_{i+1}$ für das Potential $\phi$ und $\mathbf{b}_{i+1}$ für die Inhomogenität

$$\mathbf{\Phi}_{i+1} = \left( \begin{array}{c} \phi_{i+1} (z_1) \\... ...z_n) ) \left( \frac{{h_{n-1}} + {h_{n}}}{2} \right) \end{array} \right) \qquad,$$

wobei n die Zahl der Stützstellen ist, läßt sich die eindimensionale Poisson-Gleichung (1.10) für den Iterationsschritt $i+1$ als lineares Gleichungssystem in Matrixform schreiben

$$\epsilon_0 \underline{\mathbf{A}} \mathbf{\Phi}_{i+1} = -\mathbf{b}_{i+1} \qquad.$$ (1.65)

Die Matrix $\underline{\mathbf{A}}$ hat eine tridiagonale Bandform und wird durch Anwenden des Gauß-Algorithmus [Pre92] invertiert. Dadurch kann zu jeder gegebenen Ladungsträgerverteilung $\mathbf{\rho}$ ein Potential $\mathbf{\phi}$ bestimmt werden.

Um schnell, also mit minimaler Anzahl an Iterationsschritten, eine konvergierte Lösung zu erhalten, wird das Newton-Verfahren eingesetzt (Abschnitt 1.5.1). Für den neuen Differential-Operator $\tilde \Delta$ in (1.45) wird die Ableitung der Ladungsträgerdichten $\rho$ nach dem Potential $\phi$ benötigt. Aus (1.29), (1.30), (1.34) und (1.35) ergeben sich

$$\left. \frac{\partial n^{3d}}{\partial \phi} \right\vert _{\phi ( z ) } = e \beta N^{3d}_c (z) F_{-1/2} ( \beta [E_{Fn}(z) - E_c(z)] )$$ (1.66)

$$\left. \frac{\partial p^{3d}}{\partial \phi} \right\vert _{\phi ( z ) } = e \beta N^{3d}_v (z) F_{-1/2} ( \beta [E_v(z) - E_{Fp}(z)] )$$ (1.67)

$$\left. \frac{\partial N^+_D}{\partial \phi} \right\vert _{\phi ( z ) } = - e \beta \frac{N_D (z) 2 exp(\beta[E_{Fn}(z) - E_D(z)])}{(1 + 2 exp(\beta[E_{Fn}(z) - E_D(z)]))^2}$$ (1.68)

$$\left. \frac{\partial N^-_A}{\partial \phi} \right\vert _{\phi ( z ) } = e \beta \frac{N_A (z) 4 exp(\beta[E_A(z) - E_{Fp}(z)])}{(1 + 4 exp(\beta[E_A(z) - E_{Fp}(z)]))^2} \qquad,$$ (1.69)

und daraus

$$\left. \frac{\partial \rho}{\partial \phi} \right\vert _{\phi ( z... ...) } + \left. \frac{\partial N_D^+}{\partial \phi}\right\vert _{\phi ( z ) } ] ,$$

wobei die Rechnung zur Ableitung des Fermi-Dirac-Integrals in Anhang C zu finden ist.

Unter Verwendung der Jakobi-Matrix $\underline{\mathbf{J}}$ [Wüs95]

$$\underline{\mathbf{J}} = \left( \begin{array}{ccccc} \left. \frac... ...al \rho_{i+1}}{\partial \phi} \right\vert _{\phi_i (z_n ) } \end{array} \right)$$ (1.70)

wird die diskretisierte Poisson-Gleichung (1.66) für das Newton-Verfahren modifiziert

$$[\epsilon_0 \underline{\mathbf{A}} + \underline{\mathbf{J}}] \mat... ...}_{i+1} = - \mathbf{b}_{i+1} + \underline{\mathbf{J}} \mathbf{\Phi}_{i} \qquad.$$ (1.71)

Es müssen nun noch die Randbedingungen verwendet werden. Aus der Dirichlet-Randbedingung für die Poisson-Gleichung (1.8) folgt

$$\phi(z_1) = \phi_L$$ (1.72)

wobei $\phi_L$ der feste Wert des Potentials am linken Rand ist. Damit und aus (1.65) für $\phi(z_m) = \phi(z_2)$ folgt

$$\epsilon_0 \epsilon (z_{\frac{3}{2}}) \frac{\phi_{L}}{h_{1}} - \e... ...{3})}{h_{2}} = \rho ( \phi (z_2) ) \left( \frac{{h_{1}} + {h_{2}}}{2} \right) .$$

Die ersten beiden Zeilen des Gleichungssystems (1.66) schreiben sich für die Dirichlet-Randbedingung damit um zu

$$\epsilon_0 \left( \begin{array}{ccccc} 1 & 0 & 0 & 0 & \dots\... ..._{i+1}(z_1) \\ \phi_{i+1}(z_2) \\ \vdots \end{array}\right) =$$

$$- \left( \begin{array}{c} \phi_L \\ \rho_{i+1} ( \phi_i (z_2... ...rac{3}{2}}) \phi_L}{h_{1}}\\ \vdots \end{array}\right) \qquad.$$ (1.73)

Analog folgt für die letzten beiden Zeilen

$$\epsilon_0 \left( \begin{array}{ccccc} \vdots & \vdots & \vdo... ...\ \phi_{i+1}(z_{n-1}) \\ \phi_{i+1}(z_n) \end{array}\right) =$$

$$- \left( \begin{array}{c} \vdots \\ \rho_{i+1} ( \phi_i (z_{... ..._{n+\frac{1}{2}}) \phi_{R}}{h_{n}}\\ \phi_R \end{array}\right)$$ (1.74)

mit $\phi_R$ dem festen Wert des Potentials am rechten Rand.

Abbildung 1.5: Angenommener Leitungsbandkantenverlauf erzeugt durch eine Oberflächenladung $\sigma$ am linken Rand $L$ bei $z_1$.

Bei der Neumann-Randbedingung (1.9) ist die Ableitung des Potentials, also das Feld E, am Rand durch die Oberflächenladung $\sigma$ gegeben. Da man dort nur eine Stützstelle zur Konstruktion der Ableitung hat, betrachtet man zuerst den linksseitigen und rechtsseitigen Differenzenquotienten

$$\phi'_{links} = \frac{\phi(z_1) - \phi(z_0)}{h_0} = E_{au} \qquad \land \qquad \phi'_{rechts} = \frac{\phi(z_2) - \phi(z_1)}{h_1} = E_{in}\qquad.$$ (1.75)

Unter der Annahme, daß am Ort L der Potentialverlauf spiegelsymmetrisch ist (siehe Abbildung 1.5, vergleiche gegebenenfalls mit Abbildung 3.3), folgt für das innere Feld $E_{in}$ und das äußere Feld $E_{aus}$

$$E = E_{in} = - E_{aus}$$ (1.76)

und damit

$$\frac{\phi(z_1) - \phi(z_0)}{h_0} \overset{!}{=} - \frac{\phi(z_2) - \phi(z_1)}{h_1} = \frac{\phi(z_1) - \phi(z_2)}{h_1} \qquad.$$ (1.77)

Durch stückweise Integration der Poisson-Gleichung (1.10) über den linken Rand findet man

$$\int_{z_{1-1/2}}^{z_{1+1/2}} \epsilon_0 \partial_z \lbrack \epsilon( z ) \partial_z \phi ( z ) \rbrack \ dz = - \int_{z_{1-1/2}}^{z_{1+1/2}} \sigma ( \phi ( z ) ) \delta (z - z_1) \ dz$$ (1.78)

$$\Rightarrow \epsilon_0 \epsilon(z_1) \left[ \partial_z \left. \phi \right\vert _{z_{3/2}} - \partial_z \left. \phi \right\vert _{z_{1/2}} \right] = - \sigma ( \phi ( z_1 ) )$$ (1.79)

$$\Rightarrow \epsilon_0 \epsilon(z_1) \left[ E_{aus} - E_{in} \right] = - \sigma ( \phi ( z_1 ) )$$ (1.80)

$$\overset{(\ref{FeldNeumannRand})}{\Leftrightarrow} 2 \epsilon_0 \epsilon(z_1) E = - \sigma ( \phi ( z_1 ) ) \qquad.$$ (1.81)

Das und verwenden von (1.80) ergibt

$$\phi(z_1) \frac{2 \epsilon_0 \epsilon (z_{1}) }{h_1} - \phi(z_2) \frac{2 \epsilon_0 \epsilon (z_{1}) }{h_1} = - \sigma ( \phi ( z_1 ) ) \qquad.$$ (1.82)

Somit schreiben sich die ersten beiden Zeilen des Gleichungssystems (1.66) für die Neuman-Randbedingung um zu

$$\epsilon_0 \left( \begin{array}{ccccc} \frac{2 \epsilon_0 \ep... ..._{i+1}(z_1) \\ \phi_{i+1}(z_2) \\ \vdots \end{array}\right) =$$

$$- \left( \begin{array}{c} \sigma_{i+1} ( \phi_{i} ( z_1 ) ) \... ...}) \phi_{i+1}(z_1)}{h_{1}}\\ \vdots \end{array}\right) \qquad.$$ (1.83)

Ein Verfahren zur Lösung der Strom-Gleichungen ist in Abschnitt 1.5.2 ausführlich besprochen worden. Zusammen mit diesem Algorithmus und der Definition des Integrals auf diskreten Stützstellen (1.60) können damit die Quasi-Fermi-Niveaus berechnet werden.