18. Exchange coupling in a double quantum dot in FD-SOI—Part 2: Exact diagonalization

18.1. Requirements

18.1.1. Software components

  • QTCAD

  • Gmsh

18.1.2. Geometry file

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

18.1.3. Python script

  • qtcad/examples/tutorials/exchange_2.py

18.1.4. References

18.2. Briefing

In Exchange coupling in a double quantum dot in FD-SOI—Part 1: Perturbation theory, we considered a simple analytic expression of the exchange interaction strength which was derived from the Fermi–Hubbard Hamiltonian using perturbation theory. This expression was used to estimate the exchange interaction strength from two parameters that may be obtained from single-electron properties: the double dot tunnel coupling \(\Omega\) and the on-site Coulomb interaction strength \(U\).

In this tutorial, we go beyond this approximate perturbative treatment by using the exact diagonalization method to calculate the exchange interaction strength, as explained in the Exchange interaction through exact diagonalization theory section. Contrary to previous tutorials, which may be executed on a laptop within minutes, this tutorial is much more computationally expensive. In particular, computing the Coulomb integrals may be quite time consuming. Therefore, to execute the current tutorial, we recommend using a workstation that contains enough CPU cores to run at least 16 threads in parallel.

18.4. Setting up the device

Again, we start by setting up the paths for the input mesh and for the output files.

# Set up paths
script_dir = pathlib.Path(__file__).parent.resolve()
path_mesh = script_dir / "meshes"
path_out = script_dir / "output"

We then load the outcome of the perturbative calculation of Exchange coupling in a double quantum dot in FD-SOI—Part 1: Perturbation theory for future comparison with the exact diagonalization approach taken below.

# Load the perturbative value of the exchange coupling from a previous
# calculation
path_exchange_perturbation = path_out / "exchange_perturbation.txt"
exchange_perturbation = np.loadtxt(path_exchange_perturbation)

We then load the mesh, again using the final refined mesh from the adaptive calculation done in Tunnel coupling in a double quantum dot in FD-SOI—Part 1: Plunger gate tuning.

# Load the mesh
scaling = 1e-9
path_file = str(path_mesh / "refined_dqdfdsoi.msh2")
mesh = Mesh(scaling, path_file)

In this tutorial, as in Exchange coupling in a double quantum dot in FD-SOI—Part 1: Perturbation theory, we will use the bias configuration of Tunnel coupling in a double quantum dot in FD-SOI—Part 2: Tuning the barrier gate that gave the strongest tunnel coupling. This is to maximize the exchange interaction strength. Indeed, we expect larger values of exchange to be computed with lower relative error and to therefore be able to compute a reasonably accurate value of exchange even for coarse meshes and small basis sets. Therefore, due to the smaller computational resource requirements, this large exchange configuration is more appropriate for a tutorial.

We use the same gate bias configuration as in this tutorial, except that we do not include the numerical offset of plunger gate 2 that was computed in Tunnel coupling in a double quantum dot in FD-SOI—Part 1: Plunger gate tuning for the weakest tunnel coupling. This is not necessary here, since the artificial double dot detuning that results from mesh asymmetry is negligible compared with the strongest value of tunnel coupling found in Tunnel coupling in a double quantum dot in FD-SOI—Part 2: Tuning the barrier gate. For weaker tunneling, however, including a numerical offset may be necessary.

# Define the gate biases
back_gate_bias = -0.5
barrier_gate_1_bias = 0.5
plunger_gate_1_bias = 0.59
barrier_gate_2_bias = 0.57
plunger_gate_2_bias = 0.59
barrier_gate_3_bias = 0.5

After setting variable values for the gate biases, we define the Device using get_double_dot_fdsoi and create a SubMesh for the double quantum dot region.

# Define the device
dvc = get_double_dot_fdsoi(mesh, back_gate_bias, barrier_gate_1_bias,
    plunger_gate_1_bias, barrier_gate_2_bias, plunger_gate_2_bias,
    barrier_gate_3_bias)

# List of regions forming the double quantum dot region
dot_region_list = ["oxide_dot", "gate_oxide_dot", "buried_oxide_dot",
                    "channel_dot"]

# Create the submesh object for the dot region
submesh = SubMesh(dvc.mesh, dot_region_list)

18.5. Electrostatics and single-electron eigenstates

As usual, the next step is to compute the electric potential throughout the device using the non-linear Poisson solver. To do so, we load the potential that was obtained in Tunnel coupling in a double quantum dot in FD-SOI—Part 1: Plunger gate tuning into the device, configure the non-linear Poisson solver parameters, instantiate the Poisson Solver object, and run it using its solve method with the optional argument initialize set to False to use the electric potential stored in the device as an initial guess.

# Load the electric potential from the previous run to speed up Poisson
# calculation
phi = io.load(path_out / "tunnel_coupling_final_potential.hdf5",
    var_name="phi")
dvc.set_potential(phi)

# Configure the non-linear Poisson solver
params_poisson = PoissonSolverParams()
params_poisson.tol = 1e-3 # Convergence threshold (tolerance) for the error

# Instantiate Poisson solver
poisson_slv = PoissonSolver(dvc, solver_params=params_poisson)

# Solve Poisson's equation
poisson_slv.solve(initialize=False)

We next configure the single-electron Schrödinger solver, set the confinement potential from the electric potential obtained from the non-linear Poisson solution, create a SubDevice object for the double quantum dot region, instantiate a Schrödinger Solver object, call its solve method, and print the resulting eigenenergies.

# Configure the Schrödinger solver
params_schrod = SchrodingerSolverParams()
params_schrod.tol = 1e-6      # Tolerance on energies in eV

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

# Create a subdevice for the dot region
subdvc = SubDevice(dvc, submesh)

# Create a Schrodinger solver, and run it
print("Solving the single-electron Schrödinger equation")
schrod_solver = SchrodingerSolver(subdvc, solver_params=params_schrod)
schrod_solver.solve()
subdvc.print_energies()

We get

Energy levels (eV)
[0.0433603  0.04362118 0.04834174 0.04860898 0.05153343 0.0549579
0.0551495  0.05543713 0.05654119 0.06071358]

18.6. Exact diagonalization of the two-electron Hamiltonian

We first create a many-body SolverParams object to configure the exact diagonalization method, and use it in the instantiation of a many-body Solver object.

# Solving the many-body hamiltonian
# Instantiate many-body solver
manybody_solver_params = ManyBodySolverParams()
manybody_solver_params.n_degen = 2
manybody_solver_params.num_states = 8
manybody_solver_params.num_particles = [2]
manybody_solver_params.overlap = True
manybody_solver = ManyBodySolver(subdvc,
    solver_params=manybody_solver_params)

We use n_degen = 2 to account for spin degeneracy. In addition, we set num_states = 8, meaning that the first 8 spin-degenerate single-electron eigenstates are used as a basis set to construct the many-body Hamiltonian (bringing the total number of basis states to \(8\times2 = 16\)). We found that for the current gate bias configuration, smaller basis sets are not sufficiently accurate. Larger basis sets would be more accurate, but also lead to much longer computation times. In addition, we set num_particles = [2] to indicate that we are only interested in the two-particle subspace of the many-body Hamiltonian, and set overlap = True to indicate that we will keep all the terms in the four-dimensional Coulomb integrals matrix. Setting overlap to False would greatly speed up calculations, but we found this to also lead to results that are incorrect by several orders of magnitude.

Having set up the many-body solver, we run its solve method.

# solve the many-body hamiltonian
t0 = timer()
manybody_solver.solve()
print("Two-electron solution found in {} s".format(timer()-t0))

In addition, we use the default timer of the timeit module to measure the computation time. When writing this tutorial, we used roughly 70 threads of a 64-core AMD Ryzen Threadripper 3990X workstation, and completed this part of the calculation in approximately 45 minutes. Because this exchange calculation allows us to concentrate on the two-electron subspace, construction and diagonalization of the many-body Hamiltonian is not a computational bottleneck. In fact, almost all the computation time is spent in the calculation of Coulomb integrals, which are well parallelized and for which significant speedups can be otained simply by using a large number of cores.

18.7. Exchange interaction strength

Once the many-body problem is solved, we may readily calculate the exchange interaction strength.

We first extract the two-electron eigenenergies from the above calculations by calling the get_N_particle_subspace method of the Device class. After printing the first eight eigenenergies, we save all the energies in a text file.

# Get the two-electron eigenenergies
two_electron_energies = subdvc.get_N_particle_subspace(2).eigval/ct.e
print("Two-electron energies")
print(two_electron_energies[0:8])
np.savetxt(path_out/"exchange_two_electron_energies.txt", two_electron_energies)

For the first eight eigenenergies, we get:

Two-electron energies
[0.09096447 0.09096742 0.09096742 0.09096742 0.09568591 0.09568591
0.09568591 0.09569444]

As expected, we find a singlet (non-degenerate) ground state, and a triply-degenerate (triplet) first excited state.

The exchange interaction strength is given simply by the energy difference between the singlet and triplet states.

# Calculate the exchange interaction strength
exchange = two_electron_energies[1] - two_electron_energies[0]

Finally, we print the exchange interaction strength in MHz, and compare it with the value obtained from perturbation theory in the previous tutorial.

# Print the output, and compare it with perturbation theory
print ('Exchange coupling (exact diagonalization): {:.3f} MHz'\
    .format(exchange*ct.e/ct.h/1e6))
print ('Exchange coupling (perturbation theory): {:.3f} MHz'\
    .format(exchange_perturbation/ct.h/1e6))

We get:

Exchange coupling (exact diagonalization): 713.546 MHz
Exchange coupling (perturbation theory): 922.082 MHz

For the current gate bias configuration, solver parameters, and mesh, perturbation theory is in agreement with the exact diagonalization approach within a relative error of roughly 20 %.

18.8. Convergence analysis

Since the perturbative approach is approximate, we cannot expect it to agree with exact diagonalization, in general. Conditions under which the perturbative approach is expected to be accurate in principle are given in Perturbation theory of exchange coupling in a double quantum dot. Beyond these parametric regimes, for this method to be accurate, finite-element simulations employed to calculate the tunnel coupling and the Coulomb interaction strength must also be converged with respect to mesh refinements.

In addition, despite its name, the exact diagonalization approach is only perfectly accurate for infinitely fine meshes and infinitely large single-electron basis sets. In practice, an increasingly fine mesh and an increasingly large number of basis states are observed to be required for convergence as the exchange interaction strength is decreased.

Therefore, here, we show a convergence analysis of both the exact diagonalization and perturbative approaches with respect to the mesh characteristic length and the number of orbital basis states. The results of this convergence analysis are summarized in the plots below.

Exchange convergence with mesh.

Fig. 18.8.1 Convergence of the exchange interaction strength with mesh characteristic length. In these calculations, we set the num_states parameter of the many-body SolverParams object to 10 (compared with 8, above), and tune the h_dot parameter of the non-linear Poisson SolverParams object over five values between 0.8 and 0.4. This enables to refine the mesh within the quantum dot region only.

Exchange convergence with basis set.

Fig. 18.8.2 Convergence of the exchange interaction strength with the dimension of the single-electron basis set. In all these simulations, we set the h_dot parameter of the non-linear Poisson SolverParams object to 0.4 (compared with 0.8, above) to employ a relatively fine mesh in the double quantum dot region, and tune the num_states parameter of the many-body SolverParams object over 4 values between 4 and 10. For num_states = 4, we find that the ground state is triply degenerate and the first excited state is a singlet, which is a clear sign that the basis set is insufficiently large. Indeed, we expect the exchange interaction to lower the singlet state energy with respect to the triplet states in all realistic gated quantum dots.

We observe that agreement between perturbation theory and exact diagonalization is significantly better for fine meshes and large basis sets. However, while both techniques appear to converge to a fixed value, these individual converged values do not overlap with one another even for the largest meshes and basis sets. This is to be expected, since some errors are introduced when using the Fermi–Hubbard model and perturbation theory that are independent of the discretization of the problem.

18.9. Full code

__copyright__ = "Copyright 2023, Nanoacademic Technologies Inc."

import numpy as np
import pathlib
from timeit import default_timer as timer
from device import constants as ct
from device import io
from device.mesh3d import Mesh, SubMesh
from device import materials as mt
from device import Device, SubDevice
from device.poisson import SolverParams as PoissonSolverParams
from device.poisson import Solver as PoissonSolver
from device.schrodinger import SolverParams as SchrodingerSolverParams
from device.schrodinger import Solver as SchrodingerSolver
from device.many_body import Solver as ManyBodySolver
from device.many_body import SolverParams as ManyBodySolverParams
from examples.tutorials.double_dot_fdsoi import get_double_dot_fdsoi

# Set up paths
script_dir = pathlib.Path(__file__).parent.resolve()
path_mesh = script_dir / "meshes"
path_out = script_dir / "output"

# Load the perturbative value of the exchange coupling from a previous
# calculation
path_exchange_perturbation = path_out / "exchange_perturbation.txt"
exchange_perturbation = np.loadtxt(path_exchange_perturbation)

# Load the mesh
scaling = 1e-9
path_file = str(path_mesh / "refined_dqdfdsoi.msh2")
mesh = Mesh(scaling, path_file)

# Define the gate biases
back_gate_bias = -0.5
barrier_gate_1_bias = 0.5
plunger_gate_1_bias = 0.59
barrier_gate_2_bias = 0.57
plunger_gate_2_bias = 0.59
barrier_gate_3_bias = 0.5

# Define the device
dvc = get_double_dot_fdsoi(mesh, back_gate_bias, barrier_gate_1_bias, 
    plunger_gate_1_bias, barrier_gate_2_bias, plunger_gate_2_bias, 
    barrier_gate_3_bias)

# List of regions forming the double quantum dot region
dot_region_list = ["oxide_dot", "gate_oxide_dot", "buried_oxide_dot", 
                    "channel_dot"]

# Create the submesh object for the dot region
submesh = SubMesh(dvc.mesh, dot_region_list)

# Load the electric potential from the previous run to speed up Poisson 
# calculation
phi = io.load(path_out / "tunnel_coupling_final_potential.hdf5",
    var_name="phi")
dvc.set_potential(phi)

# Configure the non-linear Poisson solver
params_poisson = PoissonSolverParams()
params_poisson.tol = 1e-3 # Convergence threshold (tolerance) for the error

# Instantiate Poisson solver
poisson_slv = PoissonSolver(dvc, solver_params=params_poisson)

# Solve Poisson's equation
poisson_slv.solve(initialize=False)

# Configure the Schrödinger solver
params_schrod = SchrodingerSolverParams()
params_schrod.tol = 1e-6      # Tolerance on energies in eV

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

# Create a subdevice for the dot region
subdvc = SubDevice(dvc, submesh)

# Create a Schrodinger solver, and run it
print("Solving the single-electron Schrödinger equation")
schrod_solver = SchrodingerSolver(subdvc, solver_params=params_schrod)
schrod_solver.solve()
subdvc.print_energies()

# Solving the many-body hamiltonian
# Instantiate many-body solver
manybody_solver_params = ManyBodySolverParams()
manybody_solver_params.n_degen = 2
manybody_solver_params.num_states = 8
manybody_solver_params.num_particles = [2]
manybody_solver_params.overlap = True
manybody_solver = ManyBodySolver(subdvc, 
    solver_params=manybody_solver_params)

# solve the many-body hamiltonian
t0 = timer()
manybody_solver.solve()
print("Two-electron solution found in {} s".format(timer()-t0))

# Get the two-electron eigenenergies
two_electron_energies = subdvc.get_N_particle_subspace(2).eigval/ct.e
print("Two-electron energies")
print(two_electron_energies[0:8])
np.savetxt(path_out/"exchange_two_electron_energies.txt", two_electron_energies)

# Calculate the exchange interaction strength
exchange = two_electron_energies[1] - two_electron_energies[0]

# Print the output, and compare it with perturbation theory
print ('Exchange coupling (exact diagonalization): {:.3f} MHz'\
    .format(exchange*ct.e/ct.h/1e6))
print ('Exchange coupling (perturbation theory): {:.3f} MHz'\
    .format(exchange_perturbation/ct.h/1e6))