Technical Questions

1. How to customize material properties in QTCAD?

You can customize materials in QTCAD by creating instances of the Material class with your desired properties or by using built-in materials. Key steps include:

  • Create a custom material: Instantiate Material with a dictionary of properties like permittivity (eps), bandgap (Eg), electron affinity (chi), and effective masses (mc, mv).

    from qtcad.device import materials as mt

custom_mat = mt.Material({
    'name': 'CustomMaterial',
    'eps': 12.0,
    'Eg': 1.1,
    'chi': 4.0,
    'mc': [0.26],
    'mv': [0.45]
})

  • Use predefined materials: QTCAD provides predefined materials like silicon (mt.Si) and silicon oxide (mt.SiO2) with default properties which can be modified.

  • Define alloy materials: Use fitting functions to model alloy compositions:

  • Assign materials to device regions: Use the device’s Device.new_region method to specify materials and doping for each region.

    device.new_region('barrier', mt.SiO2)
device.new_region('confined', mt.Si, pdoping=0, ndoping=0)

  • Additional properties: You can also customize advanced parameters such as effective mass tensors, band deformation potentials, and degeneracies, listed in the Material documentation.

Furthermore, in QTCAD Atoms, the tight-binding and Keating valence force-field model parameters can also be modified in the Atoms Materials module.

2. Does QTCAD support magnetic fields?

QTCAD can simulate the effects of magnetic fields on electron or hole spins. However, it does not solve for magnetic fields from a given magnetization distribution.

To obtain magnetic field maps for use in QTCAD simulations (e.g., for spin qubit devices), users may utilize other software tools such as:

These tools can generate magnetic field distributions that can be imported into QTCAD for further analysis.

3. How Can I include strain effects in QTCAD simulations?

QTCAD’s Schrödinger solver supports strain effects using the Bir-Pikus formalism (see Strain theory). Since QTCAD uses \(\mathbf{k}\cdot\mathbf{p}\) or effective-mass theory under the envelope function approximation, it focuses on band edges of parabolic bands rather than full band structures.

To include strain in your QTCAD simulation:

  • Deformation potentials: Strain effects are parametrized by deformation potentials set via material attributes:

    • cond_band_def_pot (for conduction band / electron calculations)

    • vlnce_band_def_pot (for valence band / hole calculations)

    These can be adjusted using the qtcad.device.materials.Material.set_param() method. See the qtcad.device.materials.Material documentation for details.

  • Strain tensor: The strain tensor can be applied to the device using the qtcad.device.Device.set_strain method. This method accepts:

    • A uniform strain tensor as a 3×3 NumPy array

    • A spatially varying strain tensor as a function returning a 3×3 array

    For examples and usage, refer to this tutorial.

Furthermore, our atomisitc modeling capabilities allow for strain effects to be simulated at the atomic level, providing a more detailed understanding of strain-induced phenomena in quantum devices.

If you have further questions, reach out to us by:

4. How to eliminate unphysical charge in my device?

Problem description In typical gated quantum dot devices, electrons or holes are confined within a host semiconductor (such as silicon or gallium arsenide), while a wider bandgap semiconductor or oxide forms a potential barrier separating the confined carriers from electrostatic gates used to tune the quantum dots. When strong gate potentials are applied, the conduction band edge (for electrons) or valence band edge (for holes) near the gate surfaces may cross the Fermi level, causing the non-linear Poisson solver to assign a significant charge density in that region. However, this charge accumulation may be unphysical at cryogenic temperatures because carriers cannot tunnel through sufficiently wide potential barriers from doped regions to the gate surfaces. This unphysical charge accumulation can prevent convergence of the Poisson solver.

Illustration

  1. Generate the device mesh by running:

python qtcad/examples/practical_application/1-devicegen.py

  1. Modify the applied gate potentials in Practical Application #2 example by setting:

# Applied potentials
Vtop_1 = -1.0
Vtop_2 = -1.0
Vtop_3 = -1.0
Vbottom_1 = -1.0
Vbottom_2 = -1.0
Vbottom_3 = -1.0

#(Changing from -0.5 V to -1.0 V)

  1. Instead of running the Schrödinger solver, plot the band diagram and hole density below top gate 2:

# Plot band diagram along a vertical linecut below top gate 2
x, y, z = d.mesh.glob_nodes.T
an.plot_bands(d, (235e-9, 500e-9, np.max(z)), (235e-9, 500e-9, np.min(z)))

# Plot hole density along the same linecut
an.plot_linecut(mesh, d.p/1e6, (235e-9, 500e-9, np.max(z)), (235e-9, 500e-9, np.min(z)),
    title="Hole density (cm$^{-3}$)")

  1. Observe that the Poisson solver error saturates (around 0.02 V after ~18 iterations), failing to converge.

  2. The band diagram shows strong valence band edge bending above the Fermi level near the gate interface, causing large, likely unphysical hole accumulation detrimental to solver convergence.

Solution: Suppressing Unwanted Charges

Immediately after defining device regions, declare the cap region as an ideal insulator by running:

# Set the cap to be an ideal insulator, in which:
## 1) Classical charge density is forced to zero
## 2) Carrier wave functions vanish in the Schrödinger solver
d.new_insulator("cap")

In this ideal insulator region:

  • Classical charge density is set exactly to zero

  • Carrier wave functions are forced to zero

This prevents charge accumulation near the gate interface, enabling the Poisson solver to converge. Band bending disappears near the gate, and hole density plots confirm the absence of unphysical charge accumulation.

Example code snippet

import pathlib
import numpy as np

from device import constants as ct
from device.mesh3d import Mesh, SubMesh
from device import materials as mt
from device import Device
from device.solver_params import SolverParams
from device.poisson import Solver as PoissonSolver
from device import analysis as an

# Load 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 and set temperature
d = Device(mesh)
d.set_temperature(0.1)

# Define doping and material regions
doping = 3e18 * 1e6  # In SI units

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)

# Enable ideal insulator behavior for the cap region to suppress charge
# Uncomment the following line to activate:
# d.new_insulator("cap")

# Align bands using GaAs as reference
d.align_bands(mt.GaAs)

# Set applied potentials (strong gate voltages)
Vtop_1 = -1.0
Vtop_2 = -1.0
Vtop_3 = -1.0
Vbottom_1 = -1.0
Vbottom_2 = -1.0
Vbottom_3 = -1.0

barrier = 0.834 * ct.e  # Schottky barrier height (n-type)

# Define 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")

# Configure and run non-linear Poisson solver
solver_params = SolverParams({"tol": 1e-3, "maxiter": 100}, problem="poisson")
poisson_solver = PoissonSolver(d, solver_params=solver_params)

# Create submesh for quantum dot region and set dot region
submesh = SubMesh(mesh, ["spacer_dot", "two_deg_dot", "substrate_dot"])
d.set_dot_region(submesh)

# Solve Poisson equation
poisson_solver.solve()

# Plot results
x, y, z = d.mesh.glob_nodes.T
an.plot_bands(d, (235e-9, 500e-9, np.max(z)), (235e-9, 500e-9, np.min(z)))
an.plot_linecut(mesh, d.p / 1e6, (235e-9, 500e-9, np.max(z)), (235e-9, 500e-9, np.min(z)),
                title="Hole density (cm$^{-3}$)")

5. I encountered an error when running the NEGF solver. What should I do?

If you encounter an error when running the NEGF solver, it is often due to issues with the mesh or boundary conditions. Here are some steps to troubleshoot:

  • Ensure that the mesh is extruded along the channel direction and that the boundaries are properly defined.

  • Verify that the boundary conditions are correctly set for the NEGF solver. In particular, two quantum leads—source and drain—must be specified using the Device.new_quantum_lead_bnd() method. These must be defined on the main Device object before constructing the NEGF SubDevice.

For more details, see the NEGF solver documentation: NEGF solver API documentation

6. My atomistic tight-binding simulations ran out of memory. What can I do?

Atomistic tight-binding simulations can require large amounts of memory, especially for large atomic structures or large basis sets. The QTCAD Atoms Schrödinger Solver supports two different algorithms with different memory footprints. The algorithm is selected through the memory_mode attribute of the solver’s SolverParams:

  • the default algorithm, set through memory_mode="robust", which provides the highest numerical stability, and is recommended for problems that may be numerically challenging;

  • the low-memory algorithm, set through memory_mode="low", which uses around 40% less memory for large atomic structures, and may be used when the memory available on the system is insufficient to run the default algorithm.

The low-memory algorithm is less numerically robust than the default algorithm, in the sense that it may occasionally produce unphysical states. Such situations can occur, for example, in some systems with very large atomic structures or rapidly varying confinement potentials. To circumvent this issue, two strategies are implemented in QTCAD:

  • a small random perturbation is added to the tight-binding Hamiltonian, which improves the numerical robustness of the algorithm; its magnitude is controlled by the noise attribute of SolverParams;

  • the solver’s internal data structures are checked at a relatively early stage of the computation to determine whether unphysical states are likely to be produced; if so, the solver is automatically restarted with a new random perturbation of the tight-binding Hamiltonian; the maximum number of restarts is set by the max_retries attribute of SolverParams.

The default values of noise and max_retries should ensure that physical states are produced in most cases. If not, numerical stability may be improved by:

  • reducing computational complexity, for example by reducing the size of the simulated atomic structures or by using an atomistic tight-binding basis set with fewer orbitals per atom;

  • ensuring that the potential applied on the atomic structure properly confines charge carriers;

  • setting the energy_target input of the solver’s solve method within the bandgap and sufficiently far (ideally several hundred meV away) from the band edge of interest (valence band maximum for holes and conduction band minimum for electrons);

  • reducing the number of threads through the num_threads attribute of SolverParams.

The QTCAD Atoms Schrödinger Solver also supports writing internal data structures to disk to save memory, as controlled by the max_memory attribute of SolverParams. In this case, the computational bottleneck becomes writing and reading data to and from the disk. This can significantly increase runtime, especially on slow storage devices, but allows for simulations that could not otherwise fit in memory.