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.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).

../../../_images/point_charge_ground.png

Fig. 19.5.1 Ground-state wavefunction localized under plunger gate 1.

../../../_images/point_charge_Excited.png

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
../../../_images/pc_refinement.png

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).

../../../_images/DeltaE_oxide.png

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.

../../../_images/DeltaE_oxide_complete.png

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
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

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

# 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


# 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

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

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

    # Load the data
    data = np.loadtxt(E_file)
    x = data[:, 0]
    y = data[:, 1]
    gap = data[:, 3] - data[:, 2]

    # 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()