Overview - Basics 1

Note: This documentation is nearly 10 years old and might not be related any more to the actual version of the nextnano³ source code.
It is still displayed here for historical reasons until a decent manual for nextnano³ will be available.
Interested readers should consult the diploma and PhD theses of S. Hackenbuchner, M. Sabathil, T. Zibold (nextnano++) and T. Andlauer (nextnano++) available from the publications website.

You can dowload this page as a pdf file (basics.pdf) but most of it is written in German and this file won't be updated any more. Be sure to check the online documentation for the latest version.

Contents

1 Basics 1 -Overview
1.1 Main equations
1.1.1 Boundary conditions
          Boundary conditions for the Schrödinger equations
          Boundary conditions for the Poisson and current continuity equation
1.1.2 Calculation of the potential and the densities at fixed quasi-Fermi level
        1. Purely classical calculation of the density
        2. Purely quantum mechanical calculation of the density
        3. Mixed classical/quantum mechanical calculation of the density whereby
            only bound states contribute to the quantum mechanical density
        4. Mixed classical/quantum mechanical calculation of the density whereby
            in quantum regions the density is calculated only quantum mechanically
1.1.3 Calculation of the quasi-Fermi level at fixed potential
1.2 General structure
1.3 Generation of the grid
1.4 Control of the program flow via the input file
1.5 Technical details for the states in quantum regions
1.5.1 Separation of classical and quantum mechanical states
1.6 Models for density calculation in 1D, 2D and 3D
1.6.1 Classical densitites
1.6.2 Density for 1-band quantum states
1.6.2.1 3D states
1.6.2.2 2D states (homogeneous in z direction)
1.6.2.3 1D states (homogeneous in (x, y) direction)
1.6.3 Derivatives of Densities
1.6.4 k.p densities
1.6.4.1 Evaluation of the density integration
1.6.4.2 1D integration
1.6.4.3 2D integration, discretization of k space
1.6.4.4 Density of states integration
1.6.4.5 k.p dispersion without density integration
1.7 Details for the interior program flow
1.7.1 Rough structure
1.7.2 Schrödinger-Poisson equation
1.7.3 Calculation of multi-band k.p states
1.7.4 Limitations
2 Basics 2 - Definition of k.p Hamiltonian
2.1 Wurtzite
2.1.1 Spin-orbit interaction
2.1.2 Crystal field splitting
2.1.3 Strain effects
2.1.4 Transformation to rotated coordinate system
2.1.5 Special wurtzite parameters
2.2 Zincblende
2.2.1 Spin-orbit interaction
2.2.2 Strain effects
2.2.3 Transformation to rotated coordinate system
2.3 The rotation problem
2.3.1 1D simulation
2.4 Discretization

1 Basics 1 - Overview

1.1 Main equations

First, the main equations are summarized. To keep things simple, we restrict ourselves to electrons and donors and do not treat recombination terms. Holes and acceptors will be treated analogously. Furthermore the following equations refer to 1D device geometries. For the solution of the system we take a quasi-Fermi level E_Fn(z) which varies over space.

For the density of the ionized donors and for the current hereby holds:

N_D(z) is the donor density, E_D is the inonization energy, g_D is the degeneracy factor of the energy levels and µ_n(z) is the electron mobility.

Equation (1.1) is the Poisson equation. Hereby n stands for the electron density and N_D⁺ is the density of the ionized donors. The 1-band Schrödinger equation (1.2) gives the density of states for the quantized states (that have the wave functions Psi_n(z) and the corresponding eigenvalues E_n) whereby V(z) = E_c0(z) - e Phi(z) is the band edge function. The electron density n(z) which goes into the Poisson equation is obtained
- by occupying the quantized states Psi_n(z) using the quasi-Fermi level E_Fn(z)
- or by the Thomas-Fermi approximation, respectively (in case of a classical calculation).
More details follow below...

The current continuity equation (1.3) provides the connection between the electron density and the quasi-Fermi level. Here again, more details follow below...

One realizes that the self-consistent solution of this system of equations is uniquely determined by the potential Phi(z) and the quasi-Fermi level E_Fn(z).

Methodically, therefore, our approach is divided into two blocks.

• Block 1: Poisson block
  The Schrödinger and Poisson equations are solved self-consistently at fixed, space dependent quasi-Fermi level resulting in the new potential and the new quantized states.
• Block 2: Current block
  The new quasi-Fermi level is obtained by solving the current continuity equation at fixed potential and fixed states.

This cycle will be repeated until self-convergence is reached.

Both parts will be discussed below in more detail.

1.1.1 Boundary conditions

Boundary conditions for the Schrödinger equation

For the Schrödinger equation, Dirichlet boundary conditions are used, i.e.

and Neumann boundary conditions, i.e.

respectively. At first, the effect of these boundary conditions is studied at a shallow, one-dimensional potential function on a region of 10 nm length (see Fig. 1).

With Dirichlet boundary conditions, the eigenfunctions have cosine shape as it is the case at a potential well with infinitely high walls. If one takes Neumann boundary conditions, then one obtains the same eigenvalues but the corresponding eigenfunctions are shifted by a phase of pi/2. Additionally, the eigenvalue zero occurs with a constant eigenfunction.

Fig.1: Eigenvalues of a shallow potential of 0 eV with Dirichlet and Neumann boundary conditions. The lowest Neumann eigenvalue has the value zero and thus lies exactly onto the bottom of the potential.

If the potential isn't shallow but sloped in a way, then we have the situation as in Fig. 2.

Fig. 2: Eigenvalues at a sloped potential with Dirichlet and Neumann boundary conditions

If one solves the Schrödinger equation only in an inner region of the device, then one can either take Dirichlet or Neumann boundary conditions. Eventually, the states are normalized to one, i.e.

This means that a fully occupied state is equivalent to exactly one electron. At contacts we are making the idealistic assumption that there exists charge neutrality and that the density of states is the same as in bulk material. With pure Dirichlet boundary conditions, however, the density would decrease to zero at contacts whereas with Neumann boundary conditions the density would increase. Thus, it is reasonable to assume something like a mixed state. Precisely, this means that we solve the Schrödinger equation once with Dirichlet and once with Neumann boundary conditions and to normalize the states to 1/2:

This can be made plausible by looking at Fig. 2. For the same eigenvalue, one obtains there once cosine and once a sine function. The sum over both is constant. As both these eigenfunctions are normalized to 1/2, the fully occupation of these mixed states again corresponds exactly to one electron.

Boundary conditions for the Poisson and current continuity equation

The solution of the Schrödinger equation (and of the Thomas-Fermi approximation at classical calculations respectively) gives the density of states for the device. This density of states will be occupied using the quasi-Fermi level. To guarantee charge neutrality, we are first solving the Poisson and Schrödinger equation self-consistently for the case of thermodynamic equilibrium with constant quasi-Fermi level E_F0 whereby we require as boundary conditions at contacts:

     (Neumann boundary condition)

This corresponds to the assumption of charge neutrality at Ohmic contacts. Then, according to Gauss' theorem the total charge is zero, i.e. the charge of the free charge carriers is equal to the charge of the ionized donors and acceptors. From this, one obtains the built-in potential. The values of this potential at the contacts Phi_l (left contact) and Phi_r (right contact) are then taken as Dirichlet boundary conditions for the calculations in nonequilibrium.

Precisely, this means by applying a voltage U:

This leads to the result that the charge carrier density at an Ohmic contact at applied voltage should be equal to the charge carrier density at equilibrium.

For the calculation of the built-in potential for a Schottky contact, the band edge will be pinned above the Fermi level by a Schottky barrier. This corresponds to a Dirichlet boundary condition. The calculation at applied voltage then is analogous to the one for an Ohmic contact.

1.1.2 Calculation of the potential and the densities at fixed quasi-Fermi level

We assume that there exists a quasi-Fermi level E_Fn(z) which varies over space and is valid for all states. First, we are calculating the densities where the potential Phi(z) is fixed and thus the function of the band edge E_c0(z) - e Phi(z) is fixed as well. Hereby four cases have to be considered:

1. Purely classical calculation of the density

In this case, the charge carrier density is calculated via the Thomas-Fermi approximation which deals with the fully three-dimensional density which is, however, homogeneous in the x and y direction:

(1.10)

Hereby, g_valley stands for the degeneracy factor of the valley and m^* for the effective mass of the electron. The Fermi function F_1/2 is defined as:

(1.11)

 

2. Purely quantum mechanical calculation of the density

By solving the Schrödinger equation over the whole device one gets the quantized states where E_n is the energy and Psi_n(z) the wave function of the n^th eigenstate. Eventually these states will be occupied with the local quasi-Fermi levels. The three-dimensional charge carrier density n(z) which is homogeneous in x and y direction is obtained from:

Hereby, g_valley stands for the degeneracy factor of the valley and m^* for the effective mass of the electron. Here, the quantum region 'touches' both contacts. Thus, we are using for the solution of the Schrödinger equation mixed Dirichlet and Neumann boundary conditions as was described above.
Note that this formula leads to discontinuities in the quantum mechanical charge density if m*(z) has discontinuities (e.g. at material interfaces). A solution is to average m*(z) over all parallel masses in the quantum cluster weighted by the probability of the wave functions.

 

3. Mixed classical/quantum mechanical calculation of the density whereby only bound states contribute to the quantum mechanical density

Often one wants to calculate quantum mechanically only in certain regions of the device. For this, we are dividing the device into classical regions (CR) and quantum regions (QR).

In the classical regions, the density is calculated via the Thomas-Fermi approximation.

In quantum regions, first, the Schrödinger equation is solved by using Dirichlet boundary conditions, then the separation energy (E_separation) up to which the quantum mechanical states are considered, is determined. The density of the bound (localized) states is called n_bound. Then it holds:

The density for energies higher than E_separation (continuum density n_cont) is calculated via Thomas-Fermi approximation:

 

Hereby, F_1/2 stands for the incomplete Fermi integral:

	(1.12)
The total density in the quantum region is then given by:

 

4. Mixed classical/quantum mechanical calculation of the density whereby in quantum regions the density is calculated only quantum mechanically

As in the previous paragraph, we are dividing the device into classical (CR) and quantum (QR) regions.

In the classical regions, the density is calculated via the Thomas-Fermi approximation.

In quantum regions, first, the Schrödinger equation is solved by using mixed Dirichlet and Neumann boundary conditions. The density in the quantum regions is calculated by occupying these states over the local quasi-Fermi level. If one uses only Neumann boundary conditions, one obtains a density of states which increases at the boundary of the quantum region resulting in substantial discontinuities in the density function at the transition from classical to quantum regions. It holds for the quantum mechanical density:

The densities are then calculated via a function of the potential (and thus of the quantum mechanical states). With these terms for the densities, the Poisson and Schrödinger equation are then solved self-consistently, e.g. by a predictor-corrector approach. Thus for the fixed quasi-Fermi level a new potential and new quantum states are obtained.

1.1.3 Calculation of the quasi-Fermi level at fixed potential

Starting from the Boltzmann equation we get for the current density:

Hereby µ_n(z) stands for the mobility of the electrons. From the solution of the current continuity equation

div j = 0

one eventually gets the quasi-Fermi level E_Fn(z)

Hereby, the charge carrier density goes into the current continuity equation. Here again, as in the previous chapters, we have to consider different cases.

- For the case of purely classically or purely quantum mechanically calculated density the whole density goes into the current continuity equation.

- If the density is calculated mixed classically/quantum mechanically, there are two possibilities.

- If one assumes that the bound states do not contribute to the current flow, in the current continuity equation only the continuum density will be considered:
n(z) = n_cont(z)

- Often, it is not possible to distinguish between bound (localized) and delocalized states, especially in higher dimensions. Then the whole density will be used in the current continuity equation implying that the current will change substantially.

1.2 General structure

The general structure is described here.

1.3 Generation of the grid

The generation of the grid is described here.

1.4 Control of the program flow via the input file

The control of the program flow (flow-scheme) via the input file is described here.

1.5 Technical details for the states in quantum regions

1.5.1 Separation of classical and quantum mechanical states

The separation model (quantum and continuum states) is described here.

1.6 Models for density calculation in 1D, 2D and 3D

1.6.1Classical densitites

One can derive the Thomas-Fermi density as follows:
The density of states can be written as:

(1.14)

The wavefunctions for propagating waves are:

With

(eq. 1.14)

one obtains for the electron density:

We assume a parabolic dispersion:

Summation can be transformed to integration by:

Then the electron density reads:

(1.15)

With

one gets:

(1.16)

We further allow for nonparabolicities by introducing a nonparabolicity factor. The nonparabolicity factor in the database has following definition:

It has by definition the units [alpha_eV]=1/eV. If we carry out the following transformation:

alpha is dimensionless and the density now reads:

(1.17)

In case we have a combination of quantum and classical states one only has to calculate the classical density above an energy level (e.g. the highest bound state), which we call E_edge. As a consequence, we only have to integrate from E_edge (instead of E_c) to infinity in (eq. 1.14). Thus we get incomplete Fermi integrals:

 

		(1.18)

For holes one only has to replace   (E_F - E_c) / kT   by   -(E_F - E_v) / kT  .

1.6.2 Density for 1-band quantum states

1.6.2.1 3D states

With the three dimensional eigenstates from the Schrödinger equation the electron density reads:

(1.19)

1.6.2.2 2D states (homogeneous in z direction)

The quantum states are the eigenstates of the one-band Schrödinger equation. We will derive the according density for a quantum region which is homogeneous in the z direction:

Due to homogeneity in z direction the wavefunctions can be separated:

Then for the energies the relation holds:

In general it is

(eq. 1.14),

so the electron density reads:

where f(E_n) denotes the Fermi function. Due to the separation the summation can be written as:

Thus one obtains:

Finally the density is (provided one integrates to infinity)

(1.20)

In general we have an energy edge E_edge which separates quantum mechanical and classical states. Here we have to integrate for the bound states from the eigenenergy of the considered state to E_edge, so:

This can be written as:

1.6.2.3 1D states (homogeneous in (x, y) direction)

With one dimensional wavefunctions the density is

(1.22)

1.6.3 Derivatives of Densities

In order to calculate the Newton correction for the nonlinear Poisson equation, derivatives of the densities with respect to the potential have to be calculated. Here, one makes use of the relation

Note that for incomplete Fermi integrals the boundaries may depend on the potential, so one has to make use of the formula:

(1.23)

It is

for

Thus after partial integration one gets:

		(1.24)

The Gamma function has the values:

	(1.25)
	(1.26)
	(1.27)
	(1.28)

More information on this topic can be found here.

1.6.4 k.p densities

(More details about the kp densities can be found here as well.)

The wave functions are given by:

The density of states is

It holds

One can get rid of the functions u_iu_j* in the following way. When calculating the density at position (x, y, z), we mean the averaged density in the box at position (x, y, z). Furthermore, we assume that the functions u_iu_j* vary very rapidly, so we can separate the integral:

As boundary conditions for the k.p states, Dirichlet and periodic boundary conditions are implemented.

1.6.4.1 Evaluation of the density integration

The evaluation of the density integration is describe here (method of Brillouin zone integration).

1.6.4.2 1D integration

If the energy is isotropic inside the 2D Brillouin zone, i.e. E(k_x,k_y)=E(k_||) with k_||²=k_x²+k_y², then the integration can be reduced to a one-dimensional integration:

1.6.4.3 2D integration, discretization of k space

The integration is written as follows:

1.6.4.4 Density of states integration

Another possibility is to integrate over the density of states:

1.6.4.5 k.p dispersion without density integration

The program chooses the maximum k vector up to which the states are calculated automatically (see above). However, if one wants to know the dispersion up to a given k vector, one has to specify the value in the input file under the keyword $output-kp-data and its specifier 'k-par-disp-range' [1/Angstrom]. Then the k.p states will be calculated again in the end and put out.

More details about the kp densities can be found here.

1.7 Details for the interior program flow

1.7.1 Rough structure

The rough program flow is described here and more detailed here.

1.7.2 Schrödinger-Poisson equation

During this part of the program the quasi-Fermi levels that are given in the module 'fermi_level' don't change. Thus, the densities are obtained by occupying the states over these quasi-Fermi levels. For the self-consistent solution of the Schrödinger and Poisson equation we are using a predictor-corrector method. The program flow is the following:

First, for the current valid potential the quantum mechanical states are determined. From these a predictor density can be deduced. In essence, one uses a perturbation ansatz in first order, i.e. one assumes that the wave functions remain unchanged at variations of the potential and that the eigenvalues are shifted according to the potantial changes. With this predictor density, the nonlinear Poisson equation will be solved and thereby a new (corrected) potential is obtained. From this, one can immediately determine the new eigenstates aso. This procedure will be repeated until convergence.

The outer cycle is situated in the subroutine 'non_linear_poisson'. From there, the subroutines
- 'calculate_eigenstates' for the calculation of the current eigenstates and
- 'newton' for the calculation of the nonlinear Poisson equation (including the predictor density)
are called alternately

More information is available for the nonlinear Poisson equation (section inside the Poisson block) and for the calculation of the eigenstates.

1.7.3 Calculation of multi-band k.p states

The calculation of the eigenstates is controlled centrally by the sub program 'calculate_eigenstates'. This in turn calls for kp regions the routine 'update_kp'. 'update_kp' is the superior program for the k.p problem.
First the eigenstates for k_||=0 are determined. This occurs by calling 'setup_kp_mat' (setup_kp_mat1D) (setup of the matrix) and then calling 'calculate_kp_eigenvalues' (calculation of eigenstates).
In 'get_range_k' is is decided up to which maximum value of k_|| the Brillouin zone has to be discretised. Then the values for k_|| are being determined ('get_k_par2D' (1D)) and finally the eigenstates for these are determined again.

More details can be found here.

1.7.4 Limitations

In this version, nextnano³ is treating all material types in the same way. This means specifically that the Input parser requires the full parameters even for metals and air and even if they don't contribute to the calculation (e.g.  the metal could be defined additionally as a Poisson region which represents an Ohmic contact and thus appearing in the Poisson and Schrödinger equation only as a contact).

For the k.p Hamiltonian up to now only Dirichlet and periodic boundary conditions are implemented. If one wants to calculate the current with full quantum mechanical density, one should use the 1-band Schrödinger equation instead. Mobility models are implemented but no recombination terms so far. Particularly, for a nonequilibrium calculation, contacts have to be specified on both ends in order to allow current flow. (If not, no current can flow because of current conservation.)

For the solution of eigenvalue problems, LAPACK library routines are used which are specified in the MODULE numerical_control1D, as well as iterative methods, like e.g. ARPACK.

Please go to Basics 2 to continue reading!