Analytical and numerical solution strategies

1. The analytical route: five steps to a determinant#

For a closed system, such as an electron in a finite well, Lecture 6.1 and earlier lectures of the course followed five steps:

1. write the solution in each region of constant potential as a superposition of two exponentials (plane waves where \(E > U\), growing and decaying exponentials where \(E < U\)): \(2N\) unknown coefficients for \(N\) regions;
2. require the wave function to vanish at \(x \to \pm\infty\), which removes 2 unknowns;
3. require \(\psi\) and \(d\psi/dx\) to be continuous at each of the \(N - 1\) interfaces: \(2N - 2\) homogeneous linear equations for the \(2N - 2\) remaining unknowns;
4. a nonzero solution exists only where the determinant of the coefficient matrix vanishes; find those energies graphically or numerically;
5. normalize \(\int|\psi|^2dx = 1\) to fix the last coefficient.

For a well of width \(W = 6\) nm and depth \(V_0 = 0.4\) eV (mass \(m = 0.067\,m_0\) everywhere), the three regions give four unknowns \((D, A_+, A_-, E')\) in \(\psi = De^{\gamma x}\) (left), \(A_+e^{ikx} + A_-e^{-ikx}\) (well), \(E'e^{-\gamma x}\) (right), with \(k = \sqrt{2mE}/\hbar\) and \(\gamma = \sqrt{2m(V_0 - E)}/\hbar\), and four matching equations at \(x = 0\) and \(x = W\). Step 4 is done here numerically: the determinant of the 4×4 matrix is evaluated on an energy grid and its sign changes are refined with a root finder.

The determinant crosses zero at two energies, 0.0782 and 0.2870 eV, the two bound states of the well. (Its real part vanishes identically, to \(10^{-13}\) of its magnitude, so the imaginary part is plotted.) The same energies follow from the familiar transcendental equations \(\tan(kW/2) = \gamma/k\) and \(-\cot(kW/2) = \gamma/k\) (Section 4 checks them against the discretized Hamiltonian). The method is exact, but each new structure needs its own ansatz and its own matrix, and the size of the matrix grows with the number of regions.

2. The transfer matrix and its limits#

For open systems Lecture 6.1 replaced the determinant by the transfer matrix: one 2×2 matrix per interface, multiplied in order. Klimeck calls it "very appealing", and it works well "for simple structures like these, where you have flat bands, not too many barriers, not too thick of barriers", but "it doesn't work all that well for realistically large devices". Two reasons can be seen directly.

Flat potentials only. The plane-wave ansatz needs a constant potential in each layer. A barrier under an applied bias has a sloped top; the transfer matrix must cut it into \(n\) flat slices. Here a 5 nm, 0.3 eV barrier drops by 0.1 eV across its width (the right contact is at \(-0.1\) eV), and the sliced transfer matrix is compared with open-nemo, which simply puts the sloped potential on its grid.

The error falls as \(1/n^2\) (a factor of four per doubling) until it reaches the discretization error of open-nemo's 0.01 nm grid, about \(2\cdot10^{-5}\), at 64 slices. Slicing works, but every non-flat region of a device (band bending, applied bias, graded alloys) multiplies the number of layers.

Loss of precision. Inside a barrier the transfer matrix carries a growing and a decaying exponential, \(e^{\pm\gamma x}\). Starting from the transmitted wave on the right and multiplying towards the left is stable: the growing exponential is the one that matters. Starting instead from the incident and reflected waves on the left, \([A, B] = [1, r]\), and propagating to the right, the wave function on the right should be the transmitted wave alone, \([t, 0]\). In floating-point arithmetic \(r\) carries a relative error of \(\epsilon \approx 10^{-16}\), and this error feeds the exponential that grows towards the right, by \(e^{2\gamma L}\).

The amplitude of the wave incident from the right, which should be zero, grows exponentially with the barrier width and follows \(\epsilon e^{2\gamma L}\). At a width of 30 nm (\(\gamma L \approx 18\)) the spurious wave reaches 20 % of the transmitted one, and the wave function computed this way is meaningless. Lecture 6.1 showed the same loss in \(\det M\). Stable formulations (scattering matrices, or the Green's functions with open boundary conditions that open-nemo uses) avoid multiplying growing exponentials; Klimeck points to his nanoHUB course for them (Further reading).

3. Discretizing the Schrödinger equation#

The fully numerical route needs no ansatz. Space is replaced by a grid of points \(x_i = x_0 + ia\) with spacing \(a\), and the wave function by its values \(\psi_i = \psi(x_i)\). Adding the Taylor expansions of \(\psi(x_0 + a)\) and \(\psi(x_0 - a)\) cancels the odd terms:

The neglected term is of order \(a^2\). For a plane wave \(e^{ikx}\) the three-point formula gives \(-(2/a^2)(1 - \cos ka)\) instead of \(-k^2\), so the kinetic energy on the grid is

a cosine band. This is the dispersion of a tight-binding chain with one s orbital per site and nearest-neighbour coupling \(-t_0\); Klimeck remarks that the natural grid spacing is therefore the atomic spacing of the crystal, which is what atomistic tight-binding models use. Below, the grid dispersion is compared with the parabola \(\hbar^2k^2/2m\) for three spacings, including one GaAs monolayer (0.2825 nm, half the cubic lattice constant).

The grid underestimates the kinetic energy, by \((ka)^2/12\) to lowest order: at 0.1 eV the error is 0.01 % for \(a = 0.1\) nm, 0.12 % for one monolayer and 0.37 % for 0.5 nm. In this effective-mass model the electron energies of interest correspond to small \(ka\), and the grid is accurate. The grid cannot represent wave numbers beyond \(\pi/a\); its band ends at the zone edge with energy \(4t_0\) (9.1 eV for \(a = 0.5\) nm, the curve that bends over in the left panel).

4. The tridiagonal Hamiltonian#

Inserted into \(-\frac{\hbar^2}{2m}\psi'' + U(x)\psi = E\psi\), the three-point formula gives one equation per grid point,

with \(\psi_0 = \psi_{N+1} = 0\) for a closed system: the terms that would reach outside drop out of the first and last equations. In matrix form \(H\psi = E\psi\), with \(H\) an \(N\times N\) tridiagonal matrix: \(2t_0 + U_i\) on the diagonal and \(-t_0\) on the two neighbouring diagonals.

open-nemo builds the same matrix for a Heterostructure1D. The check below builds \(H\) by hand for the 6 nm, 0.4 eV well inside 20 nm barriers and compares it element by element, then solves the eigenvalue problem (with SciPy's tridiagonal solver, which stores only the two diagonals).

The hand-built matrix and open-nemo's are identical. The eigenvalues converge to the zeros of the determinant of Section 1 and the transcendental equations, and the error falls by about a factor of four each time the spacing is halved: the \(a^2\) error of the three-point formula. Even the coarse 0.4 nm grid is within about 2 meV. The barriers are 20 nm thick so that the hard walls at the ends do not matter: the upper state decays into the barrier with \(\gamma = 0.45\) nm$^{-1}$, and with 10 nm barriers the walls would shift it by about 0.01 meV, more than the grid error at the finest spacing.

5. Closed systems are wrong for open devices#

The eigenvalue problem assumes closed boundaries, \(\psi = 0\) at both ends. For a device from which electrons leak into contacts, this gives "just plain wrong" answers, in Klimeck's words. A double barrier (3 nm barriers of 0.3 eV around a 5 nm well) with GaAs contacts of length \(L_c\) on each side is solved both ways: as a closed box of total length \(2L_c + 11\) nm, and as an open system with nemo.resonances, which attaches semi-infinite contacts through self-energies.

The closed box has a discrete spectrum whose energies depend on the arbitrary length of the contacts: most eigenvalues are standing waves of the contact regions, falling as \(1/L_c^2\) and appearing in nearly degenerate pairs (left and right contact). The state of the well appears as an eigenvalue only where it is not mixed with a contact state, and it is mixed more often the longer the contacts are: with 40 nm contacts there are several eigenvalues near the resonance and none of them is "the" well state. The open system gives one answer, nearly independent of \(L_c\): a resonance at 0.0885 eV with a width of 4.6 to 4.8 meV, a lifetime \(\hbar/\Gamma\) of 0.14 ps, and a continuous spectrum of contact states at every energy. These are the quantities that a transport calculation needs (Lecture 6.3), and they are invisible in the closed box.

6. Exploration: a position-dependent mass#

In a heterostructure the effective mass changes at the interfaces, for example from GaAs (\(0.067\,m_0\)) to AlGaAs (about \(0.092\,m_0\) for 30 % aluminium). Current conservation then requires \(\psi\) and \(\psi'/m\) to be continuous (BenDaniel–Duke), and the kinetic operator becomes \(-\frac{\hbar^2}{2}\frac{d}{dx}\frac{1}{m(x)}\frac{d}{dx}\). Its three-point discretization puts the mass on the bonds between sites, \(t_{i,i+1} = \hbar^2/(2m_{i+1/2}a^2)\), and keeps \(H\) symmetric; this is what open-nemo does. The naive alternative, \(-\frac{\hbar^2}{2m_i}\psi''\) with the mass on the sites, gives a non-symmetric matrix. Both are compared with the exact levels of a 6 nm well of depth 0.3 eV, whose matching condition is \(\tan(kW/2) = m_w\gamma/(m_bk)\) (even states).

The mass-on-bonds discretization converges to the exact levels; at 0.025 nm both states are within 0.01 meV. The naive discretization converges too, but to the wrong answer: 11.8 meV too high for the ground state and 14 meV for the excited state, independent of the grid. It describes a well with \(\psi'\) instead of \(\psi'/m\) continuous, which does not conserve current. Refining the grid does not fix an inconsistent discretization.

Summary#

• Closed systems can be solved analytically: ansatz, boundary conditions, a determinant set to zero, normalization. Every new geometry needs its own matrix.
• The transfer matrix solves open systems with flat layers. Sloped potentials must be sliced (error \(\propto 1/n^2\)), and propagating through thick barriers in the unstable direction amplifies round-off by \(e^{2\gamma L}\).
• On a grid of spacing \(a\) the Schrödinger equation becomes \(H\psi = E\psi\) with a tridiagonal \(H\): \(2t_0 + U_i\) on the diagonal, \(-t_0\) beside it, \(t_0 = \hbar^2/(2ma^2)\). Its error is of order \((ka)^2/12\), and its dispersion \(2t_0(1 - \cos ka)\) is that of a one-orbital tight-binding chain.
• Closed boundaries turn the continuum of the contacts into box states whose energies depend on the contact length. Devices with contacts need open boundary conditions.
• A position-dependent mass must be discretized on the bonds; the naive form converges to wrong energies.

Exercises#

1. Extend coefficient_matrix to a double well (five regions, 8×8 matrix) and find the splitting of the ground state for a 2 nm barrier between two 6 nm wells.
2. Use the stable (backward) propagation of Section 2 to compute the 30 nm wave function, and compare it with open-nemo's spectral_left as in Lecture 6.1.
3. Show that the eigenvalues of the hard-wall box with \(N\) sites are exactly \(E_n = 2t_0(1 - \cos(n\pi/(N+1)))\), and check this against eigh_tridiagonal.
4. Repeat Section 5 with thicker barriers (5 nm). How long must the contacts be before the closed box shows several eigenvalues within one resonance width?

Further reading#

• Lecture 6.1: the transfer matrix method and Lecture 7.1: band structure from a periodic potential.
• Green's functions from scratch and Bound states in open devices: the open boundary conditions that replace the closed box.
• G. Klimeck, "Nanoelectronic Modeling: From Quantum Mechanics and Atoms to Realistic Devices" (2010), nanoHUB resource 8086: the course cited on slide 11.