2. Creating the Tunnel Falls Device

2.1. Requirements

2.1.1. Software components

  • QTCAD

2.1.2. Python script

  • qtcad/examples/practical_application/Tunnel_Falls_DAPS/double_dot_tunnel_falls.py

2.1.3. References

2.2. Briefing

This is the second part of the Tunnel Falls practical application. The module presented in this section is an importable support module rather than a script that is meant to be run directly. Its main entry point, get_double_dot_tunnel_falls, loads the mesh generated in Mesh generation using QTCAD Builder and creates the corresponding QTCAD Device object for the later Poisson and Schrödinger calculations.

The device model follows the simplified center-region geometry introduced in the previous section. It has five finger gates B4, P5, B5, P6, and B6, a buried screening gate SG, and the center-screening gate CS. The B finger gates are barrier gates, while the P finger gates are plunger gates. The material stack is based on the Tunnel Falls \(\mathrm{Si}/\mathrm{SiGe}\) heterostructure of [MEB+25], which uses relaxed \(\mathrm{Si}_{0.7}\mathrm{Ge}_{0.3}\) barriers and a \(4.6\,\mathrm{nm}\) strained \(\mathrm{Si}_{0.972}\mathrm{Ge}_{0.028}\) quantum well. In the model, the quantum well is approximated as pure silicon. This simplification avoids introducing a separate material parameter set for strained \(\mathrm{SiGe}\).

2.3. Setup

2.3.2. Default values

We then define some default values used by the different scripts in this practical application.

# Default values --------------------------------------------------------------
nm = 1e-9
Ge_fraction = 0.30  # Ge fraction in the SiGe layers.

The Mesh is generated in nanometers, so nm is used as the length-scaling factor when loading the Mesh. The value Ge_fraction = 0.30 is taken from the relaxed \(\mathrm{Si}_{0.7}\mathrm{Ge}_{0.3}\) barriers shown in the Tunnel Falls heterostructure cross-section of [MEB+25].

Next, we define the parameters used to determine the location of the slices we will use for visualization.

# Visualization parameters
buffer_thickness = 50.0
qw_thickness = 4.6
normal = (0.0, 0.0, 1.0)
origin = (0.0, 0.0, (buffer_thickness + 0.5 * qw_thickness) * nm)

The slices will be taken normal to the growth direction in the center of the quantum well.

Finally, we define the default gate voltages used by the device helper.

# Default gate voltages
V_B4 = 1.0
V_P5 = 2.5
V_B5 = 1.2
V_P6 = 2.5
V_B6 = 1.0
V_SG = -0.5
V_CS = -0.5

The labels follow the Tunnel Falls layout: B gates control tunnel barriers, P gates act as plungers, SG is the buried screening gate, and CS is the center-screening gate separating the qubit side from the opposite side of the device. The plunger defaults set P5 and P6 equal to one another; the detuning script later raises one plunger and lowers the other around this common value.

These voltages should be understood as an illustrative bias configuration for the simplified electrostatic model, rather than as calibrated experimental gate settings. The positive plunger voltages lower the electron potential energy under P5 and P6 enough to form the two dot minima. The barrier-gate voltages are kept lower than the plunger voltages to maintain confinement between and around the dots, with B5 set slightly higher than B4 and B6 to avoid over-isolating the two dots from one another. The negative SG and CS voltages help screen and deplete regions away from the active double-dot channel. This combination provides a stable, nearly symmetric starting configuration for the examples, while leaving the plunger voltages available for the detuning sweep.

2.4. Visualization helpers

The first helper saves a two-dimensional slice of a quantity defined over the full Device or a SubDevice.

def save_slice(
    device: Device | SubDevice,
    data: np.ndarray,
    file_name: str | Path,
    label: str,
) -> None:
    """Helper function to save a slice of the device data.

    Args:
        device: The Device object containing the data.
        data: The data array to be sliced and saved.
        file_name: The filename to save the slice plot.
        label: The label for the colorbar axis.
    """
    plot_slice(
        device.mesh,
        data,
        normal=normal,
        origin=origin,
        cb_axis_label=label,
        path=file_name,
        show_figure=False,
    )

The second helper converts an eigenfunction of a Device or a SubDevice into a scalar probability density.

def state_probability_density(device: Device | SubDevice, state: int) -> np.ndarray:
    """Return the scalar probability density associated with an eigenstate.

    Args:
        device: Device containing previously computed eigenfunctions.
        state: Single-particle state index for which to compute the density.

    Returns:
        Scalar probability density evaluated on the device mesh.
    """
    density = np.abs(device.eigenfunctions[:, state]) ** 2
    if density.ndim > 1:
        density = np.sum(density, axis=-1)
    return density

For a single-band calculation, this is simply \(|\psi_i(\mathbf{r})|^2\), where \(\psi_i(\mathbf{r})\) is the \(i^{\mathrm{th}}\) energy eigenfunction. For a multicomponent calculation, such as the explicit two-valley model used later in this practical application, the final axis indexes valley components. Summing over that final axis gives the total position-dependent density of the selected eigenstate.

2.5. Creating the device

Finally, we define the function that creates and returns the Tunnel Falls Device model.

def get_double_dot_tunnel_falls(
    mesh_file_name: str = "tunnel_falls_double_dot.msh",
    V_B4: float = V_B4,
    V_P5: float = V_P5,
    V_B5: float = V_B5,
    V_P6: float = V_P6,
    V_B6: float = V_B6,
    V_SG: float = V_SG,
    V_CS: float = V_CS,
) -> tuple[Device, list[str], list[str]]:
    """Create a Device object for the simplified Tunnel Falls double-dot mesh.

    Args:
        mesh_file_name: Name of the mesh file stored in the local ``meshes``
            directory.
        V_B4: Potential applied to gate ``B4``.
        V_P5: Potential applied to gate ``P5``.
        V_B5: Potential applied to gate ``B5``.
        V_P6: Potential applied to gate ``P6``.
        V_B6: Potential applied to gate ``B6``.
        V_SG: Potential applied to gate ``SG``.
        V_CS: Potential applied to gate ``CS``.

    Returns:
        Device, left-dot region labels, and right-dot region labels.
    """

The first lines inside the function load the Mesh and instantiate the Device object.

script_dir = Path(__file__).parent.resolve()
# Load the mesh.
mesh_file = script_dir / "meshes" / mesh_file_name
mesh = Mesh(nm, mesh_file)

# Create the device.
device = Device(mesh, conf_carriers="e")

The nm scaling factor converts the Mesh coordinates from nanometers to SI units. The default mesh file is generated in Mesh generation using QTCAD Builder. Later scripts may instead pass different mesh files defined on the same geometry. For example, a mesh that has been refined through an adaptive Poisson Solver can be used. The conf_carriers="e" argument tells QTCAD that the confined carriers are electrons.

Next, we define the two semiconductor material models used in the \(\mathrm{Si}/\mathrm{SiGe}\) heterostructure.

# Generate materials for quantum-well and barrier regions.
sige_barrier = mt.SiGe_DFT.copy()
sige_barrier.set_alloy_composition(Ge_fraction)

strained_well = mt.Si_strained_on_SiGe.copy()
strained_well.set_alloy_composition(Ge_fraction)

The relaxed SiGe barriers use the SiGe_DFT material with the \(\mathrm{Ge}\) fraction set to 0.30, following [MEB+25]. The quantum well uses Si_strained_on_SiGe with the same alloy composition because the \(\mathrm{Si}_{0.972}\mathrm{Ge}_{0.028}\) well is strained by the surrounding relaxed \(\mathrm{Si}/\mathrm{SiGe}\). Here, for simplicity, we approximate \(\mathrm{Si}_{0.972}\mathrm{Ge}_{0.028}\), which is \(97.2\%\) \(\mathrm{Si}\), by pure \(\mathrm{Si}\). The bandgap and band alignment parameters stored in the SiGe_DFT and Si_strained_on_SiGe Materials were obtained via atomistic first-principles calculations based on the RESCU package [VPP+26, NanoacademicTInc26].

We then assign materials to the main physical regions created by Builder.

# Create the regions.
device.new_region("relaxed_buffer", sige_barrier)
device.new_region("quantum_well", strained_well)
device.new_region("upper_barrier", sige_barrier)
device.new_region("si_cap", mt.Si)
device.new_region("gate_oxide_sio2", mt.SiO2)
device.new_region("gate_oxide_hfo2", mt.HfO2)
device.new_region("screening_ild", mt.SiO2)

Each call to new_region maps a physical volume name from the Mesh to a Material. These names must match the physical volume names written by the Builder script.

The dot subregions are assigned the same materials as the layers they overlap.

device.new_region("relaxed_buffer.dot_left", sige_barrier)
device.new_region("quantum_well.dot_left", strained_well)
device.new_region("upper_barrier.dot_left", sige_barrier)
device.new_region("relaxed_buffer.dot_right", sige_barrier)
device.new_region("quantum_well.dot_right", strained_well)
device.new_region("upper_barrier.dot_right", sige_barrier)

Next, we define the electrostatic gate boundaries. The gate work function is chosen to correspond to the midgap energy of the strained silicon well. In practice, the gate work function depends on details of the fabrication process and can differ from the value chosen here.

# Set up boundary conditions
gate_work_function = strained_well.Eg / 2 + strained_well.chi
gate_biases = {
    "B4": V_B4,
    "P5": V_P5,
    "B5": V_B5,
    "P6": V_P6,
    "B6": V_B6,
    "SG": V_SG,
    "CS": V_CS,
}
for gate_name, gate_bias in gate_biases.items():
    device.new_gate_bnd(gate_name, gate_bias, gate_work_function)

The left and right dot region lists are then assembled.

# Create the double quantum dot region.
left_dot_region = [
    "relaxed_buffer.dot_left",
    "quantum_well.dot_left",
    "upper_barrier.dot_left",
]
right_dot_region = [
    "relaxed_buffer.dot_right",
    "quantum_well.dot_right",
    "upper_barrier.dot_right",
]

These lists are returned together with the full Device. The linear Poisson Solver uses the full Device, while the Schrödinger workflow uses these lists to make a SubDevice restricted to the left and right dot regions.

Finally, the function returns the constructed objects.

return device, left_dot_region, right_dot_region

The valley physics is intentionally not set in this helper. In the final detuning-spectrum script, it is introduced explicitly through a position-dependent two-valley \(\mathbf{k} \cdot \mathbf{p}\) model which will be generated in the next section of this practical application.

2.6. Full code

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

from pathlib import Path
import numpy as np
from qtcad.device import Device, SubDevice
from qtcad.device import materials as mt
from qtcad.device.analysis import plot_slice
from qtcad.device.mesh3d import Mesh

# Default values --------------------------------------------------------------
nm = 1e-9
Ge_fraction = 0.30  # Ge fraction in the SiGe layers.

# Visualization parameters
buffer_thickness = 50.0
qw_thickness = 4.6
normal = (0.0, 0.0, 1.0)
origin = (0.0, 0.0, (buffer_thickness + 0.5 * qw_thickness) * nm)

# Default gate voltages
V_B4 = 1.0
V_P5 = 2.5
V_B5 = 1.2
V_P6 = 2.5
V_B6 = 1.0
V_SG = -0.5
V_CS = -0.5

def save_slice(
    device: Device | SubDevice,
    data: np.ndarray,
    file_name: str | Path,
    label: str,
) -> None:
    """Helper function to save a slice of the device data.

    Args:
        device: The Device object containing the data.
        data: The data array to be sliced and saved.
        file_name: The filename to save the slice plot.
        label: The label for the colorbar axis.
    """
    plot_slice(
        device.mesh,
        data,
        normal=normal,
        origin=origin,
        cb_axis_label=label,
        path=file_name,
        show_figure=False,
    )

def state_probability_density(device: Device | SubDevice, state: int) -> np.ndarray:
    """Return the scalar probability density associated with an eigenstate.

    Args:
        device: Device containing previously computed eigenfunctions.
        state: Single-particle state index for which to compute the density.

    Returns:
        Scalar probability density evaluated on the device mesh.
    """
    density = np.abs(device.eigenfunctions[:, state]) ** 2
    if density.ndim > 1:
        density = np.sum(density, axis=-1)
    return density

def get_double_dot_tunnel_falls(
    mesh_file_name: str = "tunnel_falls_double_dot.msh",
    V_B4: float = V_B4,
    V_P5: float = V_P5,
    V_B5: float = V_B5,
    V_P6: float = V_P6,
    V_B6: float = V_B6,
    V_SG: float = V_SG,
    V_CS: float = V_CS,
) -> tuple[Device, list[str], list[str]]:
    """Create a Device object for the simplified Tunnel Falls double-dot mesh.

    Args:
        mesh_file_name: Name of the mesh file stored in the local ``meshes``
            directory.
        V_B4: Potential applied to gate ``B4``.
        V_P5: Potential applied to gate ``P5``.
        V_B5: Potential applied to gate ``B5``.
        V_P6: Potential applied to gate ``P6``.
        V_B6: Potential applied to gate ``B6``.
        V_SG: Potential applied to gate ``SG``.
        V_CS: Potential applied to gate ``CS``.

    Returns:
        Device, left-dot region labels, and right-dot region labels.
    """

    script_dir = Path(__file__).parent.resolve()
    # Load the mesh.
    mesh_file = script_dir / "meshes" / mesh_file_name
    mesh = Mesh(nm, mesh_file)

    # Create the device.
    device = Device(mesh, conf_carriers="e")

    # Generate materials for quantum-well and barrier regions.
    sige_barrier = mt.SiGe_DFT.copy()
    sige_barrier.set_alloy_composition(Ge_fraction)

    strained_well = mt.Si_strained_on_SiGe.copy()
    strained_well.set_alloy_composition(Ge_fraction)

    # Create the regions.
    device.new_region("relaxed_buffer", sige_barrier)
    device.new_region("quantum_well", strained_well)
    device.new_region("upper_barrier", sige_barrier)
    device.new_region("si_cap", mt.Si)
    device.new_region("gate_oxide_sio2", mt.SiO2)
    device.new_region("gate_oxide_hfo2", mt.HfO2)
    device.new_region("screening_ild", mt.SiO2)

    device.new_region("relaxed_buffer.dot_left", sige_barrier)
    device.new_region("quantum_well.dot_left", strained_well)
    device.new_region("upper_barrier.dot_left", sige_barrier)
    device.new_region("relaxed_buffer.dot_right", sige_barrier)
    device.new_region("quantum_well.dot_right", strained_well)
    device.new_region("upper_barrier.dot_right", sige_barrier)

    # Set up boundary conditions
    gate_work_function = strained_well.Eg / 2 + strained_well.chi
    gate_biases = {
        "B4": V_B4,
        "P5": V_P5,
        "B5": V_B5,
        "P6": V_P6,
        "B6": V_B6,
        "SG": V_SG,
        "CS": V_CS,
    }
    for gate_name, gate_bias in gate_biases.items():
        device.new_gate_bnd(gate_name, gate_bias, gate_work_function)

    # Create the double quantum dot region.
    left_dot_region = [
        "relaxed_buffer.dot_left",
        "quantum_well.dot_left",
        "upper_barrier.dot_left",
    ]
    right_dot_region = [
        "relaxed_buffer.dot_right",
        "quantum_well.dot_right",
        "upper_barrier.dot_right",
    ]

    return device, left_dot_region, right_dot_region