1. Multiscale simulations of a SiGe–Si–SiGe quantum dot—Part 1: Finite-element method for device electrostatics

1.1. Requirements

1.1.1. Software components

  • QTCAD

  • Gmsh

  • Paraview

1.1.2. Geometry file

  • qtcad/examples/tutorials/meshes/multiscale.geo

1.1.3. Python script

  • qtcad/examples/tutorials/multiscale_1.py

1.1.4. References

1.2. Briefing

Spin qubit devices typically incorporate a quantum dot (QD) formed in a semiconductor heterostructure and controlled by a set of gate electrodes. These electrodes and the QD are separated by one or more layers of semiconductor and/or insulating materials that help to confine electrons or holes within the QD. As a result, the overall physical dimensions of spin qubit devices generally range from about one hundred nanometers up to a few micrometers. In contrast, the wavefunctions of the confined electrons or holes often extend over only a few nanometers to a few tens of nanometers along at least one confinement direction.

At these small length scales, atomistic effects can significantly influence device behavior. For instance, random alloy fluctuations and surface roughness can impact valley splitting in Si and Ge quantum dots [PWLK+22], which is an important figure of merit for device performance. Moreover, this random atomistic disorder can lead to substantial device-to-device variability, posing a key challenge for the scalability of spin-qubit-based quantum computing technologies. To accurately capture these phenomena, multiscale simulation approaches are essential.

In this series of tutorials, we apply a multiscale simulation approach to model quantum control in a SiGe–Si–SiGe QD. First, in the present tutorial, device electrostatics are resolved over the full spin qubit device scale using the finite element method (FEM), specifically the linear Poisson solver. The second tutorial will focus on the atomistic modeling of the QD using the tight-binding (TB) model, while the third tutorial will focus on harnessing the FEM and TB simulation results to model quantum control of the QD.

The goals of this tutorial are as follows.

  • Resolve the device electrostatics and solve the Schrödinger equation using the FEM.

  • Extract the electric potential in the region of the QD. This will be used to set the potential energy term for the TB Hamiltonian in the next tutorial.

  • Estimate the size of the QD ground state using the FEM. This will be used to construct an atomic structure of appropriate size for the TB simulations of the next tutorial.

  • Compute the derivative of the confinement potential with respect to the barrier gate voltage. This will be used to compute the Rabi frequency of qubit transitions induced by an oscillating barrier gate voltage.

Note

This tutorial involves a finite-element mesh containing around one million nodes. We expect the simulation presented here to run in a few minutes on a relatively advanced desktop workstation (as described in Software requirements). Less powerful machine may have memory constraints when running this simulation.

1.3. Device geometry

The device geometry is defined in the Gmsh geometry file qtcad/examples/tutorials/meshes/multiscale.geo. The heterostructure stack of the device consists of the following layers:

  • a \(5\;\mathrm{nm}\)-thick \(\mathrm{HfO}_2\) layer;

  • a \(30\;\mathrm{nm}\)-thick SiGe “spacer” layer;

  • a \(3\;\mathrm{nm}\)-thick Si “well” layer in which resides the QD;

  • a \(100\;\mathrm{nm}\)-thick SiGe “buffer” layer.

Together with the SiGe spacer layer, the SiGe buffer layer forms the SiGe barrier that confines electrons within the Si well layer. This heterostructure stack closely mimics that of the device described in Ref. [HHC+21].

The appropriate .msh and files may be generated by executing the following command in the qtcad/examples/tutorials/meshes folder.

gmsh - multiscale.geo

1.4. Setting up the device, solving the Poisson equation, and solving the Schrödinger equation

1.4.2. Loading the mesh and defining the device

Next, we load the mesh of the device.

# Load the mesh
script_dir = pathlib.Path(__file__).parent.resolve()
path_mesh = str(script_dir / "meshes" / "multiscale.msh")
scaling = 1e-9
mesh = Mesh(scaling, path_mesh)

This mesh may then be used to define the device.

# Create device
d = Device(mesh, conf_carriers="e")
d.set_temperature(1.6)

Note that we set the temperature to \(1.6\;\mathrm{K}\) [HHC+21].

Our TB simulations assume a valence-band offset of \(0.68\;\mathrm{eV}\) between Si and Ge, a value which slightly differs from that assumed by default in the FEM simulations. To ensure consistency between the FEM and TB simulations, we adjust the valence-band offset in the FEM simulation.

# Ensure band alignment between Ge and Si matches that assumed in TB
# calculations
VBO = 0.68 * ct.e # valence band offset from Niquet et al.
mt.Ge.set_param("chi", mt.Si.chi - VBO + mt.Si.Eg - mt.Ge.Eg)

We then define the SiGe alloy to be incorporated in the device. Specifically, we consider an alloy with a \(30\%\) Ge concentration and a \(70\%\) Si concentration.

# SiGe alloy
x = 0.3 # Ge fraction in SiGe alloy
SiGe = mt.SiGe
SiGe.set_alloy_composition(x)

From here, the various layers of the heterostructure may be defined. We note that the device is completely undoped.

# Define material stack in heterostructure
d.new_region("oxide", mt.HfO2)
d.new_region("cap_top", SiGe)
d.new_region("cap_bot", SiGe)
d.new_region("cap_qd", SiGe)
d.new_region("well", mt.Si)
d.new_region("well_qd", mt.Si)
d.new_region("buffer_top", SiGe)
d.new_region("buffer_qd", SiGe)
d.new_region("buffer_bot", SiGe)

Finally, we define the gate boundary conditions for the global top gate and the barrier gates. Here, we assume that the workfunction of the gate metal (\(\mathrm{TiN} \left(100\right)\)) is \(2.9\;\mathrm{eV}\) [CC20].

# Applied potentials
V_confinement = -1.9
V_plunger = -0.4

# Set up boundary conditions
Ew = 2.9 * ct.e # workfunction of TiN (100)
d.new_gate_bnd("gate_confinement", phi=V_confinement, Ew=Ew)
d.new_gate_bnd("gate_plunger", phi=V_plunger, Ew=Ew)

Two gate electrodes are present in the device. First, a confinement gate, on which a voltage of \(-1.9\;\mathrm{V}\) is applied, is used to confine electrons within the QD. Second, a plunger gate, on which a voltage of \(-0.4\;\mathrm{V}\) is applied, is used ensure the energies of the ground state and first few excited states lie below the Fermi energy.

1.4.3. Solving the linear Poisson equation

Next, we solve the linear Poisson equation.

# Set up and run adaptive Poisson solver
poisson_solver = PoissonSolver(d=d)
poisson_solver.solve()

We now save various quantities of interest. First, the electric potential is saved into a .hdf5 file for use in the next tutorial of this series.

# Save the Poisson equation solution in .hdf5 format
path_phi_hdf5 = str(script_dir / "output" / "multiscale_phi.hdf5")
arrays_dict = {"phi": d.phi}
io.save(path_phi_hdf5, arrays_dict)

Similarly, the conduction band edge is saved as .vtu files for visualization in Paraview.

# Save the conduction band edge and electron density in .vtu format
path_CB_vtu = str(script_dir / "output" / "multiscale_CB.vtu")
CB = d.cond_band_edge() / ct.e
out_dict = {"EC (eV)": CB}
io.save(path_CB_vtu, out_dict, d.mesh)

We also save slices of these quantities as .png files. The slices are taken normal to the heterostructure growth axis, \(1\;\mathrm{nm}\) below the interface between the SiGe spacer and the Si QD.

# Save slices of the conduction band edge and electron density
normal = (0., 0., 1.) # the plane of the slice is normal to the z axis
# the slice is just below the interface between the SiGe spacer and the Si QW
origin = (0., 0., -36. * scaling)
label_CB = "Conduction band edge (eV)"
path_CB_slice = str(script_dir / "output" / "multiscale_CB_slice.png")
plot_slice(d.mesh, CB, normal=normal, origin=origin, cb_axis_label=label_CB,
    path=path_CB_slice, show_figure=False)
Slice of the conduction band edge in the QW and QD.

Fig. 1.4.3 Slice of the conduction band edge in the QW and QD.

Here, we see that the barrier gates effectively confine electrons within the QD. In addition, the conduction band edge in the center of the QD is slightly below \(0\;\mathrm{eV}\), a value which is low enough to ensure that the QD ground state energy is below \(0\;\mathrm{eV}\) (as we will see later).

1.4.4. Solving the Schrödinger equation

Having resolved the device electrostatics, we may now solve the Schrödinger equation. To do so, we first compute the potential energy in the device.

# Get the potential energy from the band edge for usage in the Schrodinger
# solver
d.set_V_from_phi()

Since the QD is only a small part of the device, we restrict the Schrödinger solver to the QD region. This is done by defining an appropriate subdevice.

# Create a subdevice for the dot region
qd = ["cap_qd", "well_qd", "buffer_qd"] # quantum dot regions
submesh = SubMesh(d.mesh, qd)
subdvc = SubDevice(d, submesh)

Note that the dot region includes thin layers (\(2\;\mathrm{nm}\) thick) of the SiGe spacer and SiGe buffer in the direct vicinity of the Si layer. Indeed, in general, wavefunctions have tails that slightly penetrate into the barrier material.

We also save the potential energy in the QD for future reference.

# Save confinement potential in subdevice for later use
subdvc.set_V_from_phi()
V_0 = subdvc.get_Vconf()

We may now solve the Schrödinger equation.

# Set up and run Schrodinger solver
schrod_solver = SchrodingerSolver(d=subdvc)
schrod_solver.solve()

The eigenenergies and wavefunction statistics may now be printed using the following commands.

# Print energies and analyze ground state
subdvc.print_energies()
analyze_dot(subdvc, eigenstate=0, verbose=True)

This should result in the following output.

Energy levels (eV)
[-0.01577973 -0.00597006 -0.00596592 0.00390579 0.00390866 0.00397118 0.01384522 0.01387222 0.01390747 0.01398089]
--------------------------------------------------------------------------------
Geometric properties of the quantum dot:
--------------------------------------------------------------------------------
position        :  [-2.68226488e-12  1.56254636e-12 -3.56817108e-08]
std             :  [4.80501884e-09 4.80516151e-09 1.09137233e-09]
size            :  [1.92200754e-08 1.92206460e-08 4.36548931e-09]
--------------------------------------------------------------------------------

Note that only three states are below the Fermi energy. In addition, the ground state has a size of around \(19.2\;\mathrm{nm}\) along the \(x\) and \(y\) axes, and \(4.4\;\mathrm{nm}\) along the \(z\) axis (the heterostructure growth axis). This information will be relevant in the next tutorial to construct an atomic structure for the QD.

1.4.5. Saving the solution of the Schrödinger equation

We save the wavefunctions of the ground state, first excited state, and second excited state in a .vtu file.

# Save the first few eigenstates in .vtu format
path_wf_vtu = str(script_dir / "output" / "multiscale_wf.vtu")
wf0 = subdvc.eigenfunctions[:,0]
wf1 = subdvc.eigenfunctions[:,1]
wf2 = subdvc.eigenfunctions[:,2]
out_dict = {"Ground state": wf0, "1st excited": wf1, "2nd excited": wf2}
io.save(path_wf_vtu, out_dict, submesh)

In addition, we save a slice of the ground state wavefunction.

# Save slice of the ground state
path_wf0_slice = str(script_dir / "output" / "multiscale_wf0_slice.png")
label_wf = "Ground state wavefunction"
plot_slice(submesh, wf0, normal=normal, origin=origin,
    cb_axis_label=label_wf, title="Ground state", path=path_wf0_slice,
    show_figure=False)
Slice of the ground state wavefunction in the QD.

Fig. 1.4.4 Slice of the ground state wavefunction in the QD.

1.5. Computing the derivative of the confinement potential with respect to the barrier gate voltage

Finally, we compute the derivative of the confinement potential with respect to the barrier gate voltage. The derivative calculation is performed numerically. Specifically, we increment the barrier gate voltage by a small increment of \(0.01\;\mathrm{V}\) and solve the linear Poisson equation again. The confinement potential obtained via this procedure is then stored in a variable V_1, which is subtracted from V_0 and divided by the small gate voltage increment V_epsilon to obtain the derivative.

# Compute derivative of quantum dot confinement potential with respect to
# barrier gate voltage
V_epsilon = 0.01 # small increment in voltage to numerically compute derivative
V_confinement += V_epsilon
d.new_gate_bnd("gate_confinement", phi=V_confinement, Ew=Ew)
d.new_gate_bnd("gate_plunger", phi=V_plunger, Ew=Ew)
poisson_solver.solve()
subdvc = SubDevice(d, submesh)
subdvc.set_V_from_phi()
V_1 = subdvc.get_Vconf()
derivative = (V_1 - V_0) / V_epsilon
path_der_hdf5 = str(script_dir / "output" / "multiscale_der.hdf5")
arrays_dict = {"der": derivative}
io.save(path_der_hdf5, arrays_dict) # save the derivative in .hdf5 format
path_der_slice = str(script_dir / "output" / "multiscale_der_slice.png")
plot_slice(submesh, derivative/ct.e, normal=normal, origin=origin,
    cb_axis_label="Derivative (meV/V)", path=path_der_slice,
    show_figure=False) # save slice of the derivative

Above, note that we have saved the derivative in the .hdf5 and .png formats.

Slice of the derivative of the confinement potential with respect to the barrier gate voltage in the QD.

Fig. 1.5.4 Slice of the derivative of the confinement potential with respect to the barrier gate voltage in the QD.

We note that the slice shown in the figure above is restricted to the QD region (namely, it excludes the QW surrounding the QD). We see that the derivative is around \(-0.475\;\mathrm{meV/V}\) in the center of the QD, while it is around \(-0.650\;\mathrm{meV/V}\) on the periphery of the QD. This is expected, since the barrier gates are physically closer to the periphery of the QD than to its center, thereby leading to stronger gate control on the periphery.

1.6. Full code

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

import pathlib
from qtcad.device import constants as ct
from qtcad.device import materials as mt
from qtcad.device.mesh3d import Mesh, SubMesh
from qtcad.device import Device, SubDevice
from qtcad.device.poisson_linear import Solver as PoissonSolver
from qtcad.device.schrodinger import Solver as SchrodingerSolver
from qtcad.device import io
from qtcad.device.analysis import plot_slice, analyze_dot

# Load the mesh
script_dir = pathlib.Path(__file__).parent.resolve()
path_mesh = str(script_dir / "meshes" / "multiscale.msh")
scaling = 1e-9
mesh = Mesh(scaling, path_mesh)

# Create device
d = Device(mesh, conf_carriers="e")
d.set_temperature(1.6)

# Ensure band alignment between Ge and Si matches that assumed in TB
# calculations
VBO = 0.68 * ct.e # valence band offset from Niquet et al.
mt.Ge.set_param("chi", mt.Si.chi - VBO + mt.Si.Eg - mt.Ge.Eg)

# SiGe alloy
x = 0.3 # Ge fraction in SiGe alloy
SiGe = mt.SiGe
SiGe.set_alloy_composition(x)

# Define material stack in heterostructure
d.new_region("oxide", mt.HfO2)
d.new_region("cap_top", SiGe)
d.new_region("cap_bot", SiGe)
d.new_region("cap_qd", SiGe)
d.new_region("well", mt.Si)
d.new_region("well_qd", mt.Si)
d.new_region("buffer_top", SiGe)
d.new_region("buffer_qd", SiGe)
d.new_region("buffer_bot", SiGe)

# Applied potentials
V_confinement = -1.9
V_plunger = -0.4

# Set up boundary conditions
Ew = 2.9 * ct.e # workfunction of TiN (100)
d.new_gate_bnd("gate_confinement", phi=V_confinement, Ew=Ew)
d.new_gate_bnd("gate_plunger", phi=V_plunger, Ew=Ew)

# Set up and run adaptive Poisson solver
poisson_solver = PoissonSolver(d=d)
poisson_solver.solve()

# Save the Poisson equation solution in .hdf5 format
path_phi_hdf5 = str(script_dir / "output" / "multiscale_phi.hdf5")
arrays_dict = {"phi": d.phi}
io.save(path_phi_hdf5, arrays_dict)

# Save the conduction band edge and electron density in .vtu format
path_CB_vtu = str(script_dir / "output" / "multiscale_CB.vtu")
CB = d.cond_band_edge() / ct.e
out_dict = {"EC (eV)": CB}
io.save(path_CB_vtu, out_dict, d.mesh)

# Save slices of the conduction band edge and electron density
normal = (0., 0., 1.) # the plane of the slice is normal to the z axis
# the slice is just below the interface between the SiGe spacer and the Si QW
origin = (0., 0., -36. * scaling)
label_CB = "Conduction band edge (eV)"
path_CB_slice = str(script_dir / "output" / "multiscale_CB_slice.png")
plot_slice(d.mesh, CB, normal=normal, origin=origin, cb_axis_label=label_CB,
    path=path_CB_slice, show_figure=False)

# Get the potential energy from the band edge for usage in the Schrodinger
# solver
d.set_V_from_phi()

# Create a subdevice for the dot region
qd = ["cap_qd", "well_qd", "buffer_qd"] # quantum dot regions
submesh = SubMesh(d.mesh, qd)
subdvc = SubDevice(d, submesh)

# Save confinement potential in subdevice for later use
subdvc.set_V_from_phi()
V_0 = subdvc.get_Vconf()

# Set up and run Schrodinger solver
schrod_solver = SchrodingerSolver(d=subdvc)
schrod_solver.solve()

# Print energies and analyze ground state
subdvc.print_energies()
analyze_dot(subdvc, eigenstate=0, verbose=True)

# Save the first few eigenstates in .vtu format
path_wf_vtu = str(script_dir / "output" / "multiscale_wf.vtu")
wf0 = subdvc.eigenfunctions[:,0]
wf1 = subdvc.eigenfunctions[:,1]
wf2 = subdvc.eigenfunctions[:,2]
out_dict = {"Ground state": wf0, "1st excited": wf1, "2nd excited": wf2}
io.save(path_wf_vtu, out_dict, submesh)

# Save slice of the ground state
path_wf0_slice = str(script_dir / "output" / "multiscale_wf0_slice.png")
label_wf = "Ground state wavefunction"
plot_slice(submesh, wf0, normal=normal, origin=origin,
    cb_axis_label=label_wf, title="Ground state", path=path_wf0_slice,
    show_figure=False)

# Compute derivative of quantum dot confinement potential with respect to
# barrier gate voltage
V_epsilon = 0.01 # small increment in voltage to numerically compute derivative
V_confinement += V_epsilon
d.new_gate_bnd("gate_confinement", phi=V_confinement, Ew=Ew)
d.new_gate_bnd("gate_plunger", phi=V_plunger, Ew=Ew)
poisson_solver.solve()
subdvc = SubDevice(d, submesh)
subdvc.set_V_from_phi()
V_1 = subdvc.get_Vconf()
derivative = (V_1 - V_0) / V_epsilon
path_der_hdf5 = str(script_dir / "output" / "multiscale_der.hdf5")
arrays_dict = {"der": derivative}
io.save(path_der_hdf5, arrays_dict) # save the derivative in .hdf5 format
path_der_slice = str(script_dir / "output" / "multiscale_der_slice.png")
plot_slice(submesh, derivative/ct.e, normal=normal, origin=origin,
    cb_axis_label="Derivative (meV/V)", path=path_der_slice,
    show_figure=False) # save slice of the derivative