10. Transport in the Coulomb blockade regime

In this section, we describe the basic equations implemented in QTCAD to model transport under strong quantum confinement (such as through a quantum dot). The most important effect that is captured by the QTCAD junction simulator is the Coulomb blockade. Through QTCAD transport simulations, it is possible to calculate Coulomb peaks and charge-stability diagrams, which are used experimentally to demonstrate the single-electron regime.

The junction concept

Fig. 10.1 illustrates the junction concept, which is implemented in the Junction class of the QTCAD transport simulator.

A junction consists of a Device or SubDevice object connected to two Lead objects, as illustrated schematically in real space in Fig. 10.1 (a). The device is always modeled explicitly in real space. In contrast, the leads may either be modeled phenomenologically using the Lead class, in which case no explicit geometric representation is considered, or explicitly using the WKBLead class, in which case the leads are explicitly associated with certain physical regions of a Device object.

The electrostatics and energy-level structure of the source, drain and dot may be calculated with the QTCAD Device simulator, accounting for full device geometry. Once stored inside a Junction object, these energy levels may be shifted by phenomenological gate potentials, accounting for lever arms that are either input from experimental data, or calculated from the device geometry using QTCAD’s Poisson and Schrödinger solvers (see Lever arm theory, Calculating the lever arm of a quantum dot, and Quantum transport—Charge stability diagram of a double quantum dot). This approach enables to save time because it avoids having to solve the Poisson and Schrödinger equations every time the gate potential is changed.

Fundamental tunneling Hamiltonian

In QTCAD’s transport solver, we consider a many-body description of the device, including Coulomb interactions between electrons. In addition, we consider that coupling between the device and the leads arises from tunneling events. We thus use the following Hamiltonian for the device and leads in the Fock basis corresponding to single-electron eigenfunctions of the source, drain, and leads:

(10.1)\[\begin{split}H &= \sum_{i\sigma} \epsilon_{i\sigma}\hat c^\dagger_{i\sigma}c_{i\sigma} + \frac 12\sum\limits_{ijkl,\sigma \sigma '} {{V_{ijlk}}} \hat c_{i\sigma }^\dagger \hat c_{j\sigma '}^\dagger {\hat c_{k\sigma '}}{\hat c_{l\sigma }}\\ &\qquad + \sum_{\mathbf kL}\epsilon_{\mathbf kL} c^\dagger_{\mathbf kL} c_{\mathbf kL} + \sum_{\mathbf kL,i\sigma}\left(t^\ast_{\mathbf kL,i\sigma} c^\dagger_{\mathbf kL} c_{i\sigma} + \textrm{H.c.}\right).\end{split}\]

In Eq. (10.1), the first line corresponds to the many-body Hamiltonian introduced in Eq. (8.5), in which \(\epsilon_{i\sigma}\) is the energy of the single-electron eigenstate with orbital index \(i\) and degeneracy index \(\sigma\) (arising from, e.g., spin or valley), and \(V_{ijlk}\) represents matrix elements of the Coulomb interaction defined in Eq. (8.6). The first term of the second line is the free Hamiltonian of the source and drain, with \(\epsilon_{\mathbf kL}\) the energy of a lead \(L\in \{\textrm{S}, \textrm{D}\}\) (source or drain) eigenstate with wave vector \(\mathbf k\). Finally, the last term of the equation captures tunneling between the device and the leads, with \(t_{\mathbf kL,i\sigma}\) being the tunneling matrix element between \((i\sigma)\) and \(\mathbf k\) single-electron eigenstates of the device and lead \(L\), respectively.

The sequential tunneling master equation

In the QTCAD transport simulator, the approach is to first diagonalize the device Hamiltonian using the exact diagonalization technique introduced in Many-body solver, and then treat coupling to the leads as a weak perturbation.

Within the sequential tunneling regime, we model transport through the device using the master equation [BF04]

(10.2)\[\dot p_m = -p_m\sum_{n\neq m}\gamma_{nm} +\sum_{n\neq m}p_n\gamma_{mn}=0,\]

where we have assumed steady-state conditions. In (10.2), we have introduced the indices \(n\) and \(m\), which run over all the many-body eigenstates of the device. In addition, we have introduced the occupation probability \(p_m\) of state \(m\) (note that occupation probabilities are subject to the condition \(\sum_m p_m=1\)) and the total transition rates

(10.3)\[\gamma_{mn}=\Gamma^{\textrm S}_{n\rightarrow m} + \Gamma^{\textrm D}_{n\rightarrow m}\]

from state \(n\) to state \(m\), where \(\Gamma^{\textrm S(\textrm D)}_{n\rightarrow m}\) is the contribution to the transition rate induced by the source (drain) lead.

The tunneling Hamiltonian introduced in Eq. (10.1) may cause two types of transitions between many-body eigenstates. Representing the \(\alpha\)-th many-body eigenstate within the \(N\)-electron subspace as \(|\alpha_N\rangle\), tunneling may lead to transitions of the form \(|\alpha_N\rangle\rightarrow |\beta_{N+1}\rangle\) (the device absorbs an electron from the lead) or \(|\alpha_N\rangle\rightarrow |\beta_{N-1}\rangle\) (the device releases an electron into the lead). The rates for these transitions are given by Fermi’s golden rule:

(10.4)\[\begin{split}{\Gamma ^L}({\alpha _N} \to {\beta _{N + 1}})\\ &\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!= {n_F}\left(\frac{{E_\beta^{N+1} } - {E_\alpha^N } - {\mu _L}}{k_B T}\right) \sum\limits_{\mu \nu } {\left\langle {{\beta _{N + 1}}} \right|} c_\mu ^\dagger \left| {{\alpha _N}} \right\rangle \Gamma _{\mu \nu }^L({E_\beta^{N+1} } - {E_\alpha^N })\left\langle {{\alpha _N}} \right|c_\nu ^{} \left| {{\beta _{N + 1}}} \right\rangle,\\ {\Gamma ^L}({\alpha _N} \to {\beta _{N - 1}})\\ &\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!= \left[ {1 - {n_F}\left(\frac{{E_\alpha^N } - {E_\beta^{N-1} } - {\mu _L}}{k_B T}\right)} \right] \sum\limits_{\mu \nu } {\left\langle {{\alpha _N}} \right|} c_\mu ^\dagger \left| {{\beta _{N - 1}}} \right\rangle \Gamma _{\mu \nu }^L({E_\alpha^N } - {E_\beta^{N-1} })\left\langle {{\beta _{N - 1}}} \right|c_\nu ^{} \left| {{\alpha _N}} \right\rangle,\end{split}\]

where \(L\;\in\;\{\textrm S, \textrm D\}\) is the lead index (source or drain), \(\mu\) and \(\nu\) are collective indices that simultaneously label the orbital-state and degeneracy indices of the single-electron basis states (i.e., \(\mu\to i,\sigma\); \(\nu \to j,\sigma'\)), \(E^N_\alpha\) is the energy of the many-body eigenstate \(|\alpha_N\rangle\), \(\mu_L\) is the chemical potential of the lead \(L\), and \(n_F(x)\) is the Fermi–Dirac distribution

(10.5)\[n_F(x) = \frac{1}{1+\textrm e^x}.\]

Importantly, in Eq. (10.4), we have introduced the broadening function

(10.6)\[\Gamma _{\mu \nu }^L(E) = \frac{2\pi}\hbar \sum\limits_{{\mathbf k_L}} {t_{{\mathbf k_L},\mu }^*t_{{\mathbf k_L},\nu}} \;\delta (E - {\epsilon_{{\mathbf k_L}}}),\]

where \(\mathbf k_L\) labels single-electron eigenstates of lead \(L\).

As will be seen below, there are two ways in which the broadening function may be evaluated in QTCAD:

through a featureless approximation,
using the WKB approximation.

Broadening function: featureless approximation

In the featureless approximation, the broadening function \(\Gamma^L_{\mu\nu}(E)\) is simply taken to be a constant independent of \(E\), \(\mu\), and \(\nu\) for each lead. By default, this constant is \(1\;\textrm{Hz}\) in QTCAD, but may be modified to arbitrary values through the set_source_broad_func and set_drain_broad_func methods of the Junction class.

See Calculating the charge-stability diagram of a quantum dot for an example of transport simulation in which these methods are used.

Broadening function: WKB calculation

The Wentzel–Kramers–Brillouin (WKB) approximation is a famous method used in the calculation of one-dimensional quantum tunneling problems.

To use the WKB approximation to calculate the broadening function, we make the following assumptions:

The Schrödinger equation of the relevant lead is separable in the \(x\), \(y\), and \(z\) directions: the lead potential is given by \(V_L(x,y,z) = V^L_x(x) + V^L_y(y) + V^L_z(z)\).
The source, device, and drain are aligned along a single Cartesian axis; in this theory section, we choose this axis to be the \(x\) axis.
Within the relevant lead, there is confinement along the \(y\) and \(z\) directions, leading to lead subbands associated with quantum numbers \(n_y\) and \(n_z\), respectively.
A single effective mass \(m^\ast\) can be defined in the transport direction.
A continuum (density-of-states) approximation may be taken for the lead eigenstates along \(x\).

Under these assumptions, the broadening function is given by

(10.7)\[\begin{split}\Gamma_{\mu\nu}^L(E)&=\frac{\ell}{\hbar^2} \sum_{n_y,n_z}\sqrt{\frac{2m^\ast}{E-E_{n_y}-E_{n_z}}}\; t^\ast_{n_y,n_z,\mu}(E)t_{n_y,n_z,\nu}(E)\theta_{n_yn_z}(E),\\ \theta_{n_yn_z}(E)&=\left\{ \begin{array}{ll} 1 &\textrm{if}\;E>E_{n_y}+E_{n_z},\\ 0 &\textrm{otherwise}. \end{array} \right.\end{split}\]

where \(E_{n_y}\) and \(E_{n_z}\) are the eigenenergies of the lead Schrödinger equations along the transverse directions \(y\) and \(z\), and \(\ell\) is the length of the lead in the \(x\) direction.

The tunneling matrix elements \(t_{n_y,n_z,\mu}(E)\) entering Eq. (10.7) are given by the transition matrix elements [Rei95]

(10.8)\[H_\textrm{LD} = \int_{x_1}^\infty dx \int_{-\infty}^\infty dy\int_{-\infty}^\infty dz \;\psi_{\textrm D}^\ast(\mathbf{r}) \left[V_{\textrm D}(\mathbf{r})- V_0 - V^L_y(y) - V^L_z(z)\right]\psi_L(\mathbf{r}),\]

where \(\psi_{\textrm D}(\mathbf{r})\) is a device single-particle envelope function and \(\psi_L(\mathbf{r}) = \psi^{WKB}_L(x)\psi^{y}_L(y)\psi^{z}_L(z)\) is a lead wave function (recall the separability assumption). As illustrated below, we have also introduced the barrier height \(V_0\), the quantum dot potential \(V_{\textrm D}(\mathbf{r})\), and the barrier position \(x_1\). In the figure below, we have also introduced \(V^{\textrm D}_x(x)\), a linecut of \(V_{\textrm D}(\mathbf{r})\) along the \(x\) direction.

Potentials involved in the transition matrix elements. — Fig. 10.2 Potentials involved in the calculation of the transition matrix elements.

The quantum dot wave function is easily calculated by numerically solving Schrödinger’s equation in the dot region (see, e.g., Simulating bound states in a quantum dot using the Poisson and Schrödinger solvers). In contrast, the lead wave function along \(x\) is obtained analytically using the WKB approximation

(10.9)\[\psi^{WKB}_L(x) = \frac{D}{\sqrt{P(x)}}\exp \left[-\frac1\hbar\int_{x_0}^x dx'P(x')\right],\]

where \(D\) is a normalization constant and

(10.10)\[P(x) = \sqrt{|2m^\ast[\epsilon - V^L_x(x)]|}.\]

In the above equations, we have introduced \(\epsilon = E - E_{n_y} - E_{n_z}\), the energy of the lead electron reserved for motion along \(x\). We have also introduced \(x_0\), the classical turning point of the lead potential.

Note that in Eq. (10.8), we have assumed that the quantum dot wave function is negligible for \(x<x_1\). This is a crucial approximation since in QTCAD, quantum dot wave functions are calculated within dot regions whose boundaries loosely coincide with barrier locations.

Note

In QTCAD, \(x_1\) is chosen automatically as the edge of the SubDevice representing the quantum system (SubDevice over which the Schrödinger equation is solved to find bound states) along the transport direction. For an example, see see Quantum transport—WKB approximation.

Note

The WKB feature is currently limited to single-valley electron calculations. In other words, it is currently not possible to use the WKB feature in combination with the Schrödinger solver for holes, or with the multi-valley effective-mass theory solver.

Sequential tunneling current

The charge current associated with electrons entering the device from the source is

(10.11)\[{I^\textrm S} = - e\sum\limits_{\alpha \beta N} {{p _{\alpha _N}} \left[ {{\Gamma^\textrm S}({\alpha _N} \to {\beta _{N + 1}}) - {\Gamma^\textrm S}({\alpha _N} \to {\beta _{N - 1}})} \right]}.\]

Also, introducing the charge current \(I^\textrm D\) associated with electrons entering the device from the drain, we have

(10.12)\[I^\textrm S + I^\textrm D = 0\;\Rightarrow\; I^\textrm D = -I^\textrm S\]

due to charge conservation.

Therefore, the current flowing through the device can easily be calculated from the occupation probabilities (obtained by solving the master equation) and the transition rates introduced above.

Particle addition spectrum

When the source–drain bias is sufficiently small, the quantum dot system is approximately at thermodynamic equilibrium with the leads. Since we are ultimately interested in the charge stability diagram of the quantum dot system, a central quantity of interest is the total particle number \(N=\sum_i n_i\), with \(n_i\) the number of particles in orbital \(i\) introduced in Eq. (8.16). At thermal equilibrium with a heat bath at constant temperature \(T\) and chemical equilibrium with leads at chemical potential \(\mu\), the average particle number is

(10.13)\[\langle N \rangle = k_B T\frac{\partial \log \mathcal Z}{\partial \mu},\]

where we have introduced the grand partition function

(10.14)\[\mathcal Z = \mathrm{Tr}\left[ \exp\left( -\frac{H-\mu N}{k_B T} \right) \right].\]

Transitions in the total particle number are associated with peaks in the particle addition spectrum defined by [HFJ+17]

(10.15)\[\frac{\partial \langle N\rangle}{\partial \mu} = \frac{\langle N^2\rangle - \langle N\rangle^2}{k_B T}.\]

Computing this particle addition spectrum as a function of gate bias configurations gives diagrams in which non-zero values indicate configurations that allow the total particle number to change as a function of source and drain chemical potentials, and thus for the steady-state current flowing through the quantum dot system to change. In many cases, this can allow to identify transitions in charge stability diagrams more efficiently (see, e.g., the Quantum transport—Charge stability diagram of a double quantum dot tutorial).