Electrostatics#

Solving the Poisson equation self-consistently with the open-boundary Schrödinger equation is a key part of quatrex's quantum transport formalism. By inducing a non-equilibrium charge density in the device, the electrostatic potential will change, which in turn changes the charge density. This feedback loop is solved self-consistently until convergence.

Excess Charge Approach#

We employ an excess charge approach in the electrostatic model implemented in quatrex. The excess charge density \(\rho(\mathbf{r})=n(\mathbf{r}) - p(\mathbf{r})\) is computed from the occupied/depleted electronic states around a charge neutrality level \(E_{CNL}(\mathbf{r})\). For NEGF we can write

For the wave function transport formalism we get a very similar formula, where instead of the lesser/greater Green's functions we have the contact density of states multiplied by the respective Fermi-Dirac occupancies.

In either case, the electrostatic potential \(\phi(\mathbf{r})\) is then obtained via the Poisson equation

where \(\varepsilon(\mathbf{r})\) is the dielectric permittivity. We iterate between computing transport and electrostatics until reaching convergence.

Non-linear Poisson Equation#

Since the excess charge density is a non-linear functional of the electrostatic potential, \(\rho[\phi(\mathbf{r})]\), the Poisson problem we solve has a nonlinear (or more precisely a semilinear) form

To tackle this nonlinearity, we employ a Root-Finding approach (i.e. Newton-Raphson, or also predictor-corrector scheme in this context). We look for the roots of the functional

via the recursion relation

At every step of this scheme we need to evaluate \(\rho[\phi_n(\mathbf{r})]\) and \(\frac{\delta \rho[\phi_{n}(\mathbf{r})]}{\delta \phi_{n}(\mathbf{r})}\). Since evaluating this accurately using the coupled transport formalism would be very costly, we opt to model the charge density's response to potential variations by making a few approximations.

Density Response Models#

By approximating the electronic structure around our charge neutrality level by single parabolic bands (one each for electrons and for holes) we can write

Here \(N_{ND}(\mathbf{r})\) is the effective density of states of an \(N\)-dimensional system, \(\eta_{n}(\mathbf{r})\) is a reduced electrochemical potential describing the distance between the band onset and the CNL, and \(\mathcal{F}_{k}(\eta)\) is the complete Fermi-Dirac integral of order \(k = N / 2 - 1\):

See also the very useful "Notes on Fermi-Dirac Integrals"¹.

Evaluating the charge density and its derivative with respect to the potential in this picture boils down to evaluating Fermi-Dirac integrals:

where we can use the fact that the derivative of a Fermi-Dirac integral of order \(k\) is another Fermi-Dirac integral of order \(k-1\).

We also need to determine a CNL that is consistent with this model. Since the since the CNL directly depends on the potential, we can do so by inverting the relationship between density and potential:

While in two-dimensional systems (\(N=2\)) the Fermi-Dirac integral, its derivative, and its inverse have an analytic form, they have to be evaluated numerically for density models with other dimensionalities.

Contact Chemical Potentials#

In quantum transport simulations, we exert control over the system through its contacts. In self-consistent Schrödinger-Poisson (SCSP) runs, determining appropriate and physically consistent contact chemical potentials is therefore very important.

In quatrex, the most appropriate way of configuring contacts for SCSP runs, is to set a contact's bias voltage \(V_b\) and a mid-gap energy (this could, for instance, be the contacts DFT Fermi level) and then have quatrex compute the contact's Fermi level \(E_F\) and in turn its chemical potential \(\mu = E_F - V_b\) from this.

Connecting Real-Space and Localized Orbital Basis Sets#

In the transport part of an SCSP run, we operate on quantities expressed in a basis of localized orbitals. The electrostatics, on the other hand are treated in (discretized) real-space. Taking the Green's function as an example, the transformation from real-space to a localized orbital basis \(\{\psi_i\}\) is

and the corresponding inverse projection is

where \(\hat{\psi}_{m} = \sum_k \left(\mathbf{S}^{-1}\right)_{mk} \psi_{k}\) are the dual basis states for the general case of an non-orthonormal basis.

Mulliken Charge Projection#

We employ a Mulliken charge analysis to construct the real-space charge density. In the case of non-orthonormal basis sets this can be found by matrix multiplication of the density matrix \(\boldsymbol{\rho}\) with the overlap matrix \(\mathbf{S}\).

Finite Element Discretization#

In quatrex we discretize the Laplacian using linear (first-order) tetrahedral finite elements as implemented in scikit-fem. By default natural boundary conditions are used. Parameters concerning regions in the real-space simulation domain and the meshing process (using gmsh) are set in the [device.geometry] section of the configuration.

The structure.xyz input file containing the atom/orbital coordinates that make up the device, serves as a reference point for all defined geometry regions. This file also informs the periodicity in transverse directions through the pbc="..." entry on the second line in the extended .xyz file format. Setting pbc="F F F", the real-space mesh will not enforce any periodicity, while pbc="T F F" will lead to periodic boundary conditions being enforced along \(x\)-direction.

Each 3D region can define properties such as relative permittivity, and doping concentrations, while 2D regions in the simulation domain act as gates/ground planes with work functions and set voltages. Orbital sites are embedded into the mesh to accomodate the Mulliken charge projection.

Dirichlet Boundary Conditions#

We model electrostatic control from gates in the simulation by imposing Dirichlet boundary conditions. The actual Dirichlet boundary condition entering quatrex's' Poisson solver is not directly the voltage parameter \(V_{\mathrm{gate}}\) set for that surface. Instead, we also have to consider the electrostatic alignment of semiconductor channel and the metallic gate contact. The actual Dirichlet boundary condition entering the Poisson problem is

where \(\Phi_\mathrm{gate}\) is the metal's work_function, \(\chi_\mathrm{channel}\) is the semiconductor channel's electron_affinity, and \(\Delta E_{\mathrm{CB}-F} = E_\mathrm{CB} - E_\mathrm{F}\) is the distance between Fermi level and conduction band in the semiconductor.

Multifreedom Constraints for Periodic Boundaries#

Since gmsh supports the creation of periodic meshes, we can make use of multifreedom constraints (MFC)⁴ to enforce periodic boundary conditions straightforwardly. We construct a map from one surface to its image and then assemble an MFC transformation matrix that couples the image degrees of freedom to the authoritative ones.