22. Symmetric meshing with Poisson and Schrödinger solvers

22.1. Requirements

22.1.1. Software components

  • QTCAD

  • Gmsh

22.1.2. 3.1.2. Geometry file

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

22.1.3. 3.1.3. Python script

  • qtcad/examples/tutorials/symmetric_meshing.py

22.1.4. 3.1.4. References

—

22.2. Briefing

In this tutorial, we introduce symmetric meshing, a feature which allows for the automatic generation of a full mesh from a half mesh by mirroring it across a specified symmetry plane. This process can be used e.g. to eliminate artificial detuning steps introduced in Tuning the double quantum dot into a symmetric configuration.

We start by using the QTCAD API to generate a symmetric mesh by mirroring an existing three-dimensional mesh, which represents half of the double quantum dot FD-SOI device introduced in A double quantum dot device in a fully-depleted silicon-on-insulator transistor, across a chosen symmetry plane. Next, we set up and solve the non-linear Poisson equation with adaptive mesh refinement. Finally, we solve the Schrödinger equation over a submesh representing the double quantum dot region to find the ground and first excited states. From these results, we extract the tunnel coupling between the dots and compare it with the asymmetric mesh case presented in Tunnel coupling in a double quantum dot in FD-SOI—Part 1: Plunger gate tuning. This comparison shows that both approaches give consistent results.

22.3. Geometry file

As shown before, the initial half mesh is generated from a geometry file qtcad/examples/tutorials/meshes/qdfdsoi.geo by using Gmsh with:

gmsh qdfdsoi.geo

This process generates a .msh file for the initial mesh and a .xao file, which is later used for the adaptive meshing step.

22.5. Symmetric meshing

In QTCAD, any mesh, whether one, two, or three dimensional, can be mirrored across the corresponding symmetry element: a point in one dimension, a line in two dimensions, or a plane in three dimensions. This is done by using the symmetric_mesh method of the Mesh class.

The workflow involves identifying the point, line, or plane of symmetry, generating the half mesh, importing it into QTCAD, and then calling symmetric_mesh to create the full symmetric mesh. In this tutorial, we mirror half of the double quantum dot FD-SOI device introduced in A double quantum dot device in a fully-depleted silicon-on-insulator transistor. The half mesh is stored in qdfdsoi.msh as follows:

scaling = 1e-9
file = path_mesh / "qdfdsoi.msh"
hmesh = Mesh(scaling, file)
hmesh.show()

We then mirror the mesh across the \(y = 32.5 \, \text{nm}\) plane. For a three-dimensional mesh, this plane is expressed in the general form \(ax + by + cz + d = 0\) with coefficients \(a = 0\), \(b = 1\), \(c = 0\), and \(d = -32.5 \, \text{nm}\).

mesh = hmesh.symmetric_mesh(
    0, 1, 0, -32.5e-9,
    mesh_file=str(path_mesh / "qdfdsoi_mirrored.msh"),
    geo_file=str(path_geo),
    verbose=True,
)
mesh.show()

This step automatically generates mirrored regions and boundaries, which are identical to those in the original mesh but with _mirrored appended to their names. If the paths are provided, the mirrored mesh and geometry will be saved to the specified locations using the mesh_file and geo_file arguments, respectively.

Note

Only the .xao format is supported for saving the mirrored geometry. The geo_file argument must specify the path to the .xao file of the half-geometry. The mirrored geometry will be saved with an _mirrored.xao suffix.

The following is an example of the mesh statistics output for the mirrored mesh:

Loading mesh from file:
qtcad\examples\tutorials\meshes\qdfdsoi.msh
Done.
================================================================================
3D MESH STATISTICS
--------------------------------------------------------------------------------
Total number of elements      83368
Triangular elements           748
Tetrahedral elements          82620
Total number of nodes         15743
Boundary physical names       barrier_gate_1_bnd, plunger_gate_1_bnd, barrier_gate_2_bnd,
                               source_bnd, back_gate_bnd,
Region physical names         gate_oxide_dot, gate_oxide, source, channel, oxide,
                              channel_dot, oxide_dot, buried_oxide, buried_oxide_dot,
================================================================================
Processing and mirroring the half mesh.
Done.
================================================================================
3D MESH STATISTICS
--------------------------------------------------------------------------------
Total number of elements      166736
Triangular elements           1496
Tetrahedral elements          165240
Total number of nodes         30576
Boundary physical names       barrier_gate_1_bnd, plunger_gate_1_bnd, barrier_gate_2_bnd,
                              source_bnd, back_gate_bnd, barrier_gate_1_bnd_mirrored,
                              plunger_gate_1_bnd_mirrored, barrier_gate_2_bnd_mirrored,
                              source_bnd_mirrored, back_gate_bnd_mirrored,
Region physical names         gate_oxide_dot, gate_oxide, source, channel, oxide,
                              channel_dot, oxide_dot, buried_oxide, buried_oxide_dot,
                              gate_oxide_dot_mirrored, gate_oxide_mirrored, source_mirrored,
                              channel_mirrored, oxide_mirrored, channel_dot_mirrored,
                              oxide_dot_mirrored, buried_oxide_mirrored,
                              buried_oxide_dot_mirrored,
================================================================================

This mirrored mesh can then be used to run static simulations as covered in Poisson and Schrödinger simulation of a nanowire quantum dot. In this tutorial, however, we focus on using adaptive meshing. Fig. 22.5.1 (a) shows the original half mesh, while Fig. 22.5.1 (b) shows the mirrored symmetric mesh.

Comparison between original and mirrored mesh

Fig. 22.5.1 Comparison of double quantum dot meshes. (a) Original half mesh, (b) Symmetric mirrored mesh.

22.6. Setting up the device

Next, we define variables for the biases applied on the gates.

# Define gate biases
back_gate_bias = -0.5
barrier_gate_1_bias = 0.5
plunger_gate_1_bias = 0.59
barrier_gate_2_bias = 0.51

We note that since the gate biases are applied symmetrically in this specific example, there is no need to define mirrored gate biases. The solver will automatically apply the same bias to the mirrored gates. We then create a Device object and define the regions for the half mesh, following the same procedure as in A double quantum dot device in a fully-depleted silicon-on-insulator transistor.

dvc = Device(hmesh, conf_carriers="e")
dvc.set_temperature(0.1)

# Create the regions
dvc.new_region("oxide", mt.SiO2)
dvc.new_region("oxide_dot", mt.SiO2)
dvc.new_region("gate_oxide", mt.HfO2)
dvc.new_region("gate_oxide_dot", mt.HfO2)
dvc.new_region("buried_oxide", mt.SiO2)
dvc.new_region("buried_oxide_dot", mt.SiO2)
dvc.new_region("channel", mt.Si)
dvc.new_region("channel_dot", mt.Si)
dvc.new_region("source", mt.Si, ndoping=1e20*1e6)
# (mirrored regions are omitted)

# Boundary conditions
Ew = mt.Si.Eg/2 + mt.Si.chi

dvc.new_gate_bnd("barrier_gate_1_bnd", barrier_gate_1_bias, Ew)
dvc.new_gate_bnd("plunger_gate_1_bnd", plunger_gate_1_bias, Ew)
dvc.new_gate_bnd("barrier_gate_2_bnd", barrier_gate_2_bias, Ew)
dvc.new_ohmic_bnd("source_bnd")
dvc.new_frozen_bnd("back_gate_bnd", back_gate_bias, mt.Si, 1e15*1e6,
   "n", 46*1e-3*ct.e)

Although we only define the regions and boundaries for the half mesh, when using adaptive meshing the mirrored regions and boundaries are automatically created by the Poisson class.

It is still possible to run a conventional simulation by manually defining the mirrored regions and boundaries (similar to Poisson and Schrödinger simulation of a nanowire quantum dot). However, in that case adaptive refinement is not supported, and only static simulations can be performed.

The dot region is then defined as before:

# Mark quantum dot region
dot_region_list = [
    "oxide_dot", "gate_oxide_dot", "buried_oxide_dot", "channel_dot"]
dvc.set_dot_region(dot_region_list)

22.7. Adaptive non-linear Poisson solver

After setting up the device, we configure the Poisson solver with adaptive refinement. Symmetric meshing is applied and automatically handled by the solver when the mirror_plane_coefs parameter is specified in the solver parameters. For a three-dimensional mesh, this argument defines the plane of symmetry and is passed as a tuple (a, b, c, d), corresponding to the plane equation \(ax + by + cz + d = 0\). The solver automatically generates and defines the mirrored regions and boundaries, adapts the mesh accordingly, and solves the Poisson equation.

# Configure the non-linear Poisson solver
params_poisson = PoissonSolverParams()
params_poisson.tol = 1e-4
params_poisson.initial_ref_factor = 0.1
params_poisson.final_ref_factor = 0.75
params_poisson.min_nodes = 60000
params_poisson.max_nodes = 1e5
params_poisson.maxiter_adapt = 30
params_poisson.refined_region = dot_region_list
params_poisson.h_refined = 0.7
params_poisson.mirror_plane_coefs = (0, 1, 0, -32.5e-9)
params_poisson.refined_mesh_filename = path_mesh / "sym_refined_dqdfdsoi.msh"

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

# Solve Poisson equation
poisson_slv.solve()

As in the adaptive tutorial, the solver refines the mesh automatically if convergence is not achieved, with additional refinements concentrated near the quantum-dot region (controlled by the refined_region and h_refined solver parameters).

Note

After performing the corresponding mesh refinement, the adaptive Poisson solver automatically generates the full symmetric mesh. Using the symmetric geometry and mesh produced by the mesh.symmetric_mesh method, users can also perform a manual adaptive meshing simulation, similar to Poisson solver with adaptive meshing tutorial.

Next, we perform a linecut analysis of the conduction band edge along the channel.

# 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_qdot = 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_qdot/ct.e)
fig.savefig(path_out / "sym_cbe.png")
plt.show()

The above code produces the conduction band profile across the symmetric channel, as shown in Figure Fig. 22.7.1.

Conduction band edge along the channel

Fig. 22.7.1 Conduction band edge along the channel for the symmetric double quantum dot FD-SOI device.

22.8. Schrödinger solution

We then solve the Schrödinger equation in the dot region:

# Define the mirrored dot regions
dot_region_mirrored_list = ["oxide_dot", "gate_oxide_dot",
"buried_oxide_dot", "channel_dot", "oxide_dot_mirrored", "gate_oxide_dot_mirrored",
"buried_oxide_dot_mirrored", "channel_dot_mirrored"]

# Create a submesh including only the dot region
submesh = SubMesh(dvc.mesh, dot_region_mirrored_list)

# Instantiate Schrodinger solver parameters
params_schrod = SchrodingerSolverParams()
params_schrod.tol = 1e-7 # Tolerance on energies in eV

# Calculating the confinement potential
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()

After solving the Schrödinger equation, we save and visualize the ground and the first excited state wavefunctions of the double quantum dot (both full visualization and linecuts).

# Plot the ground state wavefunction
an.plot_slices(submesh, subdvc.eigenfunctions[:,0], z= -2e-9,
   title="Ground state wavefunction", show_figure=False,
   path=path_out / "sym_ground_state.png")

# Linecut for ground and excited state
distance, linecut_ground = an.linecut(submesh, subdvc.eigenfunctions[:,0],
   (0, ymin, -2e-9), (0, ymax, -2e-9))
distance, linecut_excited = an.linecut(submesh, subdvc.eigenfunctions[:,1],
   (0, ymin, -2e-9), (0, ymax, -2e-9))

fig = plt.figure(figsize=(8,5))
ax = fig.add_subplot(1,1,1)
ax.set_xlabel("Distance along linecut (nm)")
ax.set_ylabel("Eigenfunctions (arb. units)")
ax.plot(distance/1e-9, linecut_ground, "ro-",
   label=f"Ground state wavefunction")
ax.plot(distance/1e-9, linecut_excited, "bx-",
   label=f"First excited state wavefunction")
ax.legend()
fig.savefig(path_out / "sym_eigenfunctions.png")
plt.show()

Figure Fig. 22.8.1 and Fig. 22.8.2 show the ground and excited state wavefunctions, respectively. We note that between different runs, the results for the excited states may differ from those shown in Fig. 22.8.2 by a global phase factor, which does not affect the physical interpretation.

Full ground state of the double quantum dot.

Fig. 22.8.1 Double quantum dot ground state.

Linecuts for ground and excited states.

Fig. 22.8.2 Linecuts for ground and first excited states along the channel. This figure shows that the ground state and first excited state are delocalized across both dots.

As the final step, we calculate the tunnel coupling from the energy difference between the first excited state and the ground state:

# Tunnel coupling calculation
energies = subdvc.energies
tunnel_coupling = np.abs(energies[1] - energies[0])/ct.e
print(f'Tunnel coupling: {tunnel_coupling} eV')

The corresponding tunnel coupling is printed to the console. For the chosen gate biases, the tunnel coupling is:

Tunnel coupling: 1.27e-06 eV

showing close agreement with the results obtained in Tunnel coupling in a double quantum dot in FD-SOI—Part 1: Plunger gate tuning.

22.9. 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 io
from qtcad.device import analysis as an
from qtcad.device.mesh3d import Mesh, SubMesh
from qtcad.device import Device
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

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


scaling = 1e-9
file = path_mesh / "qdfdsoi.msh"
hmesh = Mesh(scaling, file)
hmesh.show()

mesh = hmesh.symmetric_mesh(
    0, 1, 0, -32.5e-9,
    mesh_file=str(path_mesh / "qdfdsoi_mirrored.msh"),
    geo_file=str(path_geo),
    verbose=True,
)
mesh.show()

# Define gate biases
back_gate_bias = -0.5
barrier_gate_1_bias = 0.5
plunger_gate_1_bias = 0.59
barrier_gate_2_bias = 0.51

dvc = Device(hmesh, conf_carriers="e")
dvc.set_temperature(0.1)

# Create the regions
dvc.new_region("oxide", mt.SiO2)
dvc.new_region("oxide_dot", mt.SiO2)
dvc.new_region("gate_oxide", mt.HfO2)
dvc.new_region("gate_oxide_dot", mt.HfO2)
dvc.new_region("buried_oxide", mt.SiO2)
dvc.new_region("buried_oxide_dot", mt.SiO2)
dvc.new_region("channel", mt.Si)
dvc.new_region("channel_dot", mt.Si)
dvc.new_region("source", mt.Si, ndoping=1e20*1e6)
# (mirrored regions are omitted)

# Boundary conditions
Ew = mt.Si.Eg/2 + mt.Si.chi

dvc.new_gate_bnd("barrier_gate_1_bnd", barrier_gate_1_bias, Ew)
dvc.new_gate_bnd("plunger_gate_1_bnd", plunger_gate_1_bias, Ew)
dvc.new_gate_bnd("barrier_gate_2_bnd", barrier_gate_2_bias, Ew)
dvc.new_ohmic_bnd("source_bnd")
dvc.new_frozen_bnd("back_gate_bnd", back_gate_bias, mt.Si, 1e15*1e6, 
   "n", 46*1e-3*ct.e)

# Mark quantum dot region
dot_region_list = [
    "oxide_dot", "gate_oxide_dot", "buried_oxide_dot", "channel_dot"]
dvc.set_dot_region(dot_region_list)

# Configure the non-linear Poisson solver
params_poisson = PoissonSolverParams()
params_poisson.tol = 1e-4
params_poisson.initial_ref_factor = 0.1
params_poisson.final_ref_factor = 0.75
params_poisson.min_nodes = 60000
params_poisson.max_nodes = 1e5
params_poisson.maxiter_adapt = 30
params_poisson.refined_region = dot_region_list
params_poisson.h_refined = 0.7
params_poisson.mirror_plane_coefs = (0, 1, 0, -32.5e-9)
params_poisson.refined_mesh_filename = path_mesh / "sym_refined_dqdfdsoi.msh"

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

# Solve Poisson equation
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_qdot = 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_qdot/ct.e)
fig.savefig(path_out / "sym_cbe.png")
plt.show()

# Define the mirrored dot regions
dot_region_mirrored_list = ["oxide_dot", "gate_oxide_dot", 
"buried_oxide_dot", "channel_dot", "oxide_dot_mirrored", "gate_oxide_dot_mirrored", 
"buried_oxide_dot_mirrored", "channel_dot_mirrored"]

# Create a submesh including only the dot region
submesh = SubMesh(dvc.mesh, dot_region_mirrored_list)

# Instantiate Schrodinger solver parameters
params_schrod = SchrodingerSolverParams()
params_schrod.tol = 1e-7 # Tolerance on energies in eV

# Calculating the confinement potential
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()

# Plot the ground state wavefunction
an.plot_slices(submesh, subdvc.eigenfunctions[:,0], z= -2e-9, 
   title="Ground state wavefunction", show_figure=False, 
   path=path_out / "sym_ground_state.png")

# Linecut for ground and excited state
distance, linecut_ground = an.linecut(submesh, subdvc.eigenfunctions[:,0],
   (0, ymin, -2e-9), (0, ymax, -2e-9))
distance, linecut_excited = an.linecut(submesh, subdvc.eigenfunctions[:,1],
   (0, ymin, -2e-9), (0, ymax, -2e-9))

fig = plt.figure(figsize=(8,5))
ax = fig.add_subplot(1,1,1)
ax.set_xlabel("Distance along linecut (nm)")
ax.set_ylabel("Eigenfunctions (arb. units)")
ax.plot(distance/1e-9, linecut_ground, "ro-", 
   label=f"Ground state wavefunction")
ax.plot(distance/1e-9, linecut_excited, "bx-", 
   label=f"First excited state wavefunction")
ax.legend()
fig.savefig(path_out / "sym_eigenfunctions.png")
plt.show()

# Tunnel coupling calculation
energies = subdvc.energies
tunnel_coupling = np.abs(energies[1] - energies[0])/ct.e
print(f'Tunnel coupling: {tunnel_coupling} eV')