6. Detecting charges with an SET

6.1. Requirements

6.1.1. Software components

  • QTCAD

6.1.2. Python script

  • qtcad/examples/practical_application/FDSOI/4-charge_detection.py

6.1.3. References

6.2. Briefing

In Simulating the Coulomb peak positions of an SET, we demonstrated how to simulate the Coulomb peaks of a single-electron transistor (SET) formed in a fully depleted silicon-on-insulator (FD-SOI) device. In this part of the practical application, we build on those results to show how the positions of these peaks shift in the presence of a nearby charge, such as an electron confined in a neighboring quantum dot. Such peak shifts provide a mechanism for detecting that charge, which is a key operation in certain spin-qubit readout schemes [HKP+07].

To simulate this effect, we repeat the Coulomb-peak calculations from Simulating the Coulomb peak positions of an SET. However, before computing the chemical potentials, we now introduce an additional electron into the neighboring (qubit) quantum dot. We do this by first solving the Schrödinger equation in that dot and then using the ground-state wavefunction to define a fixed charge density in the subsequent Schrödinger–Poisson calculation for the SET dot. This additional fixed charge density modifies the electrostatic potential experienced by the electrons in the SET, thereby shifting the positions of the Coulomb peaks. Experimentally, these peak shifts can be monitored to detect the presence or absence of the electron in the qubit dot.

6.3. Setup

6.3.2. File paths

Next, we define the paths to the input files and to the output files where we will save the results of this simulation.

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

path_out             = script_dir / "output"
path_local_out       = path_out / "charge_sensing"
phi_in_file          = path_out / "coulomb_peaks" /  f"phi_vset{V_set:.2f}_N13.hdf5"
la_in_file           = path_out / "coulomb_peaks" / "lever_arm.txt"

# Visualization
wf_qubit_slice_file  = path_local_out / "wf_qubit_slice.png"
rho_qubit_slice_file = path_local_out / "rho_qubit_slice.png"
lever_arm_plot_file  = path_local_out / "lever_arm_plot.svg"

6.4. Setting the charge in the qubit quantum dot

Now, we create the device d using the get_double_dot_fdsoi function defined in double_dot_fdsoi.py and load the electrostatic potential previously computed in Simulating the Coulomb peak positions of an SET. After this potential is set, we also use the set_V_from_phi method to update the electron confinement potential, which is stored in the Device attribute V (i.e. d.V). The confinement potential is computed from the electrostatic potential phi and material properties (see Anderson’s rule for more details).

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

# Create the device
d, set_region, qd_region = get_double_dot_fdsoi(mesh_file_name="refined_dqdfdsoi.msh")
# load the potential
d.set_potential(io.load(phi_in_file, var_name="phi"))
d.set_V_from_phi()

Next, we solve the Schrödinger equation in the qubit quantum dot. We first create a SubDevice associated to the quantum-dot region and then solve the Schrödinger equation using the Schrödinger Solver.

# Qubit quantum dot -----------------------------------------------------------
# QD subdevice
qd_mesh = SubMesh(d.mesh, qd_region)
qd_device = SubDevice(d, qd_mesh)
qd_device.set_insulator_boundaries()
# Solve the Schrödinger equation
schrod_solver = SchrodingerSolver(qd_device, solver_params=schrod_params)
schrod_solver.solve()
qd_device.print_energies()

Note

The potential we loaded earlier from phi_in_file corresponds to the solution of the Poisson equation for a specific device configuration, i.e. one with a gate voltage V_set applied to the SET gate and \(N=13\) electrons confined to the SET. In practice, the exact wavefunctions of the electrons confined to the qubit quantum dot will depend on both the SET-gate bias and the electron occupancy of the SET. For simplicity, we do not model these dependencies explicitly and instead assume that the potential loaded from phi_in_file is sufficiently accurate for all relevant device configurations considered here. This approximation is reasonable provided that the qubit quantum-dot potential is mostly determined by the gate layout and applied DC biases, and the fine details of the SET electron distribution do not change the qubit-dot confinement appreciably. The validity of the approximation can be verified by repeating the calculation for different SET gate voltages and occupancies and comparing the results.

After obtaining the solution, we extract the ground-state wavefunction, compute the corresponding charge density, and set it as a fixed volume charge density over the full device, d.

# Charge density from QD ground state
rho0 = np.zeros_like(d.phi)
rho0[qd_device.mesh.nodes_in_parent] = -ct.e * np.abs(qd_device.eigenfunctions[:, 0])**2
d.set_vol_charge_density(rho0)

We begin by creating an array, rho0, defined over all nodes of the full device, d. For the nodes belonging to the qubit quantum-dot SubDevice, we assign to rho0 the charge density associated with the ground-state wavefunction, \(\psi_0(\mathbf{r})\), obtained from the Schrödinger Solver. The indices of these nodes within the full device are retrieved using the nodes_in_parent attribute of the relevant SubMesh. The SubMesh corresponding to the qubit quantum dot is accessible through the mesh attribute of the qubit-dot device, qd_device.

Once the charge density is assembled, we assign it as a fixed volume charge density in the full device d using the set_vol_charge_density method of the Device class. Physically, the Coulomb interaction between the electrons confined in the SET and the fixed charge in the nearby qubit quantum dot increases the chemical potential of the SET. This manifests experimentally as a shift of the Coulomb-blockade peaks toward higher gate voltages. In other words, the additional electron in the qubit quantum dot introduces extra electrostatic repulsion, requiring a higher gate voltage on the SET gate to realign the SET chemical potential with the Fermi level.

We can also visualize the ground-state wavefunction and the associated charge density using the save_slice function defined in double_dot_fdsoi.py.

# Visualize QD wavefunction and charge density
# (Wavefunction)
save_slice(qd_device, qd_device.eigenfunctions[:, 0], wf_qubit_slice_file,
    label="Ground-state, $\\psi_0$ [1 / m$^{3/2}$]")
# (Charge density)
save_slice(d, d.rho_0 / ct.e, rho_qubit_slice_file,
    label="$\\rho_0$ [$e$ / m$^3$]")
Slice of the ground-state wavefunction in the qubit quantum dot.

Fig. 6.4.10 Slice of the ground-state wavefunction in the qubit quantum dot, \(\psi_0(\mathbf{r})\), obtained by solving the Schrödinger equation using the potential loaded from phi_in_file. The wavefunction is plotted over the qubit quantum-dot region.

Slice of the charge density associated with the ground-state wavefunction.

Fig. 6.4.11 Slice of the charge density present in the entire Device d. The only charge density present is that associated with the ground-state wavefunction \(\rho_0(\mathbf{r}) = -e|\psi_0(\mathbf{r})|^2\). The charge density is plotted over the entire Device, showing that it is localized within the qubit quantum-dot region.

Note

The set_vol_charge_density method sets the background charge density of the Device, stored in the rho_0 attribute. This charge density is added to the full charge density rho by the Schrödinger–Poisson Solver via the solve method (see Poisson’s equation and Fig. 6.5.4, below).

6.5. Computing chemical potential of the FD-SOI SET

Then, we proceed as before (see Compute chemical potential of the FD-SOI SET) by creating a sub-device for the SET region and computing the chemical potential as a function of the gate voltage.

# Compute the chemical potential ----------------------------------------------
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 = 13
# Generate data
Mu = compute_chemical_potential(V_SET, Nmin, Nmax, d, set_device, out_path=path_local_out)

Here we set Nmin = 12 and Nmax = 13 to compute the chemical potential corresponding to the transitions between the \(N-1=12\) and the \(N=13\) ground states. We are computing a single transition here to minimize computational cost for this tutorial. The user is free to adjust these values to compute the chemical potentials for additional transitions as needed.

In addition to computing the chemical potentials for the different charge and voltage configurations, the compute_chemical_potential function will also save the electrostatic potential and charge density for each case in the output directory, path_local_out. Here we show the total charge density for the \(N=13\) case at \(V_{\mathrm{SET}} = 1.25 \, \mathrm{V}\).

Slice of the total charge density in the SET for N=13 electrons.

Fig. 6.5.4 Slice of the total charge density in the SET region for \(N=13\) electrons confined in the SET at gate voltage \(V_{\mathrm{SET}} = 1.25 \, \mathrm{V}\). The charge density is plotted over the entire Device, showing that there is a large accumulation of charge within the SET region and a smaller amount of charge localized in the qubit quantum dot due to the fixed charge density defined there earlier.

We see from Fig. Fig. 6.5.4 that the total charge density contains contributions both from the electrons confined to the SET and the qubit quantum dot. While the SET contains a large accumulation of charge due to the \(N=13\) electrons confined therein, there is also a smaller amount of charge localized in the qubit quantum dot due to the fixed charge density defined there earlier.

6.6. Data analysis

In the previous section we computed the chemical potential of the SET for the \(N-1=12\) to \(N=13\) transition at two different gate voltages, V_SET = [V_set, V_set - 0.05]. Using this data, along with the data produced in Simulating the Coulomb peak positions of an SET, we can extract a quantitative measure of the SET’s sensitivity to charges present in the qubit quantum dot.

6.6.1. Coulomb peak positions

As before (see Coulomb peak positions), we start by computing the lever arm and intercept of the linear fits to the chemical-potential data. From there, we can then extract the position of the Coulomb peak for the relevant transition. Again, this is the voltage at which the chemical potential aligns with the Fermi level, \(E_F = 0\).

# Coulomb peak position -------------------------------------------------------

# Fit a line: y = alpha*x + b
alpha, b = np.polyfit(V_SET, Mu / ct.e, 1)
print("--------------------------------------------------")
print(f"Charge Sensing: Lever Arm [eV/V] = {alpha}, Intercept [eV] = {b}")
peak_pos = -b / alpha
print(f"Charge Sensing: Coulomb Peak Position [V]: {peak_pos}")

This code snippet produces the following output:

--------------------------------------------------
Charge Sensing: Lever Arm [eV/V] = [-0.55586722], Intercept [eV] = [0.6913577]
Charge Sensing: Coulomb Peak Position [V]: [1.24374612]

Next, we load the chemical potential data computed in Simulating the Coulomb peak positions of an SET for the same transition (\(N-1=12\) to \(N=13\) electrons in the SET).

# Load lever arm data from previous Coulomb peak calculation
charge_neutral_la   = np.loadtxt(la_in_file)
lever_arm           = charge_neutral_la[:, 0]
intercept           = charge_neutral_la[:, 1]
V_peaks             = charge_neutral_la[:, 2]

From there, we can compare the Coulomb peak position computed here with the position computed previously without the extra charge in the qubit quantum dot.

# Compare peak positions
d_pos = np.abs(peak_pos[0] - V_peaks[0])
print(f"Shift in peak position [V]: {d_pos}")

We find that the presence of the extra electron in the qubit quantum dot shifts the position of the Coulomb peak by a certain amount.

Shift in peak position [V]: 0.00031764797895550423

6.6.2. Plotting

We can also visualze the chemical potential by plotting the linear fits for both cases: with and without the extra charge in the qubit quantum dot (\(Q_{\mathrm{qubit}}=1\) and \(Q_{\mathrm{qubit}}=0\), respectively).

# Plot
plt.figure(figsize=(8, 5))
# (Add data to plot)
V_line = np.linspace(np.min([V_peaks[0], peak_pos[0]]) - d_pos, np.max([V_peaks[0], peak_pos[0]]) + d_pos, 200)
plt.axhline(y=0, color="black", linewidth=3) # E_F = 0
plt.plot(
    V_line, alpha * V_line + b, color="blue", linewidth=2.0, linestyle="--", label="$Q_{\\mathrm{qubit}}=1$"
    ) # Plot data for charged qubit dot

alpha_q1 = lever_arm[0]
b_q1 = intercept[0]
mu_line = alpha_q1 * V_line + b_q1
plt.plot(
    V_line, mu_line, color="red", linewidth=2.0, linestyle="--", label="$Q_{\\mathrm{qubit}}=0$"
    ) # Plot data for empty qubit dot

# (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:

Chemical potential vs. gate voltage.

Fig. 6.6.3 Chemical potential of the SET \(\mu(N=13)\) for the \(N-1=12 \rightarrow N=13\) transition as a function of the SET gate voltage, \(V_{\mathrm{SET}}\). The red line corresponds to the case where the qubit quantum dot is empty (\(Q_{\mathrm{qubit}}=0\)), while the blue line corresponds to the case where the qubit quantum dot contains one electron (\(Q_{\mathrm{qubit}}=1\)). The horizontal black line indicates the Fermi level, \(E_F = 0\). The shift in the Coulomb peak position due to the extra charge in the qubit dot is given by the horizontal distance between the two dashed lines at \(\mu(N=13) = 0\). By monitoring the position of the SET Coulomb peaks, the charge state of the qubit dot can be measured.

6.7. 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 qtcad.device.schrodinger import Solver as SchrodingerSolver
from double_dot_fdsoi import get_double_dot_fdsoi, V_set, save_slice
from chemical_potential import compute_chemical_potential, schrod_params

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

path_out             = script_dir / "output"
path_local_out       = path_out / "charge_sensing"
phi_in_file          = path_out / "coulomb_peaks" /  f"phi_vset{V_set:.2f}_N13.hdf5"
la_in_file           = path_out / "coulomb_peaks" / "lever_arm.txt"

# Visualization
wf_qubit_slice_file  = path_local_out / "wf_qubit_slice.png"
rho_qubit_slice_file = path_local_out / "rho_qubit_slice.png"
lever_arm_plot_file  = path_local_out / "lever_arm_plot.svg"

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

# Create the device
d, set_region, qd_region = get_double_dot_fdsoi(mesh_file_name="refined_dqdfdsoi.msh")
# load the potential
d.set_potential(io.load(phi_in_file, var_name="phi"))
d.set_V_from_phi()

# Qubit quantum dot -----------------------------------------------------------
# QD subdevice
qd_mesh = SubMesh(d.mesh, qd_region)
qd_device = SubDevice(d, qd_mesh)
qd_device.set_insulator_boundaries()
# Solve the Schrödinger equation
schrod_solver = SchrodingerSolver(qd_device, solver_params=schrod_params)
schrod_solver.solve()
qd_device.print_energies()

# Charge density from QD ground state 
rho0 = np.zeros_like(d.phi)
rho0[qd_device.mesh.nodes_in_parent] = -ct.e * np.abs(qd_device.eigenfunctions[:, 0])**2
d.set_vol_charge_density(rho0)

# Visualize QD wavefunction and charge density 
# (Wavefunction)
save_slice(qd_device, qd_device.eigenfunctions[:, 0], wf_qubit_slice_file, 
    label="Ground-state, $\\psi_0$ [1 / m$^{3/2}$]")
# (Charge density)
save_slice(d, d.rho_0 / ct.e, rho_qubit_slice_file, 
    label="$\\rho_0$ [$e$ / m$^3$]")

# Compute the chemical potential ----------------------------------------------
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 = 13
# Generate data
Mu = compute_chemical_potential(V_SET, Nmin, Nmax, d, set_device, out_path=path_local_out)

# Coulomb peak position -------------------------------------------------------

# Fit a line: y = alpha*x + b
alpha, b = np.polyfit(V_SET, Mu / ct.e, 1)
print("--------------------------------------------------")
print(f"Charge Sensing: Lever Arm [eV/V] = {alpha}, Intercept [eV] = {b}")
peak_pos = -b / alpha
print(f"Charge Sensing: Coulomb Peak Position [V]: {peak_pos}")

# Load lever arm data from previous Coulomb peak calculation
charge_neutral_la   = np.loadtxt(la_in_file)
lever_arm           = charge_neutral_la[:, 0]
intercept           = charge_neutral_la[:, 1]
V_peaks             = charge_neutral_la[:, 2]

# Compare peak positions
d_pos = np.abs(peak_pos[0] - V_peaks[0])
print(f"Shift in peak position [V]: {d_pos}")

# Plot
plt.figure(figsize=(8, 5))
# (Add data to plot)
V_line = np.linspace(np.min([V_peaks[0], peak_pos[0]]) - d_pos, np.max([V_peaks[0], peak_pos[0]]) + d_pos, 200)
plt.axhline(y=0, color="black", linewidth=3) # E_F = 0
plt.plot(
    V_line, alpha * V_line + b, color="blue", linewidth=2.0, linestyle="--", label="$Q_{\\mathrm{qubit}}=1$"
    ) # Plot data for charged qubit dot

alpha_q1 = lever_arm[0] 
b_q1 = intercept[0] 
mu_line = alpha_q1 * V_line + b_q1
plt.plot(
    V_line, mu_line, color="red", linewidth=2.0, linestyle="--", label="$Q_{\\mathrm{qubit}}=0$"
    ) # Plot data for empty qubit dot

# (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)