# 4. Adaptive-mesh Poisson solver combined with the Schrödinger solver

## 4.1. Requirements

### 4.1.1. Software components

QTCAD

Gmsh

### 4.1.2. Geometry file

`qtcad/examples/tutorials/meshes/gaafet_adaptive.geo`

### 4.1.3. Python script

`qtcad/examples/tutorials/adaptive_schrodinger.py`

### 4.1.4. References

## 4.2. Briefing

In this tutorial, we revisit the GAAFET structure studied in Poisson and Schrödinger simulation of a GAAFET, Poisson solver with adaptive meshing, and Poisson solver with background charges. In Poisson solver with adaptive meshing, we saw that a relatively flat potential well (along the symmetry axis of the GAAFET) was formed below the Fermi level by the conduction band edge. Here, we will see that the conduction band edge was flat mostly because of screening effects from a classical electron gas forming inside the channel. Due to the nanometric dimensions of the GAAFET studied here, this approach is, in fact, inaccurate. This is because the statistical model used to populate the states of the GAAFET was classical, and thus neglected quantum confinement.

To study these confinement effects, we will again solve the non-linear Poisson equation self-consistently, but define a “quantum dot region” in which classical charges are set to zero, and where Schrödinger’s equation will be solved. Away from this dot region, charge carriers are still treated classically, enabling numerically efficient computations. In a later tutorial (Many-body analysis of a GAAFET), we will learn how Coulomb interactions between electrons inside the quantum dot region may be captured in a fully quantum-mechanical way.

## 4.3. Geometry file

As before, the initial mesh is generated from a
geometry file `qtcad/examples/tutorials/meshes/gaafet_adaptive.geo`

by running Gmsh with

```
gmsh gaafet_adaptive.geo
```

As in Poisson solver with adaptive meshing, this will produce both `.msh4`

and `.geo_unrolled`

files, as required by the adaptive-mesh Poisson solver for geometries
produced using the built-in Gmsh geometry kernel.

## 4.4. Setting up the device and the Poisson solver

### 4.4.1. Header and input parameters

The header and input parameters are almost the same as in the previous tutorial

```
from qtcad.device import constants as ct
from qtcad.device.mesh3d import Mesh, SubMesh
from qtcad.device import io
from qtcad.device import analysis
from qtcad.device import materials as mt
from qtcad.device import Device, SubDevice
from qtcad.device.poisson import Solver as PoissonSolver
from qtcad.device.poisson import SolverParams
from qtcad.device.schrodinger import Solver as SchrodingerSolver
import pathlib
```

However, three additional classes are loaded:

A

`SubMesh`

class which will be used to define a subset of the full mesh that will correspond to the dot region.A

`SubDevice`

class which will be used to define a subset of the full device that will correspond to the dot region.The Schrödinger

`Solver`

class.

We next define geometry parameters exactly as in the previous tutorial.

```
# Define some device parameters
Vgs = 2.0 # Gate-source bias
Ew = mt.Si.chi + mt.Si.Eg/2 # Metal work function
Ltot = 20e-9 # Device height in m
radius = 2.5e-9 # Device radius in m
```

Finally, we define paths to the mesh and raw geometry files.

```
# Paths to mesh file and output files
script_dir = pathlib.Path(__file__).parent.resolve()
path_mesh = script_dir / "meshes" / "gaafet_adaptive.msh4"
path_geo = script_dir / "meshes" / "gaafet_adaptive.geo_unrolled"
```

### 4.4.2. Loading the initial mesh and defining the device

We load the mesh and define the device exactly as in the previous tutorial.

```
# Load the mesh
scaling = 1e-9
mesh = Mesh(scaling, path_mesh)
# Create device from mesh and set temperature to 10 mK
d = Device(mesh, conf_carriers="e")
d.set_temperature(10e-3)
# Create device regions
d.new_region('channel_ox',mt.SiO2)
d.new_region('source_ox',mt.SiO2)
d.new_region('drain_ox',mt.SiO2)
d.new_region('channel', mt.Si)
d.new_region('source', mt.Si, pdoping=5e20*1e6, ndoping=0)
d.new_region('drain', mt.Si, pdoping=5e20*1e6, ndoping=0)
# Create device boundaries
d.new_ohmic_bnd('source_bnd')
d.new_ohmic_bnd('drain_bnd')
d.new_gate_bnd('gate_bnd',Vgs,Ew)
```

### 4.4.3. Defining the quantum dot region

The dot region will be defined by a union of region labels, here represented by a list.

```
# Define the dot region as a list of region labels
dot_region = ["channel", "channel_ox"]
```

It is important to include both the channel and the oxide surrounding it in the dot region because we will need to include the interface between these regions when solving Schrödinger’s equation. Indeed, this interface defines the potential barriers which confine the electrons inside the channel.

We then define the dot region in which no classical charge is allowed using
the `set_dot_region`

of the `Device`

class. This must be done
before solving Poisson’s equation—otherwise the Poisson
solver may put classical charges inside the quantum dot which will lead to
screening effects that would not be there if the quantum dot were to contain
just a few confined electrons.

```
# Set up the dot region in which no classical charge is allowed
d.set_dot_region(dot_region)
```

### 4.4.4. Setting up the adaptive-mesh Poisson solver

The next step is to set up the parameters of the adaptive-mesh Poisson solver. Here we proceed exactly as in the previous tutorial, except that we set the minimum number of nodes to 10,000 to ensure that the resolution of the final mesh is satisfactory.

In addition, we define two important solver parameters
(see the Poisson `SolverParams`

API reference for a more extensive description of these parameters):

`dot_region`

: a list of strings labeling a collection of regions in which we want to set a maximal value for the characteristic length.`h_dot`

: the maximal value for the characteristic length.

Importantly, setting the `dot_region`

solver parameter does not set classical
charge to zero (this was done above with the
`set_dot_region`

of the `Device`

class). Instead, the solver
parameters introduced here enable to force the adaptive meshing algorithm to
keep the mesh fine enough in the dot region for a satisfactory accuracy in
the solution of the Schrödinger equation which will be obtained later.

```
# Configure the non-linear Poisson solver
params_poisson = SolverParams()
params_poisson.tol = 1e-5 # Convergence threshold (tolerance) for the error
params_poisson.initial_ref_factor = 0.1
params_poisson.final_ref_factor = 0.75
params_poisson.min_nodes = 10000
params_poisson.dot_region = dot_region
params_poisson.h_dot = 0.25
```

Note

The parameter `h_dot`

is expressed in the units of length assumed in the
Gmsh geometry file, and *not* in SI units. Here, since we took
`scaling=1e-9`

when creating the mesh, this will mean that the maximum
characteristic length in the dot region is \(0.25\;\mathrm{nm}\).

After defining Poisson solver parameters, we instantiate the Poisson solver object.

```
# Create an adaptive-mesh non-linear Poisson solver
s = PoissonSolver(d, solver_params=params_poisson, geo_file=path_geo)
```

## 4.5. Solving the Poisson and Schrödinger equations

We launch the Poisson simulation with:

```
# Self-consistent solution
s.solve()
```

After convergence of the Poisson solver, we plot the band diagram along the symmetry axis of the GAAFET.

```
# Plot band diagram
analysis.plot_bands(d, (0,0,0), (0,0,Ltot))
```

We note two main differences with the band diagram obtained in Poisson solver with adaptive meshing. (i) The potential well is much deeper. (ii) The bottom of the potential well has a parabolic-like shape. Both differences are related to the fact that screening from a classical gas of electrons in the channel is now absent.

Using the solution of the non-linear Poisson equation, we then define the confinement potential energy which will be used to solve Schrödinger’s equation.

```
# Get the potential energy from the band edge for usage in the Schrodinger
# solver
d.set_V_from_phi()
```

The next step is to create a subset of the device in which Schrödinger’s equation will be solved.

We first create a `SubMesh`

object which contains only the nodes inside the dot region.

```
# Create a submesh including only the dot region
submesh = SubMesh(d.mesh, dot_region)
submesh.show()
```

The second command above displays the submesh.

This submesh is non-uniform (finer at the top and bottom boundaries), since it accounts for both the automatic mesh refinements and the maximum characteristic lengths set in the adaptive-mesh Poisson solver parameters.

Once the submesh is created, we define the subset of the device which lives over this submesh.

```
# Create a subdevice for the dot region
subdevice = SubDevice(d,submesh)
```

We then instantiate the Schrödinger equation over the subdevice, and solve.

```
# Create a Schrodinger solver
schrod_solver = SchrodingerSolver(subdevice)
# Solve Schrodinger's equation
schrod_solver.solve()
```

We print the energy levels obtained from the numerical solution of
Schrödinger’s equation (by default, the first ten levels—see the Schrödinger
`SolverParams`

API
reference to account for a different number of energy levels).

```
# Print energies
subdevice.print_energies()
```

The result is

```
Energy levels (eV)
[-0.89616905 -0.81939494 -0.73555787 -0.64530374 -0.54797826 -0.44382518
-0.33280352 -0.29090282 -0.29038706 -0.21040997]
```

We finally plot the first three eigenfunctions.

```
# Plot the first three eigenfunctions
analysis.plot_slices(submesh, subdevice.eigenfunctions[:,0],
title="Ground state")
analysis.plot_slices(submesh, subdevice.eigenfunctions[:,1],
title="First excited state")
analysis.plot_slices(submesh, subdevice.eigenfunctions[:,2],
title="Second excited state")
```

## 4.6. Full code

```
__copyright__ = "Copyright 2024, Nanoacademic Technologies Inc."
from qtcad.device import constants as ct
from qtcad.device.mesh3d import Mesh, SubMesh
from qtcad.device import io
from qtcad.device import analysis
from qtcad.device import materials as mt
from qtcad.device import Device, SubDevice
from qtcad.device.poisson import Solver as PoissonSolver
from qtcad.device.poisson import SolverParams
from qtcad.device.schrodinger import Solver as SchrodingerSolver
import pathlib
# Define some device parameters
Vgs = 2.0 # Gate-source bias
Ew = mt.Si.chi + mt.Si.Eg/2 # Metal work function
Ltot = 20e-9 # Device height in m
radius = 2.5e-9 # Device radius in m
# Paths to mesh file and output files
script_dir = pathlib.Path(__file__).parent.resolve()
path_mesh = script_dir / "meshes" / "gaafet_adaptive.msh4"
path_geo = script_dir / "meshes" / "gaafet_adaptive.geo_unrolled"
# Load the mesh
scaling = 1e-9
mesh = Mesh(scaling, path_mesh)
# Create device from mesh and set temperature to 10 mK
d = Device(mesh, conf_carriers="e")
d.set_temperature(10e-3)
# Create device regions
d.new_region('channel_ox',mt.SiO2)
d.new_region('source_ox',mt.SiO2)
d.new_region('drain_ox',mt.SiO2)
d.new_region('channel', mt.Si)
d.new_region('source', mt.Si, pdoping=5e20*1e6, ndoping=0)
d.new_region('drain', mt.Si, pdoping=5e20*1e6, ndoping=0)
# Create device boundaries
d.new_ohmic_bnd('source_bnd')
d.new_ohmic_bnd('drain_bnd')
d.new_gate_bnd('gate_bnd',Vgs,Ew)
# Define the dot region as a list of region labels
dot_region = ["channel", "channel_ox"]
# Set up the dot region in which no classical charge is allowed
d.set_dot_region(dot_region)
# Configure the non-linear Poisson solver
params_poisson = SolverParams()
params_poisson.tol = 1e-5 # Convergence threshold (tolerance) for the error
params_poisson.initial_ref_factor = 0.1
params_poisson.final_ref_factor = 0.75
params_poisson.min_nodes = 10000
params_poisson.dot_region = dot_region
params_poisson.h_dot = 0.25
# Create an adaptive-mesh non-linear Poisson solver
s = PoissonSolver(d, solver_params=params_poisson, geo_file=path_geo)
# Self-consistent solution
s.solve()
# Plot band diagram
analysis.plot_bands(d, (0,0,0), (0,0,Ltot))
# Get the potential energy from the band edge for usage in the Schrodinger
# solver
d.set_V_from_phi()
# Create a submesh including only the dot region
submesh = SubMesh(d.mesh, dot_region)
submesh.show()
# Create a subdevice for the dot region
subdevice = SubDevice(d,submesh)
# Create a Schrodinger solver
schrod_solver = SchrodingerSolver(subdevice)
# Solve Schrodinger's equation
schrod_solver.solve()
# Print energies
subdevice.print_energies()
# Plot the first three eigenfunctions
analysis.plot_slices(submesh, subdevice.eigenfunctions[:,0],
title="Ground state")
analysis.plot_slices(submesh, subdevice.eigenfunctions[:,1],
title="First excited state")
analysis.plot_slices(submesh, subdevice.eigenfunctions[:,2],
title="Second excited state")
```