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 consider a geometry similar to the one described in A double quantum dot device in a fully-depleted silicon-on-insulator transistor, but with three clavier gates instead of plunger and barrier gates. In addition to the modified gate configuration, this structure features two opposing boundaries across which we apply periodic boundary conditions. The full geometry is shown in Fig. 21.3.1.

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 apply periodic boundary conditions in QTCAD, the mesh must be periodic, meaning the node distributions on the corresponding boundary surfaces must be identical. This can be achieved in one of three ways in a 3D geometry: (1) by constructing a symmetric mesh, which is inherently periodic; (2) by using the Extrude command in Gmsh to generate the mesh along the direction connecting the centers of the periodic boundaries; or (3) by explicitly specifying periodicity in the Gmsh geometry file using the Periodic Surface keyword, 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. If these lists contain more than one element, they must be ordered such that each element in right_bnd[] is correctly paired with the element at the same index in left_bnd[], in a way that respects the specified translation.

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, boundary conditions are enforced such that the potential \(\varphi(\mathbf{r})\) at the surface \(\Sigma_j\) corresponding to the gate \(A_j\) is given by:

(21.5.1)\[\varphi(\mathbf{r}, t) = \phi_j(t) \equiv B \cos\left(\Omega t - \frac{2\pi j}{N}\right) \quad \forall \mathbf{r} \in \Sigma_j,\]

where \(B\) is the voltage amplitude, \(\Omega = 2\pi / \tau\) is the angular frequency of the time-dependent oscillation, \(j \in \{0, 1, 2\}\) is the index of the clavier gate, and \(N\) is the total number of gates. In addition we consider a back gate over which a constant bias is applied. Here, for convenience, we label the back gate as \(j=3\) and set \(\phi_{3}(t) = -0.5\) V.

# Define gate biases
back_gate_bias = -0.5
N = 3 # Number of gates
B = 1e-1 # 100 mV amplitude
tau = 1e-9 # 1 ns period
Omega = 2*np.pi/tau

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

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

    """
    return B*np.cos(Omega*t - ((2*np.pi*j)/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)

# 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 gates is turned on (\(\phi_{j}=1\)) and the potential of the other gates is set to zero (\(\phi_{k}=0\) for \(k\ne j\)). After repeating this process for all the gates, we can find the time-dependent potential, \(\varphi(\mathbf{r},t)\), using the following equation [JBF25]

(21.6.1)\[\varphi(\mathbf{r},t) = \sum_{j=0}^{N} \phi_{j}(t) \varphi_{j}(\mathbf{r}).\]

Here, \(\varphi_{j}(\mathbf{r})\) is the potential of the device when only gate \(j\) is activated and \(\phi_{j}(t)\) is the time-dependent bias applied to value of the potential enforced by the boundary condition. The indices \(j \in \{0,1,2\}\) correspond to the clavier gates \(A_j\) and \(j=3\) corresponds to the back gate. The goal is to compute the full potential \(\varphi(\mathbf{r},t)\) when \(\phi_{j}(t)\) are give by (21.5.1) for \(j \in \{0, 1, 2\}\) and \(\phi_{3}(t) = -0.5\) V.

In principle gate boundaries are appropriate to model the clavier gates, and a frozen boundary condition is appropriate for the back gate. However, to simulate the activation of each clavier gate and the back gate, we do not use new_gate_bnd or new_frozen_bnd boundary conditions, because these methods modify the applied potential at the boundary nodes according to additional physical properties. For example, new_gate_bnd includes contributions from the metal work function, as described in (1.11), and new_frozen_bnd similarly applies a modified potential based on material parameters, as described in (1.15). These modifications affect the actual potential applied to the boundary and interfere with the intended gate activation behavior.

To avoid these effects and ensure that the activation process is executed correctly, we use new_dirichlet_bnd during the activation simulations. This method directly sets the electrostatic potential at the boundary, which is necessary for the activation of individual gates with fixed voltages (e.g., 0 V or 1 V). However, we also define the input parameters required for later application of new_gate_bnd or new_frozen_bnd, which can be used once accurate modeling of gate stacks and substrate boundaries is needed.

# Define the gate biases for the first clavier gate activation
phi_A0 = 1
phi_A1 = 0
phi_A2 = 0
phi_back_gate = 0

# Apply Dirichlet boundary conditions
dvc.new_dirichlet_bnd("A0", phi_A0)
dvc.new_dirichlet_bnd("A1", phi_A1)
dvc.new_dirichlet_bnd("A2", phi_A2)
dvc.new_dirichlet_bnd("back_gate_bnd", phi_back_gate)

# Define parameters for frozen boundary to be applied later to the back gate
frozen_values = {
    'material': mt.Si,
    'doping': 1e15 * 1e6,           # Doping in cm⁻³ converted to m⁻³
    'doping_type': 'n',
    'binding_energy': 46e-3 * ct.e  # Binding energy in joules
}

# Define typical metallic work function (mid-gap) in joules
Ew = mt.Si.Eg / 2 + mt.Si.chi

Now that we have set up the device and defined the boundary conditions, we can proceed to solve the electrostatics of the device. We will use the linear Poisson solver to find the electrostatic potential 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 and the back-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), (0,1,0,0), (0,0,1,0), (0,0,0,1)]
Na = len(activate_bias)

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

# 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 {Na}')
    # Update gate biases
    dvc.set_applied_potential("A0", bias[0])
    dvc.set_applied_potential("A1", bias[1])
    dvc.set_applied_potential("A2", bias[2])
    dvc.set_applied_potential("back_gate_bnd", bias[3])

    # 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, height), (0, ymax, height), method="pyvista")
    potential_linecuts.append([distance, linecut_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 back gate 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. We then find the confinement potential \(V(\mathbf{r}, t)\) using the relation between the electrostatic potential and the confinement potential energy outlined in (2.2)

Since the gates are defined as Dirichlet boundary conditions, we need to modify the potential at the clavier gates to apply a gate boundary condition. By following Eqs. (1.11) and (1.15), we can find the correct electrostatic potential at the clavier gates and the back gate. The electrostatic potential, and similarly the confinement potential, of the device along the specified linecut at any time \(t\) can be calculated using the following function:

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

def confinement_potential(t, potential_linecuts=potential_linecuts,
    Ew=Ew, frozen_values=frozen_values):
    """
    Calculates the potential along the linecut.

    Args:
        t (float): Time in seconds.
        potential_linecuts (2d array): Pre-computed potential linecuts for
            each gate activation.
        Ew (float): Metallic work function in joules.
        frozen_values (dict): Values of the frozen boundary conditions.

    Returns:
        distance (1d array): Distance along the linecut in meters.
        potential (1d array): Electrostatic potential linecut at time t (in eV).
        confinement_pot (1d array): Confinement potential linecut at time t (in eV).
    """

    # Initialize potentials
    potential = 0
    # reference potential
    phi_ref = dvc.phi_F[0,0]  # Reference potential from the device
    for j in range(N):
        # Apply the gate boundary condition for each clavier gate
        potential += (phi_j(t, j) + phi_ref - Ew/ct.e) * potential_linecuts[j, 1, :]

    # Frozen boundary condition at the back gate
    T = 0.1  # Temperature in K
    gc = frozen_values['material'].gc
    mc = frozen_values['material'].mc
    chi = frozen_values['material'].chi

    doping = frozen_values['doping']
    binding_energy = frozen_values['binding_energy']

    x = gc * np.sqrt(mc * 3) * ((T)**1.5 * ct.C0)
    EC = binding_energy/ 2. + ct.kB * T / 2. * np.log(2 * x / doping)
    phi_bi = -(EC + chi) / ct.e

    # Apply back gate contribution
    potential += (back_gate_bias + phi_ref + phi_bi) * potential_linecuts[3, 1, :]
    confinement_pot = -ct.e *(potential-phi_ref) - chi

    # Return distance, electrostatic potential, and confinement potential
    return potential_linecuts[0, 0, :], potential, confinement_pot/ct.e

After computing the time-dependent electrostatic potential \(\varphi\) and the corresponding confinement potential \(V\), we plot them along the linecut at \(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 ($\\varphi$) [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:.2f} 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.2)\[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.3)\[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:.2f} 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. However, it is possible to avoid this approximation by analyzing the full three-dimensional confinement potential. Moreover, the position, speed, and size of the quantum dot can be computed more accurately by using the Schrödinger Solver and analyzing the resulting wavefunctions.

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

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
B = 1e-1 # 100 mV amplitude
tau = 1e-9 # 1 ns period
Omega = 2*np.pi/tau 

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

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

    """
    return B*np.cos(Omega*t - ((2*np.pi*j)/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)

# Periodic boundary conditions:
dvc.new_periodic_bnd("left_bnd", "right_bnd")

# Define the gate biases for the first clavier gate activation
phi_A0 = 1
phi_A1 = 0
phi_A2 = 0
phi_back_gate = 0

# Apply Dirichlet boundary conditions
dvc.new_dirichlet_bnd("A0", phi_A0)
dvc.new_dirichlet_bnd("A1", phi_A1)
dvc.new_dirichlet_bnd("A2", phi_A2)
dvc.new_dirichlet_bnd("back_gate_bnd", phi_back_gate)

# Define parameters for frozen boundary to be applied later to the back gate
frozen_values = {
    'material': mt.Si,
    'doping': 1e15 * 1e6,           # Doping in cm⁻³ converted to m⁻³
    'doping_type': 'n',
    'binding_energy': 46e-3 * ct.e  # Binding energy in joules
}

# Define typical metallic work function (mid-gap) in joules
Ew = mt.Si.Eg / 2 + mt.Si.chi

# 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), (0,1,0,0), (0,0,1,0), (0,0,0,1)]
Na = len(activate_bias)

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

# 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 {Na}')
    # Update gate biases
    dvc.set_applied_potential("A0", bias[0])
    dvc.set_applied_potential("A1", bias[1])
    dvc.set_applied_potential("A2", bias[2])
    dvc.set_applied_potential("back_gate_bnd", bias[3])

    # 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, height), (0, ymax, height), method="pyvista")
    potential_linecuts.append([distance, linecut_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)

def confinement_potential(t, potential_linecuts=potential_linecuts, 
    Ew=Ew, frozen_values=frozen_values):
    """
    Calculates the potential along the linecut.

    Args:
        t (float): Time in seconds.
        potential_linecuts (2d array): Pre-computed potential linecuts for 
            each gate activation.
        Ew (float): Metallic work function in joules.
        frozen_values (dict): Values of the frozen boundary conditions.

    Returns:
        distance (1d array): Distance along the linecut in meters.
        potential (1d array): Electrostatic potential linecut at time t (in eV).
        confinement_pot (1d array): Confinement potential linecut at time t (in eV).
    """

    # Initialize potentials
    potential = 0
    # reference potential
    phi_ref = dvc.phi_F[0,0]  # Reference potential from the device
    for j in range(N):
        # Apply the gate boundary condition for each clavier gate
        potential += (phi_j(t, j) + phi_ref - Ew/ct.e) * potential_linecuts[j, 1, :]

    # Frozen boundary condition at the back gate
    T = 0.1  # Temperature in K
    gc = frozen_values['material'].gc  
    mc = frozen_values['material'].mc
    chi = frozen_values['material'].chi

    doping = frozen_values['doping']
    binding_energy = frozen_values['binding_energy']

    x = gc * np.sqrt(mc * 3) * ((T)**1.5 * ct.C0)
    EC = binding_energy/ 2. + ct.kB * T / 2. * np.log(2 * x / doping)
    phi_bi = -(EC + chi) / ct.e

    # Apply back gate contribution
    potential += (back_gate_bias + phi_ref + phi_bi) * potential_linecuts[3, 1, :]
    confinement_pot = -ct.e *(potential-phi_ref) - chi

    # Return distance, electrostatic potential, and confinement potential
    return potential_linecuts[0, 0, :], potential, confinement_pot/ct.e


# 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 ($\\varphi$) [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:.2f} 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:.2f} 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()