3. Multiscale simulations of a SiGe–Si–SiGe quantum dot—Part 3: Quantum control

3.1. Requirements

3.1.1. Software components

  • QTCAD

3.1.2. Python script

  • qtcad/examples/tutorials/multiscale_3.py

3.1.3. References

3.2. Briefing

In the previous tutorial of this tutorial series, we harnessed finite-element method (FEM) simulations and tight-binding (TB) simulations to compute three parameters of a SiGe–Si–SiGe quantum dot (QD):

  • the valley splitting \(\Delta\), namely the difference between the energy levels of the two lowest-energy states of the QD, which we denote as \(\left|\nu_1\right>\) and \(\left|\nu_2\right>\);

  • the spin–orbit-coupling (SOC) Hamiltonian matrix element \(C_{\nu_1\nu_2}:=\left<\nu_2\uparrow\right|H_{\text{SOC}}\left|\nu_1\downarrow\right>\), where \(H_{\text{SOC}}\) is the SOC term of the TB Hamiltonian and \(\uparrow,\downarrow\) denote the spin-up and spin-down states;

  • the confinement potential derivative matrix element \(D_{\nu_1\nu_2}:=\left<\nu_1\right|D_{\text{BG}}\left|\nu_2\right>\), where \(D_{\text{BG}}\) is the derivative of the confinement potential with respect to barrier gate voltage.

These three values may be used to parametrize a quantum control scheme proposed in Ref. [BN18]. Loosely speaking, in this scheme, the valley states \(\left|\nu_1\right>\) and \(\left|\nu_2\right>\) are each split into two valley–spin states by applying a constant, uniform magnetic field. Near a critical magnetic field, the states \(\left|\nu_1\uparrow\right>\) and \(\left|\nu_2\downarrow\right>\) are energetically close to each other and relatively strongly coupled to each other by the SOC term of the Hamiltonian \(H_{\text{SOC}}\), leading to mixed eigenstates exhibiting an anticrossing. Rabi oscillations between these mixed states may be driven by applying a time-dependendent voltage on the barrier gate.

In this tutorial, we will not perform detailed derivations and explanations of this quantum control scheme. Rather, we highlight that FEM and TB simulations can be used to model practical quantum control schemes. Readers interested in theoretical details may refer to Ref. [BN18].

3.3. Quantum dot eigenstates

We start by importing the required modules and libraries, which we have all seen in previous tutorials.

import pathlib
import numpy as np
from matplotlib import pyplot as plt
import qtcad.device.constants as ct

Recall that we saved the three quantum control model parameters to a .txt file in the previous tutorial. We now load these parameters to a Numpy array params.

# Load the quantum control parameters from the TB simulation
script_dir = pathlib.Path(__file__).parent.resolve()
path_params = str(script_dir / "output" / "multiscale_params.txt")
params = np.loadtxt(path_params, dtype=complex)
V = np.real(params[0]) # valley splitting
C = params[1] # matrix element of SOC Hamiltonian
D = np.real(params[2]) # matrix element of derivative of potential

Note that the parameters are stored as complex numbers; however, only \(C_{\nu_1\nu_2}\) has a non-zero imaginary part. We thus take the real parts of the appropriate parameters in the .txt file.

Next, we define four variables for convenience.

# Related model parameters
E1 = - V/2 # energy level of first valley state
E2 = V/2 # energy level of second valley state
g = 2.0 # gyromagnetic factor
BA = V / (g * ct.muB) # anticrossing magnetic field

The first and second, E1 and E2, are the energy levels of the two valley states \(\left|\nu_1\right>\) and \(\left|\nu_2\right>\) in the absence of a magnetic field, which we respectively set to \(E_1=-\frac{\Delta}{2}\) and \(E_2=\frac{\Delta}{2}\). These energies are expressed in a gauge where the zero of energy is set to the middle of the valley splitting. The third, g, is the electron gyromagnetic factor, which we set to be identically equal to \(g=2\) for simplicity. The fourth, BA, is the critical magnetic field \(B_A\) at which the anticrossing between the valley–spin states occurs, as we will see.

To investigate the behaviour of the QD eigenstates as a function of magnetic field \(B\), we instantiate an array B of sampled magnetic fields.

# Sampled values of amplitude of magnetic field
Bmax = 1.5 * BA
Bnum = 1000
B = np.linspace(0.0, Bmax, Bnum)

For \(B\ll B_A\) and \(B\gg B_A\), the Zeeman-split QD eigenstates are \(\left|\nu_n\sigma\right>\) with eigenenergies \(E_n + \frac{1}{2} g\mu_B\left<\sigma\right>\), where \(n=1,2\), \(\sigma=\uparrow,\downarrow\), \(\left<\sigma\right>=+1,-1\), and \(\mu_B\) is the Bohr magneton; the spin is quantized along the direction of the magnetic field. Parenthetically, note that the second and third-lowest energy states of the QD predicted by this formula are degenerate at magnetic field \(B_A=\frac{\Delta}{g\mu_B}\). For \(B\approx B_A\), the authors of Ref. [BN18] predict the following second and third-lowest eigenenergies:

\[E_{\pm} = \frac{1}{2}\left(E_1 + E_2\right) \pm \frac{1}{2}\sqrt{\left(\Delta-g\mu_B B\right)^2 + 4\left|C_{\nu_1\nu_2}\right|^2}\,.\]

Note the dependence on \(C_{\nu_1\nu_2}\). From these formulas, we may compute the eigenenergies of the QD as a function of magnetic field.

# Parametrized energy levels of the quantum dot
E_QD_0 = E1 - 0.5 * g * ct.muB * B
E_QD_1 = 0.5 * (E1 + E2) - \
    0.5 * np.sqrt((V - g * ct.muB * B)**2 + 4 * np.abs(C)**2)
E_QD_2 = 0.5 * (E1 + E2) + \
    0.5 * np.sqrt((V - g * ct.muB * B)**2 + 4 * np.abs(C)**2)
E_QD_3 = E2 + 0.5 * g * ct.muB * B

These eigenenergies may then be plotted.

# Plot the energy levels
plt.plot(B, 1e3*E_QD_0/ct.e, label=r"$\left|0\right>$")
plt.plot(B, 1e3*E_QD_1/ct.e, label=r"$\left|1\right>$")
plt.plot(B, 1e3*E_QD_2/ct.e, label=r"$\left|2\right>$")
plt.plot(B, 1e3*E_QD_3/ct.e, label=r"$\left|3\right>$")
plt.xlabel("Magnetic field (T)")
plt.ylabel("Energy (meV)")
plt.legend()
plt.tight_layout()
plt.show()

This should result in the following figure.

Eigenenergies of the four lowest-energy states of the QD.

Fig. 3.3.1 Eigenenergies of the four lowest-energy states of the QD.

3.4. Rabi frequency

According to Ref. [BN18], the four QD eigenstates are given by:

\[\begin{split}\left|0\right> &= \beta^{\prime}\left|\nu_1\downarrow\right> - \alpha^{\prime\star}\left|\nu_2\uparrow\right>\,,\\ \left|1\right> &= \beta\left|\nu_1\uparrow\right> - \alpha^{\star}\left|\nu_2\downarrow\right>\,,\\ \left|2\right> &= \alpha\left|\nu_1\uparrow\right> + \beta\left|\nu_2\downarrow\right>\,,\\ \left|3\right> &= \alpha^{\prime}\left|\nu_1\downarrow\right> + \beta^{\prime}\left|\nu_2\uparrow\right>\,,\\\end{split}\]

where the coefficients \(\alpha,\beta,\alpha^{\prime},\beta^{\prime}\) are given by:

\[\begin{split}\alpha &= \frac{2C^{\star}_{\nu_1\nu_2}}{\sqrt{\varepsilon^2+4\left|C_{\nu_1\nu_2}\right|^2}}\,,\\ \beta &= \frac{\varepsilon}{\sqrt{\varepsilon^2+4\left|C_{\nu_1\nu_2}\right|^2}}\,,\\ \alpha^{\prime} &= -\frac{2C_{\nu_1\nu_2}}{\sqrt{\varepsilon^{\prime 2}+4\left|C_{\nu_1\nu_2}\right|^2}}\,,\\ \beta^{\prime} &= \frac{\varepsilon^{\prime}}{\sqrt{\varepsilon^{\prime 2}+4\left|C_{\nu_1\nu_2}\right|^2}}\,,\\\end{split}\]

with:

\[\begin{split}\varepsilon &= \Delta - g\mu_B B + \sqrt{\left(\Delta-g\mu_B B\right)^2 + 4\left|C_{\nu_1\nu_2}\right|^2}\,,\\ \varepsilon^{\prime} &= \Delta + g\mu_B B + \sqrt{\left(\Delta+g\mu_B B\right)^2 + 4\left|C_{\nu_1\nu_2}\right|^2}\,.\\\end{split}\]

Note that for \(B\ll B_A\), \(\left|1\right>\approx\left|\nu_1\uparrow\right>\); for \(B\gg B_A\), \(\left|1\right>\approx\left|\nu_2\downarrow\right>\). On the other hand, for all \(B\), \(\left|0\right>\approx\left|\nu_1\downarrow\right>\). Thus, for \(B\ll B_A\), the ground and first-excited state of the QD have the same valley character but distinct spin characters. On the other hand, for \(B\gg B_A\), they have the same spin character but distinct valley characters. The \(B\ll B_A\) regime may thus be referred as the “spin mode”, while the \(B\gg B_A\) regime may be referred as the “valley mode”.

The Rabi frequency \(f\) of oscillations between the \(\left|0\right>\) and \(\left|1\right>\) states induced by a time-dependent voltage on the barrier gate may then be expressed through:

\[hf = e \delta V_{\mathrm{BG}} \left|\alpha^{\star}\beta^{\prime}+\alpha^{\prime}\beta\right|\left|D_{\nu_1\nu_2}\right|\,,\]

where \(e\) is the elementary charge, \(h\) is Planck’s constant, and \(\delta V_{\mathrm{BG}}\) is the amplitude of the time-dependent voltage applied on the barrier gate.

To compute this Rabi frequency as a function of magnetic field, we harness the formulas above.

# Parametrized EDSR Rabi frequency
eps = V - g*ct.muB*B + np.sqrt((V-g*ct.muB*B)**2 + 4*np.abs(C)**2)
alpha = 2*np.conj(C) / np.sqrt(eps**2 + 4*np.abs(C)**2)
beta = eps / np.sqrt(eps**2 + 4*np.abs(C)**2)
epsp = V + g*ct.muB*B + np.sqrt((V+g*ct.muB*B)**2 + 4*np.abs(C)**2)
alphap = -np.conj(2*np.conj(C) / np.sqrt(epsp**2 + 4*np.abs(C)**2))
betap = epsp / np.sqrt(epsp**2 + 4*np.abs(C)**2)
VBG = 0.001 # amplitude of rf oscillations of barrier gate voltage
hf = VBG * np.abs(np.conj(alpha)*betap+alphap*beta) * np.abs(D)

We then plot the Rabi frequency as a function of magnetic field.

# Plot the Rabi frequency
plt.plot(B, hf/ct.h/1e6)
plt.xlabel("Magnetic field (T)")
plt.ylabel("Rabi frequency (MHz)")
plt.tight_layout()
plt.show()

This should result in the following figure.

Rabi frequency of the QD as a function of magnetic field.

Fig. 3.4.1 Rabi frequency of the QD as a function of magnetic field.

The ability to tune the QD between the “spin mode” and the “valley mode” offers the possiblity of encoding quantum information in a manner relatively robust to decoherence in the spin mode, while allowing for fast electrical manipulation in the valley mode. Overall, this tutorial thus illustrates how multiscale simulations may be used to model quantum control schemes.

3.5. Full code

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

import pathlib
import numpy as np
from qtcad.device import constants as ct
from qtcad.device.mesh3d import Mesh, SubMesh
from qtcad.device import Device
from qtcad.device import io
from qtcad.atoms.unit_cell import UnitCellZincblende
from qtcad.atoms import Atoms, SubAtoms
from qtcad.atoms.rough_surface import RoughSurface
from qtcad.atoms.keating import Solver as KeatingSolver
from qtcad.atoms.keating import SolverParams as KeatingSolverParams
from qtcad.atoms.schrodinger import Solver as SchrodingerSolver
from qtcad.atoms.schrodinger import SolverParams as SchrodingerSolverParams
from qtcad.atoms.analysis import analyze_dot, save_vtu, get_operator_element, \
    get_operator_from_fem

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

# Create device and set potential in device to result of FEM simulation
d = Device(mesh, conf_carriers="e")
path_phi_hdf5 = str(script_dir / "output" / "multiscale_phi.hdf5")
d.set_potential(io.load(path_phi_hdf5, "phi"))

# Create unit cells of Si and Ge
unit_cell_Si = UnitCellZincblende(lattice_constant=5.431*1e-10,
    atom_species_ws=np.array(["Si", "Si"]))
unit_cell_Ge = UnitCellZincblende(lattice_constant=5.658*1e-10,
    atom_species_ws=np.array(["Ge", "Ge"]))

# Plot band structures of Si and Ge
path_bs_Si = str(script_dir / "output" / "multiscale_bs_Si.png")
unit_cell_Si.plot_bandstructure(title="Si", path=path_bs_Si, show_figure=False)
path_bs_Ge = str(script_dir / "output" / "multiscale_bs_Ge.png")
unit_cell_Ge.plot_bandstructure(title="Ge", path=path_bs_Ge, show_figure=False)

# Create atomic structure for heterostructure
x = np.array([-10., 10.]) * 1e-9
y = np.array([-10., 10.]) * 1e-9
z = np.array([-50., -5.]) * 1e-9
atoms = Atoms(x=x, y=y, z=z, rng_seed=104)

# Create rough surface describing top and bottom surfaces of the well
well_top = RoughSurface(hurst=0.3, mean=-35.*1e-9, rms=5.*1e-10, rng_seed=0)
well_bot = RoughSurface(hurst=0.3, mean=-38.*1e-9, rms=5.*1e-10, rng_seed=1)

# Plot rough surface
x_plot = np.linspace(x[0], x[1], 100)
y_plot = np.linspace(y[0], y[1], 100)
path_rough_2d = str(script_dir / "output" / "multiscale_rough_2d.png")
well_top.plot(x=x_plot, y=y_plot, type="2D", path=path_rough_2d,
    show_figure=False)
path_rough_3d = str(script_dir / "output" / "multiscale_rough_3d.png")
well_top.plot(x=x_plot, y=y_plot, type="3D", path=path_rough_3d,
    show_figure=False)

# Define the material stack of the heterostructure
SiGe_alloy = np.array([["Si", 0.7], ["Ge", 0.3]])
atoms.new_layer(z_top=-5.*1e-9, z_bot=well_top, atom_species=SiGe_alloy) # cap
atoms.new_layer(z_top=well_bot, z_bot=z[0], atom_species=SiGe_alloy) # buffer

# Relax atomic coordinates
keating_solver_params = KeatingSolverParams()
keating_solver_params.verbose = True
keating_solver = KeatingSolver(atoms=atoms,
    solver_params=keating_solver_params)
keating_solver.solve()
path_xyz = str(script_dir / "output" / "multiscale.xyz")
atoms.save_xyz(path=path_xyz)

# Contruct restricted atomic structure for QD
x_QD = np.array([-8., 8.]) * 1e-9
y_QD = np.array([-8., 8.]) * 1e-9
z_QD = np.array([-40., -33.]) * 1e-9
atoms_QD = SubAtoms(parent=atoms, x=x_QD, y=y_QD, z=z_QD)

# Create Schrödinger equation solver within TB formalism
atoms_QD.set_potential_from_device(d=d)
schrodinger_solver_params = SchrodingerSolverParams()
schrodinger_solver_params.verbose = True
schrodinger_solver_params.num_states = 10
schrodinger_solver = SchrodingerSolver(atoms=atoms_QD,
    solver_params=schrodinger_solver_params)

# Solve the Schrödinger equation
energy_target = -0.1 * ct.e
schrodinger_solver.solve(energy_target=energy_target)
atoms_QD.print_energies()

# Compute valley splitting
valley_splitting = atoms_QD.energies[1] - atoms_QD.energies[0]
print("Valley splitting: " + str(valley_splitting/ct.e) + " eV")

# Store valley states in variables nu_1 and nu_2
nu_1 = atoms_QD.eigenfunctions[0]
nu_2 = atoms_QD.eigenfunctions[1]

# Compute modulus-square of projections of ground state on the
# spds* orbitals and on the chemical species in the system
# Also compute the geometric properties of the ground state wavefunction
print("Ground state (nu_1)")
nu_1.get_prob_on_orbitals(verbose=True)
nu_1.get_prob_on_chemical_species(atoms=atoms_QD, verbose=True)
analyze_dot(atoms=atoms_QD, psi=nu_1, verbose=True)

# Repeat these steps for the first excited state
print("First excited state (nu_2)")
nu_2.get_prob_on_orbitals(verbose=True)
nu_2.get_prob_on_chemical_species(atoms=atoms_QD, verbose=True)
analyze_dot(atoms=atoms_QD, psi=nu_2, verbose=True)

# Save modulus-square of projection of states on atoms to a .vtu file
prob_on_atoms_0 = nu_1.get_prob_on_atoms()
prob_on_atoms_1 = nu_2.get_prob_on_atoms()
prob_on_atoms_2 = atoms_QD.eigenfunctions[2].get_prob_on_atoms()
path_vtu = str(script_dir / "output" / "multiscale_TB_wf.vtu")
out_dict = {"Probability density (ground state)": prob_on_atoms_0,
    "Probability density (1st excited state)": prob_on_atoms_1,
    "Probability density (2nd excited state)": prob_on_atoms_2}
save_vtu(atoms=atoms_QD, out_dict=out_dict, path=path_vtu)

# Compute the matrix element of the spin-orbit coupling Hamiltonian between
# the valley states with opposite spins
H_soc = schrodinger_solver.get_hamiltonian_soc()
nu_2_up = nu_2.add_spin_to_orbital_state(spin=True, in_place=False)
nu_1_down = nu_1.add_spin_to_orbital_state(spin=False, in_place=False)
C = get_operator_element(psi_l=nu_2_up, op=H_soc, psi_r=nu_1_down)
print("Matrix element of SOC Hamiltonian: " + str(C/ct.e) + " eV")

# Compute the matrix element of the derivative of the potential with respect to
# barrier gate voltage between the valley states
qd = ["cap_qd", "well_qd", "buffer_qd"] # quantum dot regions
submesh = SubMesh(mesh, qd)
path_der_hdf5 = str(script_dir / "output" / "multiscale_der.hdf5")
der = io.load(path_der_hdf5, "der")
Dfg = get_operator_from_fem(atoms=atoms_QD, tb_orbs=nu_1.tb_orbs,
    tb_spin=nu_1.tb_spin, mesh=submesh, quantity=der)
D = get_operator_element(psi_l=nu_1, op=Dfg, psi_r=nu_2)
print("Matrix element of derivative of potential with respect to barrier " + \
    "gate voltage: " + str(D/ct.e) + " eV/V")

# Save parameters of quantum control problem to .txt file
params = np.array([valley_splitting, C, D])
path_params = str(script_dir / "output" / "multiscale_params.txt")
np.savetxt(path_params, params)