1. Quantum transport—Master equation
1.1. Requirements
1.1.1. Software components
QTCAD
1.1.2. Python script
qtcad/examples/tutorials/nanowire_diamond.py
1.1.3. References
1.2. Briefing
In this tutorial, we demonstrate how to perform a quantum transport calculation. In this calculation, we solve the master equation that describes transport from a source to a drain due to sequential tunneling through a confined (quantum dot) region. The eigenstates of the quantum dot region that appear in the master equation are obtained by exact diagonalization of the many-body Hamiltonian of conduction electrons in the truncated basis formed by a few single-particle eigenstates. These single-particle eigenstates are found by solving the Schrödinger equation in the dot region. This method yields the current flowing through the device as a function of applied source–drain and gate potentials, and thus gives us access to the position of Coulomb peaks, allowing us to draw charge-stability diagrams.
1.3. Setting up the device
1.3.1. Header
In addition to what has already been invoked in
Poisson and Schrödinger simulation of a nanowire quantum dot, here we also need to import relevant classes and
functions from the
mastereq
and junction
modules.
import numpy as np
from matplotlib import pyplot as plt
import pathlib
from progress.bar import ChargingBar as Bar
from qtcad.device import constants as ct
from qtcad.device import materials as mt
from qtcad.device.mesh3d import Mesh
from qtcad.device import io
from qtcad.device import Device
from qtcad.device.schrodinger import Solver
from qtcad.device.many_body import SolverParams
from qtcad.transport.mastereq import seq_tunnel_curr
from qtcad.transport.junction import Junction
Junction
: Class that contains all information about the transport system, i.e. the two leads (source and drain), and the quantum-dot system.seq_tunnel_curr
: Function to compute transport properties
1.3.2. Load data
Here we assume that a device structure has already been treated, and
that a single-particle confinement potential was found and saved on our
disc. For conduction electrons, this corresponds to the conduction band
edge. In this tutorial, we reuse the result from Poisson and Schrödinger simulation of a nanowire quantum dot. Indeed, the conduction band edge
(along with several other variables) was saved in .hdf5
format.
We will thus load the mesh, create an empty device, and load the only quantity that we need (the conduction band edge) for the quantum transport calculation.
# Load the mesh
script_dir = pathlib.Path(__file__).parent.resolve()
path_mesh = script_dir / "meshes/" / "nanowire.msh"
path_in = script_dir / "output/" / "nanowire.hdf5"
mesh = Mesh(1e-9, str(path_mesh))
d = Device(mesh=mesh, conf_carriers="e")
d.new_region('channel_ox',mt.SiO2)
d.new_region('source_ox',mt.SiO2)
d.new_region('drain_ox',mt.SiO2)
d.new_region('channel', mt.Si)
d.new_region('source', mt.Si, pdoping=5e20*1e6, ndoping=0)
d.new_region('drain', mt.Si, pdoping=5e20*1e6, ndoping=0)
d.V = io.load(path_in, var_name="EC") * ct.e
The last line sets the electron potential energy to be the conduction band edge calculated in the previous tutorial. Note that we multiply by the fundamental electric charge to obtain SI units, since energies obtained in the previous tutorial were saved using \(\textrm{eV}\) units.
1.3.3. Solving Schrödinger’s equation for single electrons
The first step is to solve the single-body Schrödinger’s equation for
the loaded potential energy, stored in the device attribute d.V
. The
resulting envelope functions will then be used as a basis set to solve
the many-body problem, including Coulomb interaction between electrons.
# Solve Schrodinger's equation for single electrons
s = Solver(d)
s.solve()
1.3.4. Initializing a Junction object
A Junction
object is initialized through the following interface
# Many-body solver parameters
solver_params = SolverParams()
solver_params.num_states=2 # Number of basis states to keep
solver_params.n_degen=2 # spin degenerate system
# Initialize junction
jc = Junction(d, many_body_solver_params=solver_params)
It is recommended to set the parameters of the many-body solver that lives
inside the Junction
object
using a many-body SolverParams
object that is passed as an optional argument to the
Junction
constructor. For example, parameters that are often tuned this way are:
num_states: Number of one-particle eigenstates to be included in the many-body basis set. The Fock-space dimension in the following many-body simulation is thus 2^(num_states*spin)
n_degen: The degeneracy of each level. In this case we set
n_degen = 2
to account for the spin degeneracy.
The initialization procedure may take a while since the program will set up the many-body Hamiltonian and diagonalize it. The many-body eigenenergies can be extracted via
jc.get_all_eig_energies()
After initialization one may still be able to set certain parameters through Junction-class methods, such as
jc.set_temperature(0.1) # set temperature in unit of K
jc.setVg(0.1) # set gate voltage in volts
jc.setVs(0.1) # set source voltage in volts
jc.setVd(0.1) # set drain voltage in volts
The source–drain voltages are assumed to rigidly shift the chemical potentials of corresponding reservoirs, and applying a gate voltage effectively shifts the entire energy spectrum of the device by \(-e V_g\), corresponding to an ideal lever arm \(\alpha = 1\).
Note that the electrostatics is not recalculated when we manually set
these potential parameters. However, the actual lever arm for a given
device geometry may be accounted for by setting the keyword argument
alpha
of the Junction
constructor to the chosen value (different
from the default value of 1). The single-particle eigenenergies are then
shifted by \(-e \alpha V_g\) when applying a gate potential \(V_g\).
The lever arm \(\alpha\) may easily be found by finding the single-particle
eigenenergies as a function of \(e V_g\), and fitting the result to a
straight line, as see in Calculating the lever arm of a quantum dot.
The slope of this line then gives \(-\alpha\).
Note
As explained in Many-body analysis of a nanowire quantum dot, the gate bias is expressed
with respect to the reference gate potential, i.e., the potential applied
on the gate boundary condition when solving Poisson’s equation which lead
to the confinement potential used here. For example, if the reference gate
potential was \(V_\mathrm{GS}=0.5\;\mathrm V\), the bias of
\(0.1\;\mathrm V\) applied here through the
setVg
method corresponds to
a physical gate potential of
\(0.5\;\mathrm V+0.1\;\mathrm V=0.6\;\mathrm V\).
This approach enables to produce charge stability diagrams even for
confinement potentials that were loaded from previous simulations, or
from analytic confinement potentials that are not the numerical solution
of the electrostatics of a realistic device.
1.4. Transport calculation
To calculate transport currents and statistical distribution we invoke
the seq_tunnel_curr
function
Il,Ir,prob = seq_tunnel_curr(jc)
Il: charge current from source
Ir: charge current from drain. The relationship
Il
+Ir
=0 should be respected.prob: occupation probability for each many-body state in the stationary (steady-state) regime.
We note that the currents Il
and Ir
are evaluated assuming that the
single-particle broadening functions appearing in the master equation rates
are constant over all possible pairs of single-particle eigenstates (see
The sequential tunneling master equation and Broadening function: featureless approximation sections of our
theory documentation). This has the benefit of speeding up calculations
which involve solving master equations, while giving a qualitative
understanding of different transport phenomena, e.g. Coulomb Blockade.
However, the currents themselves are not meant to be quantitatively
accurate. In this tutorial we wish to produce Coulomb diamonds.
To do this we sweep the source/drain and gate voltages over a range.
jc.set_temperature(1) # set temperature in unit of K
# Get Coulomb diamonds
v_gate_rng = np.linspace(0.45, 0.85, num=100)
v_source_rng = np.linspace(-0.08, 0.08, num=100)
diam = np.zeros((v_gate_rng.size, v_source_rng.size), dtype=float)
number_grid_points = v_gate_rng.size * v_source_rng.size # number of sampled grid points
progress_bar = Bar("Computing charge stability diagram", max = number_grid_points) # initialize progress bar
# loop over gate voltage
for i, v_gate in enumerate(v_gate_rng):
# Set gate voltage
jc.setVg(v_gate)
# loop over source/drain voltage
for j, v_source in enumerate(v_source_rng):
# Set source and drain voltages
jc.setVs(v_source)
jc.setVd(-v_source) # Drain voltage opposite in sign to source
# Compute currents
Il,Ir,prob = seq_tunnel_curr(jc) # transport calculation
diam[i,j] = Il
progress_bar.next()
progress_bar.finish()
This procedure yields a 2d array diam
which stores the current
obtained under every sampled voltage configuration. While the first
index of the array corresponds to gate potentials, the second index
corresponds to source–drain biases. Note that, to keep track of the progress
of the calculation of this charge-stability diagram, we have instantiated a
progress bar progress_bar
as implemented in the progress.bar
Python
library.
We get a quantity proportional to the differential conductance (derivative of the current with respect to source–drain bias) by taking differences of the current
# differential conductance
diam_diff_cond = np.abs(diam[:,1:] - diam[:,:-1])
Note
The convention for the sign of the currents is the following: a positive current means conventional (positive) charge tunneling from a lead to the dot. Such a current is expressed in Ampere units. Note that within this convention, when a net current is flowing through the device, individual currents associated with each lead will have opposite signs.
1.4.1. Plotting results
Finally, we plot diam_diff_cond
to visualize the Coulomb diamonds.
# plot results
fig, axs = plt.subplots()
axs.set_ylabel('$V_s(V)$',fontsize=16)
axs.set_xlabel('$V_g(V)$',fontsize=16)
diff_cond_map = axs.imshow(np.transpose(diam_diff_cond), interpolation='bilinear',
extent=[v_gate_rng[0], v_gate_rng[-1],
v_source_rng[-2], v_source_rng[0]], aspect="auto")
fig.colorbar(diff_cond_map, ax=axs, label="Differential conductance (arb. units)")
fig.tight_layout()
plt.show()
The diamonds in the above diagram manifest the Coulomb blockade regions in which the source–drain and gate voltage configuration determines a quantized electron occupation number due to Coulomb repulsion. Recall that the gate bias expressed on the \(x\) axis of this plot should be understood as given with respect to the reference gate potential applied at the boundary when solving Poisson’s equation. Therefore, for the confinement potential used here, which was loaded from the Poisson and Schrödinger simulation of a nanowire quantum dot tutorial, an offset of \(V_\mathrm{GS}=0.5\;\mathrm V\) should be added to the \(x\) axis to plot the charge stability diagram as a function of the true applied potential at the gate.
In the next tutorial (Quantum transport—WKB approximation), we will perform a more involved transport calculation using the WKB approximation. We expect this more involved calculation to be more quantitatively accurate than the procedure presented here.
1.5. Full code
__copyright__ = "Copyright 2022-2024, Nanoacademic Technologies Inc."
import numpy as np
from matplotlib import pyplot as plt
import pathlib
from progress.bar import ChargingBar as Bar
from qtcad.device import constants as ct
from qtcad.device import materials as mt
from qtcad.device.mesh3d import Mesh
from qtcad.device import io
from qtcad.device import Device
from qtcad.device.schrodinger import Solver
from qtcad.device.many_body import SolverParams
from qtcad.transport.mastereq import seq_tunnel_curr
from qtcad.transport.junction import Junction
# Load the mesh
script_dir = pathlib.Path(__file__).parent.resolve()
path_mesh = script_dir / "meshes/" / "nanowire.msh"
path_in = script_dir / "output/" / "nanowire.hdf5"
mesh = Mesh(1e-9, str(path_mesh))
d = Device(mesh=mesh, conf_carriers="e")
d.new_region('channel_ox',mt.SiO2)
d.new_region('source_ox',mt.SiO2)
d.new_region('drain_ox',mt.SiO2)
d.new_region('channel', mt.Si)
d.new_region('source', mt.Si, pdoping=5e20*1e6, ndoping=0)
d.new_region('drain', mt.Si, pdoping=5e20*1e6, ndoping=0)
d.V = io.load(path_in, var_name="EC") * ct.e
# Solve Schrodinger's equation for single electrons
s = Solver(d)
s.solve()
# Many-body solver parameters
solver_params = SolverParams()
solver_params.num_states=2 # Number of basis states to keep
solver_params.n_degen=2 # spin degenerate system
# Initialize junction
jc = Junction(d, many_body_solver_params=solver_params)
jc.set_temperature(1) # set temperature in unit of K
# Get Coulomb diamonds
v_gate_rng = np.linspace(0.45, 0.85, num=100)
v_source_rng = np.linspace(-0.08, 0.08, num=100)
diam = np.zeros((v_gate_rng.size, v_source_rng.size), dtype=float)
number_grid_points = v_gate_rng.size * v_source_rng.size # number of sampled grid points
progress_bar = Bar("Computing charge stability diagram", max = number_grid_points) # initialize progress bar
# loop over gate voltage
for i, v_gate in enumerate(v_gate_rng):
# Set gate voltage
jc.setVg(v_gate)
# loop over source/drain voltage
for j, v_source in enumerate(v_source_rng):
# Set source and drain voltages
jc.setVs(v_source)
jc.setVd(-v_source) # Drain voltage opposite in sign to source
# Compute currents
Il,Ir,prob = seq_tunnel_curr(jc) # transport calculation
diam[i,j] = Il
progress_bar.next()
progress_bar.finish()
# differential conductance
diam_diff_cond = np.abs(diam[:,1:] - diam[:,:-1])
# plot results
fig, axs = plt.subplots()
axs.set_ylabel('$V_s(V)$',fontsize=16)
axs.set_xlabel('$V_g(V)$',fontsize=16)
diff_cond_map = axs.imshow(np.transpose(diam_diff_cond), interpolation='bilinear',
extent=[v_gate_rng[0], v_gate_rng[-1],
v_source_rng[-2], v_source_rng[0]], aspect="auto")
fig.colorbar(diff_cond_map, ax=axs, label="Differential conductance (arb. units)")
fig.tight_layout()
plt.show()