19. Point charges in a double quantum dot in FD-SOI
19.1. Requirements
19.1.1. Software components
QTCAD
Gmsh
19.1.2. Geometry file
qtcad/examples/tutorials/meshes/dqdfdsoi.geo
19.1.3. Python script
qtcad/examples/tutorials/fdsoi_point_charges.py
19.1.4. References
19.2. Briefing
In this tutorial, we will investigate the effect of a point charge on the orbital levels of a quantum dot in the fully-depleted silicon-on-insulator (FD-SOI) geometry introduced in A double quantum dot device in a fully-depleted silicon-on-insulator transistor. Specifically, we will investigate how the spectrum of a quantum dot changes as a function of the position of a point charge in the device. In QTCAD, point charges are handled in accordance with the theoretical formulation presented in Point Charges.
19.3. Geometry file
As before, the initial mesh is generated from a
geometry file qtcad/examples/tutorials/meshes/dqdfdsoi.geo
by running Gmsh with
gmsh dqdfdsoi.geo
As in Tunnel coupling in a double quantum dot in FD-SOI—Part 1: Plunger gate tuning, this will produce a .msh file and a
.xao file. The .xao file format is appropriate for
geometries produced using the OpenCASCADE geometry kernel.
19.4. Header
We start by importing the necessary packages, classes, and functions.
import numpy as np
import matplotlib.pyplot as plt
import pathlib
import os
from qtcad.device.mesh3d import Mesh, SubMesh
from qtcad.device.device import Device, SubDevice
from qtcad.device.poisson import Solver, SolverParams
from qtcad.device.schrodinger import Solver as SchrodingerSolver
from qtcad.device.schrodinger import SolverParams as SchrodingerSolverParams
from qtcad.device import constants as ct
from examples.tutorials.double_dot_fdsoi import get_double_dot_fdsoi
All of these packages, classes, and functions have been introduced in previous tutorials.
19.5. Setting up the device
We start with setting up paths to various files, both inputs and outputs. The inputs include the mesh file, the geometry file, and the refined mesh file. The outputs are: a file that will contain the energy levels of the double-dot system as a function of the point charge position and a file containing a plot of the orbital level spacing as a function of the point-charge position.
# Files -----------------------------------------------------------------------
script_dir = pathlib.Path(__file__).parent.resolve()
path_mesh = script_dir / "meshes"
mesh_file = str(path_mesh / "dqdfdsoi.msh")
geo_file = str(path_mesh / "dqdfdsoi.xao")
ref_file = str(path_mesh / "refined_dqdfdsoi_oxide.msh")
path_out = script_dir / "output"
E_file = str(path_out / "fdsoi_energies_oxide_charge.txt")
Delta_file = str(path_out / "DeltaE_oxide.png")
We then setup the calculation by loading the mesh and defining the gate bias parameters.
# Setup -----------------------------------------------------------------------
# Load the mesh
scaling = 1e-9
mesh = Mesh(scaling, mesh_file)
# Define the gate bias parameters
back_gate_bias = -0.5
barrier_gate_1_bias = 0.5
plunger_gate_1_bias = 0.7
barrier_gate_2_bias = 0.5
plunger_gate_2_bias = 0.6
barrier_gate_3_bias = 0.5
Here we have detuned the plunger-gate voltages such that the ground and first excited states of the system are localized under plunger gate 1 (see Fig. 19.5.1 and Fig. 19.5.2).
Fig. 19.5.1 Ground-state wavefunction localized under plunger gate 1.
Fig. 19.5.2 Excited-state wavefunction localized under plunger gate 1.
19.6. Basic electrostatics
Next, we create a function, func_poisson that can position point charges
throughout a double-dot FD-SOI device generated via the get_double_dot_fdsoi
function and apply the adaptive non-linear Poisson solver.
The inputs of func_poisson are the positions of the point charges, the
charges of each of the point charges, and the smearing radius of the point
charges.
It returns the device over which the non-linear Poisson solver has been
applied.
The difference between this tutorial and previous ones is the use of the
add_point_charges method.
The arguments of this method are the positions of the point charges, the
charges of the point charges, the smearing radius of the point charges, and a
boolean that indicates whether the mesh should be refined around the point
charges (see Fig. 19.6.1).
This refinement is often necessary since the smearing radius is typically taken
to be smaller than the characteristic length of the mesh
elements
(\(\sigma \lesssim 0.1 \, \mathrm{nm} < 1 \, \mathrm{nm} \lesssim h\),
where \(\sigma\) is the smearing radius and \(h\) is a typical characteristic
length of a QTCAD mesh), which can lead to situations where the point charges
are not observed by the solvers.
# Poisson ---------------------------------------------------------------------
def func_poisson(pos: np.ndarray, Q: np.ndarray, r: float) -> Device:
"""Apply the adaptive non-linear Poisson solver with point charges.
Args:
pos (np.ndarray): Position of the point charges
Q (np.ndarray): Charge of the point charges
r (float): Smearing radius of the point charges
Returns:
device: Device over which the non-linear Poisson solver has been
applied (point charges have been included).
"""
# Define the device object from the function defined in the FD-SOI tutorial
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)
dvc.add_point_charges(pos, Q, r=r, refine=True)
# Configure the Non-Linear Poisson solver
solver_params = SolverParams()
solver_params.tol = 1e-3
solver_params.initial_ref_factor = 0.1
solver_params.final_ref_factor = 0.75
solver_params.min_nodes = 50000
solver_params.max_nodes = 1e5
solver_params.maxiter_adapt = 30
solver_params.maxiter = 200
solver_params.dot_region = [
"oxide_dot", "gate_oxide_dot",
"buried_oxide_dot", "channel_dot"
]
solver_params.h_dot = 0.8
solver_params.refined_mesh_filename = ref_file
# Solve Poisson
slv = Solver(dvc, solver_params=solver_params, geo_file=geo_file)
slv.solve()
return dvc
Fig. 19.6.1 Mesh produced by the adaptive meshing algorithm after the addition of a point charge at the interface between the silicon channel and the oxide separating it from the front gates, at the center of plunger gate 1, \((x, y) = (0, -17.5 \, \mathrm{nm})\). Besides the mesh being adapted to more easily solve the non-linear Poisson equation, the mesh is also refined around the point charge to ensure that it is observed by the solvers.
19.7. Schrödinger solver
We also write the function func_schrodinger that applies the Schrödinger
solver to a subdevice of the full FD-SOI device.
This subdevice is defined by the regions that form the double quantum dot.
The function takes a device object as input (i.e. the FD-SOI device generated
by func_poisson) and returns the subdevice object over which the
Schrödinger equation has been solved.
def func_schrodinger(d: Device) -> SubDevice:
""" Apply the Schrodinger solver to the double-dot region of the FD-SOI
device.
Args:
d (Device): An FD-SOI device over which the non-linear Poisson
solver has been applied.
Returns:
SubDevice: The subdevice object over which the Schödinger equation
has been solved.
"""
d.set_V_from_phi()
# List of regions forming the double quantum dot region
dot_region_list = [
"oxide_dot", "gate_oxide_dot", "buried_oxide_dot", "channel_dot"
]
# Create a submesh including only the dot region
submesh = SubMesh(d.mesh, dot_region_list)
# Create a subdevice object
subdevice = SubDevice(d, submesh)
# Configure the Schrodinger solver
solver_params = SchrodingerSolverParams()
solver_params.num_states = 10
solver_params.tol = 1e-6
# Solve Schrodinger
slv = SchrodingerSolver(subdevice, solver_params=solver_params)
slv.solve()
return subdevice
19.8. Generate Data
Now that we have created the two functions described above, we solve the Poisson and Schrödinger equations while varying the position of a single point charge to determine the effect of the point-charge position on the quantum-dot spectrum. In this example, we place a single point charge, \(Q = e\), at \(z=0\), the interface between the undoped silicon channel and the oxide that separates it from the front gates and vary its in-plane position using a for loop for each in-plane coordinate.
# Calculate the data ----------------------------------------------------------
if __name__ == '__main__':
# Loop over point-charge positions
npts_dir = 3 # Number of points per direction
x0, y0 = 0, -17.5 # Center of plunger gate 1
xmin, ymin = -20, -25 # Minimum position
x_var = np.linspace(xmin, x0, npts_dir)
x_const = x0 * np.ones(npts_dir)
y_var = np.linspace(ymin, y0, npts_dir)
y_const = y0 * np.ones(npts_dir)
x_values = np.append(x_var[0:npts_dir-1], x_const)
y_values = np.append(y_const[0:npts_dir-1], y_var)
for i, x in enumerate(x_values):
y = y_values[i]
pos = np.array([[x, y, 0]]) * scaling # position of the point charge
Q = np.array([1]) * ct.e # charge of the point charge
r = 1e-10 # radius of the point charge
dvc = func_poisson(pos, Q, r) # Solve Poisson
subdevice = func_schrodinger(dvc) # Solve Schrodinger
subdevice.print_energies() # Print energies
# Save energies
out = np.append(np.array([x, y]), subdevice.energies / ct.e)
# header
E_header = '' # set empty header
if not os.path.isfile(E_file): # checks if the file exists
# if it doesn't then add the header
E_header = "x (nm), y (nm), energies (eV)"
with open(E_file, 'a') as file:
np.savetxt(file, np.array(out)[np.newaxis, :], header = E_header)
For every calculation, we save the results in a text file, E_file, for
potential future use.
At the end of this loop, the file will contain the results of 5 calculations,
each corresponding to a different point-charge position.
This is expected to take ~20 minutes to run on a modern laptop.
19.9. Visualizing the Results
To visualize the results, we start by loading the data from the text file we
saved.
We also create the gap variable which contains the energy gap between the
ground and first excited states for every point-charge position.
# Load the data
data = np.loadtxt(E_file)
x = data[:, 0]
y = data[:, 1]
gap = data[:, 3] - data[:, 2]
We then plot the energy gap between the ground and first excited states as a function of the point-charge position.
# Plot labels
label = '$\Delta E$ (eV)'
title = 'Orbital level spacing'
xlabel = '$x$ (nm)'
ylabel = '$y$ (nm)'
# Create figure and axes
fig, ax = plt.subplots(figsize=(8, 6))
# Plot the data
scatter = ax.scatter(x, y, c=gap, cmap='viridis')
fig.colorbar(scatter, ax=ax, label=label)
# Set labels
ax.set_title(title)
ax.set_xlabel(xlabel)
ax.set_ylabel(ylabel)
ax.grid(True)
# Display the plot
plt.savefig(str(Delta_file))
plt.show()
The above code snipet generates a plot where the colorbar corresponds to the energy gap between the ground and first excited states and the x and y axes correspond to the in-plane position of the point charge (see Fig. 19.9.1).
Fig. 19.9.1 Orbital level spacing as a function of the point-charge position for the quantum dot located below plunger gate 1. The computational time for this example is ~20 minutes on a modern laptop.
The center of plunger gate 1 is located at \((x, y) = (0, -17.5 \, \mathrm{nm})\). The peak of the ground-state wavefunction and the node of the excited-state wavefunction coincides with this point (see Fig. 19.5.1 and Fig. 19.5.2). Accordingly, the energy gap between the ground and first excited states is largest when a point charge is placed at this position and is a factor of 2 larger than the gap when the point charge is far from the center.
Note
A more complete idea of the dependence of the orbital level spacing of the dot situated beneath plunger gate 1 as a function of the point charge position can be obtained by computing the spectrum at more points. For a more complete example (which took substantially more time to compute), see Fig. 19.9.2.
Fig. 19.9.2 Orbital level spacing as a function of the point-charge position for the quantum dot located below plunger gate 1. More data is presented here than in Fig. 19.9.1. This plot took substantially more time to compute.
19.10. Full code
__copyright__ = "Copyright 2024, Nanoacademic Technologies Inc."
import numpy as np
from pathlib import Path
from matplotlib import pyplot as plt
from scipy.optimize import minimize_scalar
from timeit import default_timer as timer
from qtcad.device import constants as ct
from qtcad.device import io
from qtcad.device import analysis as an
from qtcad.device.mesh3d import Mesh, SubMesh
from qtcad.device import SubDevice
from qtcad.device.poisson import Solver as PoissonSolver
from qtcad.device.poisson import SolverParams as PoissonSolverParams
from qtcad.device.schrodinger import Solver as SchrodingerSolver
from qtcad.device.schrodinger import SolverParams as SchrodingerSolverParams
from examples.tutorials.double_dot_fdsoi import get_double_dot_fdsoi
# Set up paths
script_dir = Path(__file__).parent.resolve()
path_mesh = script_dir / "meshes"
path_geo = path_mesh / "dqdfdsoi.xao"
path_out = script_dir / "output"
# Load the mesh
scaling = 1e-9
file = path_mesh / "dqdfdsoi.msh"
mesh = Mesh(scaling, file)
# Define gate biases
back_gate_bias = -0.5
barrier_gate_1_bias = 0.5
plunger_gate_1_bias = 0.59
barrier_gate_2_bias_low = 0.51
barrier_gate_2_bias_high = 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_high, 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"]
# Configure the non-linear Poisson solver
params_poisson = PoissonSolverParams()
params_poisson.tol = 1e-3 # Convergence threshold (tolerance) for the error
params_poisson.initial_ref_factor = 0.1
params_poisson.final_ref_factor = 0.75
params_poisson.min_nodes = 50000
params_poisson.max_nodes = 1e5
params_poisson.maxiter_adapt = 30
params_poisson.maxiter = 200
params_poisson.dot_region = dot_region_list
params_poisson.h_dot = 0.8
params_poisson.refined_mesh_filename = path_mesh / "refined_dqdfdsoi.msh"
# Instantiate Poisson solver
poisson_slv = PoissonSolver(dvc, solver_params=params_poisson,
geo_file=path_geo)
# Solve in the original gate bias configuration
poisson_slv.solve()
# Produce a linecut of the conduction band edge along the channel
x, y, z = dvc.mesh.glob_nodes.T
ymin = np.min(y); ymax = np.max(y)
distance, linecut_high_bias = an.linecut(dvc.mesh, dvc.cond_band_edge(),
(0, ymin, -1e-9), (0, ymax, -1e-9))
# Solve Poisson again at a lower central barrier gate bias
dvc.set_applied_potential("barrier_gate_2_bnd", barrier_gate_2_bias_low)
poisson_slv.solve(initialize=False)
# Produce another linecut of the conduction band edge along the channel
distance, linecut_low_bias = an.linecut(dvc.mesh, dvc.cond_band_edge(),
(0, ymin, -1e-9), (0, ymax, -1e-9))
# Plot the linecuts
fig = plt.figure(figsize=(8,5))
ax = fig.add_subplot(1,1,1)
ax.set_xlabel("Distance along linecut (nm)")
ax.set_ylabel("Conduction band edge, $E_C$ (eV)")
ax.plot(distance/1e-9, linecut_high_bias/ct.e, "-r",
label=f"Barrier gate 2 bias: {barrier_gate_2_bias_high} V")
ax.plot(distance/1e-9, linecut_low_bias/ct.e, "-b",
label=f"Barrier gate 2 bias: {barrier_gate_2_bias_low} V")
ax.legend()
plt.show()
# Create a submesh including only the dot region
submesh = SubMesh(dvc.mesh, dot_region_list)
# Instantiate Schrodinger solver parameters
params_schrod = SchrodingerSolverParams()
params_schrod.tol = 1e-6 # Tolerance on energies in eV
# defining the function for optimization.
def function(x, gate_bias, path_states = None, path_phi = None):
""" Tunnel splitting to be optimized.
Args:
x (float): Gate bias detuning with respect to reference configuration.
gate_bias (float): Gate bias in the reference configuration.
path_states (str): Path in which to save the ground and first
excited state wavefunctions in .vtu format. If None, do not save.
path_phi (str): Path in which to save the electric potential in .hdf5
format.
Returns:
float: The tunnel splitting (in eV).
"""
print("="*80)
print("Evaluating double-dot splitting for detuning: {} V".format(x))
print("="*80)
# Set the plunger gate bias
dvc.set_applied_potential("plunger_gate_2_bnd", gate_bias + x)
# Solve the non-linear Poisson equation
poisson_slv.solve(initialize=False)
# 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)
# Instantiate and run a Schrödinger solver
schrodinger_slv = SchrodingerSolver(subdvc, solver_params = params_schrod)
schrodinger_slv.solve()
# Store the solution as an initial guess for the next run
params_schrod.guess = (subdvc.eigenfunctions, subdvc.energies)
energies = subdvc.energies
if path_states is not None:
out_dict = {
"Ground state": subdvc.eigenfunctions[:,0],
"First excited state": subdvc.eigenfunctions[:,1]
}
io.save(path_states, out_dict, submesh)
if path_phi is not None:
out_dict = {"phi": dvc.phi}
io.save(path_phi, out_dict)
return (energies [1] - energies [0])/ct.e
# Tune the bias of plunger gate 2 to find a symmetric double dot
# configuration
t0 = timer()
result = minimize_scalar(function, bounds=(-200e-6,200e-6), args=(
np.array([plunger_gate_2_bias])), method = 'bounded',
options={"xatol": 1e-7})
# print results:
print (f'Optimum found at {result.x} V with the value: {result.fun} eV.')
print(f'Computation time: {timer()-t0} s')
# Save the optimum and detuning in a file
path_result = path_out / "detuning_and_coupling.txt"
tab = np.array([result.x, result.fun])
np.savetxt(path_result, tab)
# Save the ground and first excited state wave functions in .vtu format
# Save the electric potential in .hdf5 format
print("Solving Poisson and Schrödinger for the final detuning value.")
print("Saving eigenfunctions in .vtu format.")
path_states = str(path_out / "tunnel_coupling_wavefunctions_high_barrier.vtu")
path_phi = str(path_out / "tunnel_coupling_final_potential.hdf5")
function(result.x, np.array(plunger_gate_2_bias), path_states=path_states,
path_phi=path_phi)
# Save the ground and first excited state wave functions in .vtu format
# at high central barrier gate bias (strong tunneling)
dvc.set_applied_potential("barrier_gate_2_bnd", barrier_gate_2_bias_high)
path_states = str(path_out / "tunnel_coupling_wavefunctions_low_barrier.vtu")
function(result.x, np.array(plunger_gate_2_bias), path_states=path_states)