21. Periodic boundary condition in a quantum dot in FD-SOI

21.1. Requirements

21.1.1. Software components

  • QTCAD

  • Gmsh

21.1.2. Geometry file

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

21.1.3. Python script

  • qtcad/examples/tutorials/periodic.py

21.1.4. References

21.2. Briefing

In this tutorial, we will demonstrate how to use QTCAD to solve the Poisson equation with a periodic boundary condition. The example we will consider is a quantum dot generated using three clavier gates in a geometry with similar features to the FD-SOI structure mentioned in A double quantum dot device in a fully-depleted silicon-on-insulator transistor. By enforcing periodic boundary conditions along the channel direction and analyzing time-varying gate voltage configurations, we simulate an infinite linear array of dots and estimate (conveyer-belt) shuttling [BLP+23, JBF25, NTWM25] properties such as the position and speed of the electrons.

To simulate the time-dependent potential in the device, we use a method described in previous works [JBF25]. In this approach, the potential is first calculated with only one clavier gate active at a time. Then a linear combination of these solutions is used to approximate the full time-dependent potential accros the device.

21.3. Geometry File

In this example, we use a geometry with similar structure to the one presented in A double quantum dot device in a fully-depleted silicon-on-insulator transistor, but with three gates as shown in the following figure.

FDSOI structure

Fig. 21.3.1 Fully-depleted silicon-on-insulator (FD-SOI) structure: (a) Three- dimensional view of the structure, showing the height of the device and the channel width (5 nm). The periodic boundary condition is applied to the edge surfaces along the y-direction, labeled as left_bnd and right_bnd. Additionally, the bottom gate of the device is shown and labeled as back_gate_bnd. (b) Top view of the FD-SOI device showing the clavier gates labeled as A0, A1, and A2, along with their dimensions.

To enforce periodic boundary conditions in QTCAD, a periodic mesh is required. A periodic mesh has identical node distributions on the surfaces over which the periodic boundary condition is applied. To achieve this in a 3D geometry, we use Gmsh and specify the periodic condition for a pair of surfaces using the Periodic Surface keyword in the geometry file, as shown below:

Periodic Surface{right_bnd[]} = {left_bnd[]} Translate {0, channel_len, 0};

Here, right_bnd[] and left_bnd[] represent the surface tags of the two periodic boundaries. The Translate {0, channel_len, 0} command defines the vector connecting them, where channel_len represents the channel length.

Once this line is added, we generate the initial mesh by running:

gmsh periodic_fdsoi.geo

This will create a mesh file in .msh format and a geometry file in .xao format.

In this tutorial, we will only use the .msh file since we will be working with a static mesh to accelerate the simulation. However, adaptive meshing is also possible using an .xao raw geometry file when using periodic boundary conditions. To learn how to use adaptive meshing, please refer to the Poisson solver with adaptive meshing and Adaptive-mesh Poisson solver combined with the Schrödinger solver tutorials.

Note

For adaptive meshing with periodic boundary conditions, only the .xao format is supported. Other geometry file formats are not supported.

21.5. Setting up the device

In this step, we define paths and load the mesh.

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

# Load the mesh
scaling = 1e-9
file = path_mesh / "periodic_fdsoi.msh"
mesh = Mesh(scaling, file)

To simulate conveyor-belt shuttling, we apply a phase-modulated sinusoidal voltage to all gates arranged along the shuttling direction. Concretely, the potential applied to the gate \(A_i\) is given by:

(21.5.1)\[V_{A_i}(t) = V_0 \cos\left(\Omega t - \frac{2\pi i}{N}\right)\]

where \(V_0\) is the amplitude of the voltage, \(\Omega = 2\pi / \tau\) is the angular frequency of the time-dependent oscillation, \(i\) is the index of the clavier gate, and \(N\) is the total number of gates.

# Define gate biases
back_gate_bias = -0.5
N = 3 # Number of gates
V_0 = 1e-1 # 100 mV amplitude
tau = 1e-9 # 1 ns period
Omega = 2*np.pi/tau
# Define the gate biases for the first clavier gate simulation.
V_A0 = 1
V_A1 = 0
V_A2 = 0

# define the gate bias function:
def V_A_i(t, i):
    """
    Calculates the time-dependent gate bias of the clavier gate i.

    Args:
        t (float): Time in seconds.
        i (int): Index of the clavier gate.
    Returns:
        float: Time-dependent gate bias.

    """
    return V_0*np.cos(Omega*t - ((2*np.pi*i)/N))

Similar to previous tutorials, we use a Device object and define the regions and the boundaries of the device. We apply the periodic boundary condition on the left and right boundaries of the device, by using the new_periodic_bnd method. This method takes two boundary spcifiers as arguments, which in this tutorial are left_bnd and right_bnd.

# Define the device object
dvc = Device(mesh, conf_carriers='e')
dvc.set_temperature(0.1)

# Create the regions
dvc.new_region("oxide", mt.SiO2)
dvc.new_region("gate_oxide", mt.HfO2)
dvc.new_region("buried_oxide", mt.SiO2)
dvc.new_region("channel", mt.Si)

# Set up boundary conditions
Ew = mt.Si.Eg/2 + mt.Si.chi # Midgap
dvc.new_gate_bnd("A0", V_A0, Ew)
dvc.new_gate_bnd("A1", V_A1, Ew)
dvc.new_gate_bnd("A2", V_A2, Ew)
dvc.new_frozen_bnd("back_gate_bnd", back_gate_bias, mt.Si, 1e15*1e6,
    "n", 46*1e-3*ct.e)
# Periodic boundary conditions:
dvc.new_periodic_bnd("left_bnd", "right_bnd")

21.6. Basic electrostatics

To simulate the electrostatics of the device as a function of time, we solve the electrostatic potential when only one of the clavier gates is turned on (\(V_{A_i}=1\)) and the potential of the other gates is set to zero. After repeating this process for all the gates, we can find the time-dependent potential using the following equation [JBF25].

(21.6.1)\[\phi(x,y,z,t) = \sum_{i=0}^{N-1} V_{A_i}(t) \cdot \phi_i(x,y,z)\]

where \(\phi_i\) is the potential of the device when only the gate \(A_i\) is turned on and \(V_{A_i}(t)\) is the time-dependent voltage of the gate \(A_i\) given by Eq. (21.5.1).

The potential is simulated using the linear Poisson solver with a static mesh.

Note

The linear and non-linear Poisson solvers both support periodic boundary conditions. This holds for both static and adaptive meshing.

We do this by defining a SolverParams object for the linear Poisson Solver, by constructing the Solver itself, and by running it in the next section using the solve method.

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

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

We use the set_applied_potential method to set the potential for each clavier gate activation and then solve the electrostatics of the device to find the confinement potential. In this tutorial, we save every activation potential along a linecut in the channel direction (\(y\)-direction), \(1\;\mathrm{nm}\) below the interface between the silicon channel and the gate oxide.

# Clavier gate biases for each activation
activate_bias = [(1,0,0), (0,1,0), (0,0,1)]
# Define the linecut
x, y, z = dvc.mesh.glob_nodes.T
ymin = np.min(y);   ymax = np.max(y)
height = -1e-9  # 1 nm below the silicon channel and the gate oxid

# initialize a list for saving the linecuts
potential_linecuts = []

for idn, bias in enumerate(activate_bias):
    print(f'Calculating potential for activation {idn+1} of {N}')
    # Update gate biases
    dvc.set_applied_potential("A0", bias[0])
    dvc.set_applied_potential("A1", bias[1])
    dvc.set_applied_potential("A2", bias[2])

    # Solve
    poisson_slv.solve()

    # Calculate the confinement potential
    dvc.set_V_from_phi()

    # Save potential and confinement potential
    distance, linecut_pot = an.linecut(dvc.mesh, dvc.phi,
                    (0, ymin, -1e-9), (0, ymax, -1e-9), method="pyvista")
    distance, linecut_conf_pot = an.linecut(dvc.mesh, dvc.V/ct.e,
                    (0, ymin, -1e-9), (0, ymax, -1e-9), method="pyvista")
    potential_linecuts.append([distance, linecut_pot, linecut_conf_pot])

# save the linecuts
output_file = path_out / "potential_linecuts.npy"
np.save(output_file, np.array(potential_linecuts))

We note that an overlap between the essential boundaries (e.g., new_gate_bnd, new_frozen_bnd, …) and the periodic boundaries (e.g., new_periodic_bnd) leads to a conflict between boundary conditions that we resolve by removing the essential boundary condition from the overlapping nodes. This situation typically occurs at the edges or corners of the device where the two boundary types meet.

In this tutorial, such an overlap happens at the bottom edges of the periodic boundaries which are in contact with the back gate. Therefore, the solver will remove the overlapping nodes (bottom edges) from the frozen boundary. This process is shown as a warning message in the terminal:

UserWarning: Periodic boundary nodes overlap with essential boundary nodes.
    Essential boundary conditions are not enforced at these nodes.

We then construct the time-dependent potential using Eq. (21.6.1) along the specified linecut. Similarly, since the confinement potential \(V(x,y,z,t)\) depends linearly on the electrostatic potential \(\phi(x,y,z,t)\), we calculate it using the following equation:

(21.6.2)\[V(x,y,z,t) = \sum_{i=0}^{N-1} V_{A_i}(t) \cdot V_i(x,y,z)\]

where \(V_i\) is the confinement potential of the device when only the gate \(A_i\) is activated.

# Load the linecuts
potential_linecuts = np.load(output_file)

# Finding the time-dependent potential along the linecut
def confinement_potential(t):
    """
    Calculates the potential along the linecut.

    Args:
        t (float): Time in seconds.
    Returns:
        distance (1d array): Distance along the linecut.
        potential (2d array): Electrostatic potential linecut for the given t
        confinement_pot (1d array): Confinement potential linecut for the given t
    """
    # Calculate the time-dependent part
    V_A0 = V_A_i(t, 0)
    V_A1 = V_A_i(t, 1)
    V_A2 = V_A_i(t, 2)

    # initialize the confinement potential
    potential = 0
    confinement_pot = 0
    for i in range(N):
        potential += V_A_i(t, i)*potential_linecuts[i, 1, :]
        confinement_pot += V_A_i(t, i)*potential_linecuts[i, 2, :]

    return distance, potential, confinement_pot

We plot the electrostatic potential and confinement potential linecut for \(t=\frac{\tau}{3}\) using the following script:

# calculate the confinement potential at t=tau/3
dist, phi_t, V_t = confinement_potential(tau/3)
# Plot the potential linecut
fig, ax = plt.subplots(figsize=(8,5))
ax.set_xlabel("Distance along linecut [nm]")
ax.set_ylabel("Electrostatic potential ($\\phi$) [mV]")
ax.plot(dist/1e-9, phi_t/1e-3, "-r")
ax.axhline(phi_t[0]/1e-3, color="black", lw=0.5, ls="--", alpha=0.5)

# plot the confinement potential
ax2 = ax.twinx()
ax2.set_ylabel("Confinement-potential ($V$) [meV]", color="blue")
ax2.plot(dist/1e-9, V_t/1e-3, "-b", label="Confinement potential")
ax2.tick_params(axis="y", labelcolor="blue")
ax2.axhline(V_t[0]/1e-3, color="black", lw=0.5, ls="--", alpha=0.5)

plt.show()

This gives:

Periodic potential along the channel

Fig. 21.6.1 Periodic potential and confinement potential along the channel at \(t= \frac{\tau}{3}\). As shown in this figure, the values at the both ends of the channel is the same, indicating the periodicity of the potential.

After demonstrating the periodicity of the potential, we can now proceed to simulate a simple shuttling example by varying the gate voltages with respect to time.

We take time_steps=10 time steps, and plot the confinement potential linecut from \(t=\frac{\tau}{3}\) to \(t=\frac{2\tau}{3}\) in steps of \(\tau/10\).

# Varying the gate voltages with respect to time
time_steps = 10
time = np.linspace(tau/3, 2*tau/3, time_steps)

# plot the confinement potential linecuts
fig, ax = plt.subplots(figsize=(8, 5))
ax.set_xlabel("Distance along linecut [nm]")
ax.set_ylabel("confinement potential ($V$) [meV]")
for idt, t in enumerate(time):
    dist, phi_t, V_t = confinement_potential(t)
    ax.plot(dist/1e-9, V_t/1e-3, label=f"t = {t * 1e9:.1f} ns")
ax.annotate("",(13,30),(12,30), arrowprops=dict(facecolor='black', arrowstyle='-|>'))
ax.text(12.5, 40, "Time Progression", fontsize=10, ha='center', va='center')
ax.legend()
plt.show()

The above scripts produce the following plot:

Periodic potential time steps along the channel

Fig. 21.6.2 Periodic confinement potential along the channel in different time intervals.

Following a previous work on conveyor-belt spin shuttling in silicon [NTWM25], we can approximate the confinement potential using the following quadratic function along the linecut

(21.6.3)\[V(x) = \frac{1}{2} m^* \omega^2 (x - x_{\text{dot}})^2,\]

where \(m^*= 0.19 m_e\) is the transverse effective mass of electrons in silicon, \(\omega\) is the angular frequency of the harmonic oscillator and \(x_\mathrm{dot}\) is the position of the quantum dot.

Under this assumption, the position of the quantum dot is approximated as the minimum of the confinement potential. Additionally, the size of the quantum dot is approximated by the following equation [NTWM25].

(21.6.4)\[a = \sqrt{\frac{\hbar}{m^*\omega}},\]

We use the np.polyfit function to fit a second-order polynomial to the confinement potential linecut.

We note that the approximation of the potential as a second-order polynomial is only valid near the minima of the confinement potential. Therefore, we need to select a range around the minimum of the potential for fitting. We can do this by finding the minimum of the potential and selecting a range (e.g., 5 nm) around it. The following script shows this process:

# Fit a second-order polynomial near the potential's minimum.
fit_range = 5e-9
# Lists to store the fitting parameters.
coef_list = []
x_min_list = []
fig, ax = plt.subplots(figsize=(8, 5))
ax.set_xlabel("Distance along linecut [nm]")
ax.set_ylabel("confinement potential ($V$) [meV]")
for idt, t in enumerate(time):
    dist, _, V_t = confinement_potential(t)
    # Find the minimum of the linecut.
    min_idx = np.argmin(V_t)
    min_dist = dist[min_idx]

    # Select the range around the minimum for fitting
    fit_mask = (dist >= min_dist - fit_range) & (dist <= min_dist + fit_range)
    fit_dist = dist[fit_mask]
    fit_potential = V_t[fit_mask]

    # Fit a second-order polynomial
    coeffs = np.polyfit(fit_dist, fit_potential, 2)
    coef_list.append(coeffs)

    # create a polynomial function from the coefficients
    poly = np.poly1d(coeffs)

    # Calculate the approximated polynomial's minimum
    x_min = -coeffs[1]/(2*coeffs[0])
    x_min_list.append(x_min)    # save the minimum of the polynomial
    y_min = poly(x_min)

    # Plot the results
    line = ax.plot(dist/1e-9, V_t/1e-3, label=f"t = {t * 1e9:.1f} ns")
    ax.plot(fit_dist/1e-9, poly(fit_dist)/1e-3,'--', color=line[0].get_color())
    ax.plot(x_min/1e-9, y_min/1e-3, "ko")
ax.annotate("", (13, 30), (12, 30), arrowprops=dict(facecolor='black', arrowstyle='-|>'))
ax.text(12.5, 40, "Time Progression", fontsize=10, ha='center', va='center')
ax.legend()
plt.show()

The output of the above script is shown below, where the dashed lines represent the fitted polynomial approximation of the confinement potential. The black dots represent the approximate position of the quantum dot at each time step.

Periodic potential along the channel with polynomial approximation

Fig. 21.6.3 Confinement potential along the linecut with polynomial approximation at various time intervals shown with dashed lines.

Now we can extract the approximate position, speed, and size of the quantum dot from the fitted polynomial at each time step. We emphasize that this is a simple approximation where we have condensed the full three dimensional problem to an analysis of a one dimensional linecut of the confinement potential to highlight the features of a periodic potential in a conveyer-belt shuttling simulation.

# Extract the position, speed and size of the quantum dot.
coef_list = np.array(coef_list)
x_min_list = np.array(x_min_list)
# Shift the position to start from 0
x_min_list = x_min_list - x_min_list[0]
w = np.sqrt(2*ct.e*coef_list[:, 0]/(0.19*ct.me))
size = np.sqrt(ct.hbar/(0.19*ct.me*w)) # Size of the quantum dot
speed = np.gradient(x_min_list, time)

# Plot the position of the quantum dot over time
fig, ax = plt.subplots(figsize=(8, 5))
ax.set_xlabel("Time [ns]")
ax.set_ylabel("Position [nm]")
ax.fill_between(time*1e9, (x_min_list - size/2)/1e-9, (x_min_list + size/2)/1e-9,
                color='gray', alpha=0.3, label="Quantum Dot Size")
ax.plot(time*1e9, x_min_list/1e-9, 'o-', label="Position")
ax.legend()
plt.show()

# Plot the speed of the quantum dot over time
fig, ax = plt.subplots(figsize=(8, 5))
ax.set_xlabel("Time [ns]")
ax.set_ylabel("Speed [m/s]")
ax.plot(time*1e9, speed, 'o-')
plt.show()

# Plot the Size of the quantum dot over time
fig, ax = plt.subplots(figsize=(8, 5))
ax.set_xlabel("Time [ns]")
ax.set_ylabel("Size [nm]")
ax.plot(time* 1e9, size/1e-9, 'o-')
plt.show()

The results of the above script are shown below:

The approximated position of the quantum dot

Fig. 21.6.4 The approximated position of the quantum dot with the dot size highlighted as a shaded area around each point.

The approximated speed of the quantum dot

Fig. 21.6.5 The approximated speed of the quantum dot. The speed is calculated using the derivative of the position with respect to time.

The approximated size of the quantum dot

Fig. 21.6.6 The approximated size of the quantum dot using Eq. (21.6.4).

21.7. Full code

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

import numpy as np
from pathlib import Path
from matplotlib import pyplot as plt
from qtcad.device import constants as ct
from qtcad.device import materials as mt
from qtcad.device import analysis as an
from qtcad.device.mesh3d import Mesh
from qtcad.device import Device
from qtcad.device.poisson_linear import Solver as PoissonSolver
from qtcad.device.poisson_linear import SolverParams as PoissonSolverParams

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

# Load the mesh
scaling = 1e-9
file = path_mesh / "periodic_fdsoi.msh"
mesh = Mesh(scaling, file)

# Define gate biases
back_gate_bias = -0.5
N = 3 # Number of gates
V_0 = 1e-1 # 100 mV amplitude
tau = 1e-9 # 1 ns period
Omega = 2*np.pi/tau 
# Define the gate biases for the first clavier gate simulation.
V_A0 = 1
V_A1 = 0
V_A2 = 0

# define the gate bias function:
def V_A_i(t, i):
    """
    Calculates the time-dependent gate bias of the clavier gate i.

    Args:
        t (float): Time in seconds.
        i (int): Index of the clavier gate.
    Returns:
        float: Time-dependent gate bias.

    """
    return V_0*np.cos(Omega*t - ((2*np.pi*i)/N))

# Define the device object
dvc = Device(mesh, conf_carriers='e')
dvc.set_temperature(0.1)

# Create the regions
dvc.new_region("oxide", mt.SiO2)
dvc.new_region("gate_oxide", mt.HfO2)
dvc.new_region("buried_oxide", mt.SiO2)
dvc.new_region("channel", mt.Si)

# Set up boundary conditions
Ew = mt.Si.Eg/2 + mt.Si.chi # Midgap
dvc.new_gate_bnd("A0", V_A0, Ew)
dvc.new_gate_bnd("A1", V_A1, Ew)
dvc.new_gate_bnd("A2", V_A2, Ew)
dvc.new_frozen_bnd("back_gate_bnd", back_gate_bias, mt.Si, 1e15*1e6, 
    "n", 46*1e-3*ct.e)
# Periodic boundary conditions:
dvc.new_periodic_bnd("left_bnd", "right_bnd")

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

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

# Clavier gate biases for each activation
activate_bias = [(1,0,0), (0,1,0), (0,0,1)]
# Define the linecut
x, y, z = dvc.mesh.glob_nodes.T
ymin = np.min(y);   ymax = np.max(y)
height = -1e-9  # 1 nm below the silicon channel and the gate oxid

# initialize a list for saving the linecuts
potential_linecuts = []

for idn, bias in enumerate(activate_bias):
    print(f'Calculating potential for activation {idn+1} of {N}')
    # Update gate biases
    dvc.set_applied_potential("A0", bias[0])
    dvc.set_applied_potential("A1", bias[1])
    dvc.set_applied_potential("A2", bias[2])

    # Solve
    poisson_slv.solve()

    # Calculate the confinement potential
    dvc.set_V_from_phi()

    # Save potential and confinement potential
    distance, linecut_pot = an.linecut(dvc.mesh, dvc.phi,
                    (0, ymin, -1e-9), (0, ymax, -1e-9), method="pyvista")
    distance, linecut_conf_pot = an.linecut(dvc.mesh, dvc.V/ct.e,
                    (0, ymin, -1e-9), (0, ymax, -1e-9), method="pyvista")
    potential_linecuts.append([distance, linecut_pot, linecut_conf_pot])

# save the linecuts
output_file = path_out / "potential_linecuts.npy"
np.save(output_file, np.array(potential_linecuts))

# Load the linecuts
potential_linecuts = np.load(output_file)

# Finding the time-dependent potential along the linecut
def confinement_potential(t):
    """
    Calculates the potential along the linecut.

    Args:
        t (float): Time in seconds.
    Returns:
        distance (1d array): Distance along the linecut.
        potential (2d array): Electrostatic potential linecut for the given t
        confinement_pot (1d array): Confinement potential linecut for the given t 
    """
    # Calculate the time-dependent part
    V_A0 = V_A_i(t, 0)
    V_A1 = V_A_i(t, 1)
    V_A2 = V_A_i(t, 2)

    # initialize the confinement potential
    potential = 0
    confinement_pot = 0
    for i in range(N):
        potential += V_A_i(t, i)*potential_linecuts[i, 1, :]
        confinement_pot += V_A_i(t, i)*potential_linecuts[i, 2, :]

    return distance, potential, confinement_pot


# calculate the confinement potential at t=tau/3
dist, phi_t, V_t = confinement_potential(tau/3)
# Plot the potential linecut
fig, ax = plt.subplots(figsize=(8,5))
ax.set_xlabel("Distance along linecut [nm]")
ax.set_ylabel("Electrostatic potential ($\\phi$) [mV]")
ax.plot(dist/1e-9, phi_t/1e-3, "-r")
ax.axhline(phi_t[0]/1e-3, color="black", lw=0.5, ls="--", alpha=0.5)

# plot the confinement potential
ax2 = ax.twinx()
ax2.set_ylabel("Confinement-potential ($V$) [meV]", color="blue")
ax2.plot(dist/1e-9, V_t/1e-3, "-b", label="Confinement potential")
ax2.tick_params(axis="y", labelcolor="blue")
ax2.axhline(V_t[0]/1e-3, color="black", lw=0.5, ls="--", alpha=0.5)

plt.show()

# Varying the gate voltages with respect to time 
time_steps = 10
time = np.linspace(tau/3, 2*tau/3, time_steps)

# plot the confinement potential linecuts
fig, ax = plt.subplots(figsize=(8, 5))
ax.set_xlabel("Distance along linecut [nm]")
ax.set_ylabel("confinement potential ($V$) [meV]")
for idt, t in enumerate(time):
    dist, phi_t, V_t = confinement_potential(t)
    ax.plot(dist/1e-9, V_t/1e-3, label=f"t = {t * 1e9:.1f} ns")
ax.annotate("",(13,30),(12,30), arrowprops=dict(facecolor='black', arrowstyle='-|>'))
ax.text(12.5, 40, "Time Progression", fontsize=10, ha='center', va='center')
ax.legend()
plt.show()

# Fit a second-order polynomial near the potential's minimum.
fit_range = 5e-9
# Lists to store the fitting parameters.  
coef_list = [] 
x_min_list = [] 
fig, ax = plt.subplots(figsize=(8, 5))
ax.set_xlabel("Distance along linecut [nm]")
ax.set_ylabel("confinement potential ($V$) [meV]")
for idt, t in enumerate(time):
    dist, _, V_t = confinement_potential(t)
    # Find the minimum of the linecut.
    min_idx = np.argmin(V_t)
    min_dist = dist[min_idx]

    # Select the range around the minimum for fitting
    fit_mask = (dist >= min_dist - fit_range) & (dist <= min_dist + fit_range)
    fit_dist = dist[fit_mask]
    fit_potential = V_t[fit_mask]

    # Fit a second-order polynomial
    coeffs = np.polyfit(fit_dist, fit_potential, 2)
    coef_list.append(coeffs)

    # create a polynomial function from the coefficients    
    poly = np.poly1d(coeffs)    

    # Calculate the approximated polynomial's minimum    
    x_min = -coeffs[1]/(2*coeffs[0])
    x_min_list.append(x_min)    # save the minimum of the polynomial
    y_min = poly(x_min)

    # Plot the results
    line = ax.plot(dist/1e-9, V_t/1e-3, label=f"t = {t * 1e9:.1f} ns")
    ax.plot(fit_dist/1e-9, poly(fit_dist)/1e-3,'--', color=line[0].get_color())
    ax.plot(x_min/1e-9, y_min/1e-3, "ko")
ax.annotate("", (13, 30), (12, 30), arrowprops=dict(facecolor='black', arrowstyle='-|>'))
ax.text(12.5, 40, "Time Progression", fontsize=10, ha='center', va='center')
ax.legend()
plt.show()

# Extract the position, speed and size of the quantum dot.
coef_list = np.array(coef_list)
x_min_list = np.array(x_min_list)
# Shift the position to start from 0
x_min_list = x_min_list - x_min_list[0] 
w = np.sqrt(2*ct.e*coef_list[:, 0]/(0.19*ct.me)) 
size = np.sqrt(ct.hbar/(0.19*ct.me*w)) # Size of the quantum dot
speed = np.gradient(x_min_list, time)

# Plot the position of the quantum dot over time
fig, ax = plt.subplots(figsize=(8, 5))
ax.set_xlabel("Time [ns]")
ax.set_ylabel("Position [nm]")
ax.fill_between(time*1e9, (x_min_list - size/2)/1e-9, (x_min_list + size/2)/1e-9, 
                color='gray', alpha=0.3, label="Quantum Dot Size")
ax.plot(time*1e9, x_min_list/1e-9, 'o-', label="Position")
ax.legend()
plt.show()

# Plot the speed of the quantum dot over time
fig, ax = plt.subplots(figsize=(8, 5))
ax.set_xlabel("Time [ns]")
ax.set_ylabel("Speed [m/s]")
ax.plot(time*1e9, speed, 'o-')
plt.show()

# Plot the Size of the quantum dot over time
fig, ax = plt.subplots(figsize=(8, 5))
ax.set_xlabel("Time [ns]")
ax.set_ylabel("Size [nm]")
ax.plot(time* 1e9, size/1e-9, 'o-')
plt.show()