5. Simulating the Coulomb peak positions of an SET

5.1. Requirements

5.1.1. Software components

  • QTCAD

5.1.2. Python script

  • qtcad/examples/practical_application/FDSOI/3-coulomb_peaks.py

5.1.3. References

5.2. Briefing

In spin-qubit devices, quantum dots operating as single-electron transistors (SETs) can be used as sensitive charge sensors to read out the state of other nearby quantum dots [HKP+07]. Because SETs display discrete energy levels, in the low-bias regime, the current through a SET exhibits (Coulomb) peaks as a function of the gate voltage. These peaks occur when the chemical potential \(\mu(N)\) for a given electron number \(N\) on the SET aligns with the Fermi level (\(E_F=0 \, \mathrm{eV}\)) of an adjacent reservoir (the source in this particular example), allowing electrons to tunnel onto and off the SET island. Because the chemical potential of the SET depends on the gate voltage, the Coulomb peaks can be detected by sweeping the gate voltage and measuring the current through the SET. More details about transport through a quantum dot can be found in Quantum transport—Master equation.

The goal of this part of the practical application is to calculate the positions of the Coulomb peaks of the SET created in the previous sections, Mesh generation using QTCAD Builder and Creating the FD-SOI device. In particular, we will demonstrate how to compute the chemical potential of the SET as a function of the gate voltage using the Schrödinger–Poisson Solver introduced in Computing the chemical potential of an SET device. We will also extract certain relevant features of the SET from the computed chemical potential data, such as the lever arm, addition energies, charging energy, and gate capacitance. The workflow presented here was inspired by Ref. [BSA00].

5.3. Setup

5.3.2. File paths

Next, we define the paths to input and output files for the solver and its results:

# Files -----------------------------------------------------------------------
script_dir           = pathlib.Path(__file__).parent.resolve()

path_out             = script_dir / "output"
path_local_out       = path_out / "coulomb_peaks"
phi_in_file          = path_out / "phi" / f"phi.hdf5"

leverarm_file        = str(path_local_out / "lever_arm.txt")
lever_arm_plot_file  = str(path_local_out / "lever_arm_plot.png")

5.3.3. Compute chemical potential of the FD-SOI SET

Now, we leverage the get_double_dot_fdsoi function defined in double_dot_fdsoi.py to create the FD-SOI device over the mesh generated in Mesh generation using QTCAD Builder. We then load the potential computed in Simulating electrostatics using the linear Poisson solver and set it as the electrostatic potential of d using the set_potential method of the Device class. We can then create a sub-device for the SET region by first creating a SubMesh, set_mesh, and using this submesh to define the SubDevice, set_device. Once that is done, we use the set_insulator_boundaries method to set insulator boundary conditions on all the external surfaces of the set_device. These boundary conditions force the wavefunctions to vanish at the boundaries when solving the Schrödinger equation.

# Setup ----------------------------------------------------------------------

# Create the device
d, set_region, qd_region = get_double_dot_fdsoi()
# load the potential
d.set_potential(io.load(phi_in_file, var_name="phi"))
# Create sub-device for the SET region
set_mesh = SubMesh(d.mesh, set_region)
set_device = SubDevice(d, set_mesh)
set_device.set_insulator_boundaries()

Next, we compute the chemical potential of the SET device as a function of the gate voltage using the compute_chemical_potential function defined in chemical_potential.py.

# V_set values
V_SET = np.array([V_set, V_set - 0.05])
# Number of particles
Nmin = 12
Nmax = 14
# Generate data
Mu = compute_chemical_potential(V_SET, Nmin, Nmax, d, set_device, out_path=path_local_out)

This computation is performed for two different values of the gate voltage controlling the SET device, V_set and V_set - 0.05 V, and for a number of particles ranging from Nmin = 12 to Nmax = 14. The results are stored in the array Mu, where each row corresponds to a different gate voltage and each column corresponds to a different number of particles. The values for V_SET, Nmin, and Nmax can be adjusted by the user to generate additional data points.

The compute_chemical_potential function also saves intermediate results to the path_local_out folder. For example, a slice the charge density and ground-state wavefunction for each gate voltage and number of particles. Below we show an example of such slices for V_SET = V_set and N = 13 electrons.

Charge density slice

Fig. 5.3.2 Slice of the charge density for V_SET = 1.25 V ("plunger_gate_1_bnd" boundary) and N = 13 electrons.

Ground state wavefunction slice

Fig. 5.3.3 Slice of the ground state wavefunction for V_SET = 1.25 V ("plunger_gate_1_bnd" boundary) and N = 13 electrons.

5.4. Data analysis

Now that the chemical potential data has been generated, we can analyze it to extract useful information about the SET device.

5.4.1. Coulomb peak positions

The lever arm is a measure of the efficiency of the gate voltage in modulating the chemical potentials of the quantum dot. We can compute the lever arm by calculating the slope of the chemical potential as a function of the gate voltage. This can be done using a simple linear fit between the two gate voltage points used in the simulation.

# Coulomb peak positions ------------------------------------------------------

# Fit a line: y = alpha*x + b
lever_arm = np.zeros(Mu.shape[1])
intercept = np.zeros(Mu.shape[1])
for n in range(Mu.shape[1]):
    alpha, b = np.polyfit(V_SET, Mu[:, n] / ct.e, 1)
    lever_arm[n] = alpha
    intercept[n] = b
    print(f"n={n + Nmin + 1}: Lever Arm [eV/V] = {alpha}, Intercept [eV] = {b}")

Note

The lever arm computed here is different than the one produced by the lever arm Solver. The latter computes the lever arm based on the shift in single-particle energy levels as a function of gate voltage, whereas here we compute the lever arm based on the shift in chemical potentials.

Once we have the lever arm, we can compute the positions of the Coulomb peaks by finding the gate voltage values for which the chemical potential aligns with the Fermi level, \(E_F = 0\). We can also compute the distance between consecutive Coulomb peaks.

# Compute Coulomb peak positions
# (x-intercept)
V_peaks = -intercept / lever_arm
print("Coulomb Peak Positions [V]:")
print(V_peaks)
# (Distance between peaks)
Delta_V_peaks = np.diff(V_peaks)
print("Distance between Coulomb Peaks [V]:")
print(Delta_V_peaks)

# Save
header = "Lever Arm [eV/V]    Intercept [eV]    Coulomb Peak Position [V]"
np.savetxt(leverarm_file, np.column_stack((lever_arm, intercept, V_peaks)), fmt="%.10f", header=header)

The above code prints

n=13: Lever Arm [eV/V] = -0.552741875999999, Intercept [eV] = 0.687271847199999
n=14: Lever Arm [eV/V] = -0.5493631299999989, Intercept [eV] = 0.7003871215999989
Coulomb Peak Positions [V]:
[1.24338661 1.2749074 ]
Distance between Coulomb Peaks [V]:
[0.03152079]

We can then plot the chemical potential data along with the linear fits to visualize the lever arm.

# Plot
plt.figure(figsize=(8, 5))
# (Add data to plot)
DeltaV = np.max(np.insert(V_SET, 0, np.max(V_peaks))) - np.min(np.insert(V_SET, 0, np.min(V_peaks)))
V_line = np.linspace(min(V_SET) - DeltaV, max(V_SET) + DeltaV, 200)
plt.axhline(y=0, color="black", linewidth=3)
for j, alpha in enumerate(lever_arm):
    b = intercept[j] # Intercept for given value of N
    mu_line = alpha * V_line + b
    plt.plot(V_line, mu_line, color="red", linewidth=2.0, linestyle="--")
    plt.axvline(V_peaks[j], color="black", linestyle="--", linewidth=2.)
for j in range(Mu.shape[1]):
    plt.scatter(V_SET, Mu[:, j] / ct.e, marker="o", s=60, label=f"N = {j + Nmin + 1}")

# (Labels, legend, grid)
plt.xlabel("$V_{SET}$ [V]", fontsize=16)
plt.ylabel("$\\mu$ [eV]", fontsize=16)
plt.legend(fontsize=14)
plt.grid(True)
# (Save)
plt.tight_layout()
plt.savefig(lever_arm_plot_file, dpi=300)

The resulting plot is shown below.

Lever arm plot.

Fig. 5.4.8 The chemical potential of the SET as a function of the gate voltage. The red dashed lines represent the linear fits used to compute the lever arm. The black dashed vertical lines indicate the positions of the Coulomb peaks (intercept of the linear fits with the x axis).

Note

The linear fits were produced using only two data points for each number of particles. For more accurate results, it is recommended to compute the chemical potential at multiple gate voltage values and perform a linear regression over all the data points. The user is free to adjust Nmin, Nmax, and V_SET to generate more data points and improve the accuracy of the lever arm calculation.

5.4.2. Addition Energies

The energy required to add an extra electron to the SET island is known as the addition energy. This energy can be computed as the difference in chemical potentials for adding the N-th and (N-1)-th electrons to the island [HKP+07]:

\[E_{add}(N) = \mu(N+1) - \mu(N).\]

Using the chemical potential data computed earlier, we can easily calculate the addition energies for the SET device:

# Compute addition energies ---------------------------------------------------
addition_energies = Mu[:, 1:] - Mu[:, :-1]
print("Addition Energies [eV]:")
print(addition_energies / ct.e)

which prints

Addition Energies [eV]:
[[0.01733871]
[0.01716977]]

The addition energy can be viewed as the vertical distance between the two red dashed lines in Fig. 5.4.8.

Alternatively, within the constant interaction model [HKP+07], the addition energy can be expressed as:

\[E_{\mathrm{add}}(N) = E_C + \Delta E(N),\]

where \(E_C\) is the charging energy arising from electron-electron interactions, and

\[\Delta E(N) = \epsilon_{N+1} - \epsilon_{N}\]

is the single-particle energy spacing, with \(\epsilon_N\) denoting the energy of the \(N\)-th single-particle state. Here, each single-particle energy level is shown separately for each degenerate state (spin and valley), so that electrons occupying the same orbital but different degenerate slots are distinguished.

So far, we have considered the chemical potentials corresponding to transitions between states with \(N-1=12\), \(N=13\), and \(N+1=14\) electrons in the SET.

In silicon, where the degeneracy \(g=4\), the state with \(N-1=12\) electrons fills the first three single-particle states completely. The next electron (\(N=13\)) is the first to occupy the fourth single-particle state, and the following electron (\(N=14\)) fills the second slot of that same state. Because these two electrons occupy the same single-particle state, the single-particle energy spacing

\[\Delta E(13) = \epsilon_{14} - \epsilon_{13} = 0.\]

Thus, the addition energy for \(N=13\) is determined solely by the charging energy \(E_C\).

 # Charging energy
E_C = addition_energies

Note

In general, the charging energy \(E_C\) depends on the single-particle energy level spacing \(\Delta E\). This quantity can be computed from the single-particle energy spectrum obtained when solving the Schrödinger equation and saved in the path_local_out / "energies.txt" file (see Computing the chemical potential of an SET device).

5.4.3. Gate Capacitance

The final quantity we can extract from the Coulomb peak simulations is the gate capacitance \(C_{\mathrm{SET}}\) between the SET gate, "plunger_gate_1_bnd" and the SET island. This capacitance determines how effectively the gate voltage shifts the energy levels of the SET [HKP+07]. From Eq. (2) of [HKP+07], we have

\[\mu(N) = \alpha V_{\mathrm{SET}} + b = -\frac{E_C}{e} C_{\mathrm{SET}} V_{\mathrm{SET}} + b,\]

where \(\alpha\) and \(b\) are the lever arm and intercept computed earlier. Rearranging the above equation, we can express the gate capacitance as

\[C_{\mathrm{SET}} = -\frac{\alpha e}{E_C}.\]

Using the lever arm and charging energy computed earlier, we can compute the gate capacitance as follows:

# Capacitance between SET gate and the SET ----------------------------------
C_gate_set = -lever_arm / (E_C / ct.e) * ct.e * 1e18  # in attoFarads [aF]
print("Capacitance between SET gate and the SET [aF]:")
print(C_gate_set)

which prints

Capacitance between SET gate and the SET [aF]:
[5.10759045 5.07636812]

5.5. Full code

__copyright__ = "Copyright 2022-2025, Nanoacademic Technologies Inc."

import pathlib
import numpy as np
import matplotlib.pyplot as plt
from qtcad.device import io
from qtcad.device import constants as ct
from qtcad.device.mesh3d import SubMesh
from qtcad.device import SubDevice
from double_dot_fdsoi import get_double_dot_fdsoi, V_set
from chemical_potential import compute_chemical_potential

# Files -----------------------------------------------------------------------
script_dir           = pathlib.Path(__file__).parent.resolve()

path_out             = script_dir / "output" 
path_local_out       = path_out / "coulomb_peaks"
phi_in_file          = path_out / "phi" / f"phi.hdf5"

leverarm_file        = str(path_local_out / "lever_arm.txt")
lever_arm_plot_file  = str(path_local_out / "lever_arm_plot.png")


# Setup ----------------------------------------------------------------------

# Create the device
d, set_region, qd_region = get_double_dot_fdsoi()
# load the potential 
d.set_potential(io.load(phi_in_file, var_name="phi"))
# Create sub-device for the SET region
set_mesh = SubMesh(d.mesh, set_region)
set_device = SubDevice(d, set_mesh)
set_device.set_insulator_boundaries()

# V_set values
V_SET = np.array([V_set, V_set - 0.05])
# Number of particles
Nmin = 12
Nmax = 14
# Generate data
Mu = compute_chemical_potential(V_SET, Nmin, Nmax, d, set_device, out_path=path_local_out)

# Coulomb peak positions ------------------------------------------------------

# Fit a line: y = alpha*x + b
lever_arm = np.zeros(Mu.shape[1])
intercept = np.zeros(Mu.shape[1])
for n in range(Mu.shape[1]):
    alpha, b = np.polyfit(V_SET, Mu[:, n] / ct.e, 1)
    lever_arm[n] = alpha
    intercept[n] = b
    print(f"n={n + Nmin + 1}: Lever Arm [eV/V] = {alpha}, Intercept [eV] = {b}")

# Compute Coulomb peak positions
# (x-intercept)
V_peaks = -intercept / lever_arm
print("Coulomb Peak Positions [V]:")
print(V_peaks)
# (Distance between peaks)
Delta_V_peaks = np.diff(V_peaks)
print("Distance between Coulomb Peaks [V]:")
print(Delta_V_peaks)

# Save
header = "Lever Arm [eV/V]    Intercept [eV]    Coulomb Peak Position [V]"
np.savetxt(leverarm_file, np.column_stack((lever_arm, intercept, V_peaks)), fmt="%.10f", header=header)

# Plot
plt.figure(figsize=(8, 5))
# (Add data to plot)
DeltaV = np.max(np.insert(V_SET, 0, np.max(V_peaks))) - np.min(np.insert(V_SET, 0, np.min(V_peaks)))
V_line = np.linspace(min(V_SET) - DeltaV, max(V_SET) + DeltaV, 200)
plt.axhline(y=0, color="black", linewidth=3)
for j, alpha in enumerate(lever_arm):
    b = intercept[j] # Intercept for given value of N
    mu_line = alpha * V_line + b
    plt.plot(V_line, mu_line, color="red", linewidth=2.0, linestyle="--") 
    plt.axvline(V_peaks[j], color="black", linestyle="--", linewidth=2.)
for j in range(Mu.shape[1]):
    plt.scatter(V_SET, Mu[:, j] / ct.e, marker="o", s=60, label=f"N = {j + Nmin + 1}")     

# (Labels, legend, grid)
plt.xlabel("$V_{SET}$ [V]", fontsize=16)
plt.ylabel("$\\mu$ [eV]", fontsize=16)
plt.legend(fontsize=14)
plt.grid(True)
# (Save)
plt.tight_layout()
plt.savefig(lever_arm_plot_file, dpi=300)

# Compute addition energies ---------------------------------------------------
addition_energies = Mu[:, 1:] - Mu[:, :-1]
print("Addition Energies [eV]:")
print(addition_energies / ct.e)

 # Charging energy
E_C = addition_energies

# Capacitance between SET gate and the SET ----------------------------------
C_gate_set = -lever_arm / (E_C / ct.e) * ct.e * 1e18  # in attoFarads [aF]
print("Capacitance between SET gate and the SET [aF]:")
print(C_gate_set)