4. Calculating the Coulomb peaks of a quantum dot

4.1. Requirements

4.1.1. Software components

QTCAD

4.1.2. Mesh file

qtcad/examples/practical_application/meshes/gated_dot.msh2

4.1.3. Python script

qtcad/examples/practical_application/4-many_body.py

4.1.4. References

4.2. Briefing

In this previous tutorial, we computed the lever arm of a quantum dot. Transport in a quantum dot is characterized by a phenomenon known as Coulomb blockade. Specifically, only an integer number of electrons can reside in the quantum dot. When a chemical potential of the dot lies between the source and drain chemical potentials, the number of electrons in the dot can fluctuate between \(N-1\) and \(N\) (where \(N\) is an integer), resulting in finite conductance. On the other hand, when there are no chemical potentials between the source and drain chemical potentials, the number of electrons cannot fluctuate between two integer values, resulting is zero conductance. These chemical potentials—i.e., the difference in total energy between a quantum dot with \(N-1\) and \(N\) electrons in their ground state (see Chemical potentials)—therefore play a crucial role in transport in this device. The Coulomb peaks, i.e. local maximums in the dot conductivity with respect to an applied gate potential, are related to the chemical potentials and the lever arm (see Coulomb peak positions).

In this tutorial, we will compute the chemical potentials and Coulomb peaks of the quantum dot of the previous tutorials at a value of the plunger gates potential of \(0.3\textrm{ V}\), which approximately corresponds to the gate bias at which the quantum-dot ground state energy (which is equal to the chemical potential with \(N=1\)) crosses the Fermi level (set to zero by convention in QTCAD).

4.3. Setting up the device

4.3.1. Header

We start by importing relevant modules.

import pathlib
import numpy as np

from device.mesh3d import Mesh, SubMesh
from device import io
from device import constants as ct

from device import materials as mt
from device import Device, SubDevice
from device.schrodinger import Solver as SingleBodySolver
from device.many_body import Solver as ManyBodySolver

These modules were described in the previous tutorials with the exception of many_body, which is used to define and solve the many-body problem of \(N\) interacting electrons in a device.

4.3.2. Defining the device and its electric potential

We start by defining the device as was previously done. Note that the plunger gates potentials, Vtop_2 and Vbottom_2, are set to 0.3.

# Load the mesh
script_dir = pathlib.Path(__file__).parent.resolve()
path_mesh = str(script_dir /'meshes'/'gated_dot.msh2')
scaling = 1e-6
mesh = Mesh(scaling, path_mesh)

# Create device
d = Device(mesh)
d.set_temperature(0.1)

# # Applied potentials
Vtop_1 = -0.5
Vtop_2 = 0.3
Vtop_3 = -0.5
Vbottom_1 = -0.5
Vbottom_2 = 0.3
Vbottom_3 = -0.5

# Work function of the metallic gates at midgap of GaAs
barrier = 0.834*ct.e    # n-type Schottky barrier height

# Ionized dopant density
doping = 3e18*1e6   # In SI units

# Set up materials in heterostructure stack
d.new_region("cap", mt.GaAs)
d.new_region("barrier", mt.AlGaAs)
d.new_region("dopant_layer", mt.AlGaAs,ndoping=doping)
d.new_region("spacer_layer", mt.AlGaAs)
d.new_region("spacer_dot", mt.AlGaAs)
d.new_region("two_deg", mt.GaAs)
d.new_region("two_deg_dot", mt.GaAs)
d.new_region("substrate", mt.GaAs)
d.new_region("substrate_dot", mt.GaAs)

# Align the bands across hetero-interfaces using GaAs as the reference material
d.align_bands(mt.GaAs)

# Set up boundary conditions
d.new_schottky_bnd("top_gate_1", Vtop_1, barrier)
d.new_schottky_bnd("top_gate_2", Vtop_2, barrier)
d.new_schottky_bnd("top_gate_3", Vtop_3, barrier)
d.new_schottky_bnd("bottom_gate_1", Vbottom_1, barrier)
d.new_schottky_bnd("bottom_gate_2", Vbottom_2, barrier)
d.new_schottky_bnd("bottom_gate_3", Vbottom_3, barrier)
d.new_ohmic_bnd("ohmic_bnd")

The next step is to set the electric potential in the device. This could be done by solving the nonlinear Poisson equation, as has been done in previous tutorials. However, here, we simply load this potential from the relevant .hdf5 file generated in Calculating the lever arm of a quantum dot, which is set as the device potential saved in d.phi. Remark that the method set_potential also automatically calls set_V_from_phi to update the confinement potential accordingly.

# Load the confinement potential from previous run
path_in = str(script_dir/"output"/"electric_potential_0.300V.hdf5")
d.set_potential(io.load(path_in))

4.4. Solving Schrödinger’s equation

Next, Schrödinger’s equation is solved using the usual method.

# Create a submesh in which to solve Schrodinger's equation
submesh = SubMesh(mesh, ["spacer_dot","two_deg_dot","substrate_dot"])

# Create subdevice
subdevice = SubDevice(d, submesh)

# Solve the single-electron problem to get basis set
slv = SingleBodySolver(subdevice)
slv.solve()
subdevice.print_energies()

In this case, the energy levels are as follows.

Energy levels (eV)
[0.00048222 0.00666745 0.01287914 0.01373345 0.01921634 0.02014803
 0.02531556 0.02690116 0.02754001 0.03198617]

As a sanity check, we visualize the confining potential and the eigenfunctions of the ground state and the first two excited states.

# Plot the first few envelope functions
from device import analysis as an
an.plot_slices(submesh,subdevice.V/ct.e,title="Confinement potential")
an.plot_slices(submesh,subdevice.eigenfunctions[:,0],title="Ground state")
an.plot_slices(submesh,subdevice.eigenfunctions[:,1],title="1st excited state")
an.plot_slices(submesh,subdevice.eigenfunctions[:,2],title="2nd excited state")

This results in the following figures.

4.5. Setting up the many-body solver and calculating the Coulomb peaks

We are now ready to define a many-body Solver object, and run it with:

# Create the many-body solver and call its solve method
num_states = 3              # Number orbital basis states to consider
n_degen = 2                 # spin degenerate system
slv = ManyBodySolver(subdevice, num_states=num_states, n_degen=n_degen,
    alpha=0.036)
slv.solve()

The many-body Solver object constructor has several arguments. The first, subdevice, corresponds to the subdevice object of the quantum dot; this argument is not optional. Here, we supply three additional optional arguments. The first, num_states = num_states specifies that only the three lowest energy levels in the quantum dot should be considered. Solving the many-body problem involves the calculation of numerous Coulomb integrals; the number of these integrals scales with the fourth power of the number of energy levels. The second, n_degen = n_degen, describes the degeneracy of the dot’s eigenstates. Since our system is spin-degenerate, this argument is equal to \(2\). The third, alpha=0.036*ct.e is the lever arm which we computed in the previous tutorial.

At this point, the calculation of Coulomb integrals begins; these integrals describe the Coulomb interactions between electrons in the dot.

Computing Coulomb integrals ███∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙ 12%

This calculation may take up to a few minutes to run, after which the resulting many-body Hamiltonian is diagonalized and the results are stored in the subdevice object. After this, we can print and save the chemical potentials and Coulomb peaks.

# Output positions of chemical potentials and Coulomb peaks
path_chem_potentials = str(script_dir/"output"/"chem_potentials.txt")
path_peak_pos = str(script_dir/"output"/"peak_positions.txt")
print('chemical potentials (eV):', subdevice.chem_potentials/ct.e)
np.savetxt(path_chem_potentials,subdevice.chem_potentials/ct.e)
print('positions of Coulomb peaks relative to input value (eV):',
    subdevice.coulomb_peak_pos)
np.savetxt(path_peak_pos,subdevice.coulomb_peak_pos)

The chem_potentials attribute of the subdevice corresponds to the list of chemical potentials obtained by running the many-body solver, while the coulomb_peak_pos attribute corresponds to the Coulomb peaks. This should result in the following.

chemical potentials (eV): [0.00048222 0.00759687 0.01678584 0.02256699 0.03193739 0.03735708]
positions of Coulomb peaks relative to input value (eV): [0.01339498 0.21102411 0.46627331 0.62686094 0.88714976 1.03769655]

4.6. Full code

__copyright__ = "Copyright 2022, Nanoacademic Technologies Inc."

import pathlib
import numpy as np
from device.mesh3d import Mesh, SubMesh
from device import io
from device import constants as ct
from device import materials as mt
from device import Device, SubDevice
from device.schrodinger import Solver as SingleBodySolver
from device.many_body import Solver as ManyBodySolver

# Load the mesh
script_dir = pathlib.Path(__file__).parent.resolve()
path_mesh = str(script_dir /'meshes'/'gated_dot.msh2')
scaling = 1e-6
mesh = Mesh(scaling, path_mesh)

# Create device
d = Device(mesh)
d.set_temperature(0.1)

# # Applied potentials
Vtop_1 = -0.5
Vtop_2 = 0.3
Vtop_3 = -0.5
Vbottom_1 = -0.5
Vbottom_2 = 0.3
Vbottom_3 = -0.5

# Work function of the metallic gates at midgap of GaAs
barrier = 0.834*ct.e    # n-type Schottky barrier height

# Ionized dopant density
doping = 3e18*1e6   # In SI units

# Set up materials in heterostructure stack
d.new_region("cap", mt.GaAs)
d.new_region("barrier", mt.AlGaAs)
d.new_region("dopant_layer", mt.AlGaAs,ndoping=doping)
d.new_region("spacer_layer", mt.AlGaAs)
d.new_region("spacer_dot", mt.AlGaAs)
d.new_region("two_deg", mt.GaAs)
d.new_region("two_deg_dot", mt.GaAs)
d.new_region("substrate", mt.GaAs)
d.new_region("substrate_dot", mt.GaAs)

# Align the bands across hetero-interfaces using GaAs as the reference material
d.align_bands(mt.GaAs)

# Set up boundary conditions
d.new_schottky_bnd("top_gate_1", Vtop_1, barrier)
d.new_schottky_bnd("top_gate_2", Vtop_2, barrier)
d.new_schottky_bnd("top_gate_3", Vtop_3, barrier)
d.new_schottky_bnd("bottom_gate_1", Vbottom_1, barrier)
d.new_schottky_bnd("bottom_gate_2", Vbottom_2, barrier)
d.new_schottky_bnd("bottom_gate_3", Vbottom_3, barrier)
d.new_ohmic_bnd("ohmic_bnd")

# Load the confinement potential from previous run
path_in = str(script_dir/"output"/"electric_potential_0.300V.hdf5")
d.set_potential(io.load(path_in))

# Create a submesh in which to solve Schrodinger's equation
submesh = SubMesh(mesh, ["spacer_dot","two_deg_dot","substrate_dot"])

# Create subdevice
subdevice = SubDevice(d, submesh)

# Solve the single-electron problem to get basis set
slv = SingleBodySolver(subdevice)
slv.solve()
subdevice.print_energies()

# Plot the first few envelope functions
from device import analysis as an
an.plot_slices(submesh,subdevice.V/ct.e,title="Confinement potential")	
an.plot_slices(submesh,subdevice.eigenfunctions[:,0],title="Ground state")	
an.plot_slices(submesh,subdevice.eigenfunctions[:,1],title="1st excited state")	
an.plot_slices(submesh,subdevice.eigenfunctions[:,2],title="2nd excited state")

# Create the many-body solver and call its solve method
num_states = 3              # Number orbital basis states to consider
n_degen = 2                 # spin degenerate system
slv = ManyBodySolver(subdevice, num_states=num_states, n_degen=n_degen,
    alpha=0.036)
slv.solve()

# Output positions of chemical potentials and Coulomb peaks
path_chem_potentials = str(script_dir/"output"/"chem_potentials.txt")
path_peak_pos = str(script_dir/"output"/"peak_positions.txt")
print('chemical potentials (eV):', subdevice.chem_potentials/ct.e)
np.savetxt(path_chem_potentials,subdevice.chem_potentials/ct.e)
print('positions of Coulomb peaks relative to input value (eV):',
    subdevice.coulomb_peak_pos)
np.savetxt(path_peak_pos,subdevice.coulomb_peak_pos)