# 13. Many-body analysis of a GAAFET

## 13.1. Requirements

### 13.1.1. Software components

qtcad

gmsh

### 13.1.2. Geometry file

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

### 13.1.3. Python script

`qtcad/examples/tutorials/manybody.py`

### 13.1.4. References

## 13.2. Briefing

In this tutorial, we build on the results output in the Poisson and Schrödinger simulation of a GAAFET
tutorial (see Python script `qtcad/examples/tutorials/gaafet.py`

).
Loading the potential obtained at the end of this file, we solve the
many-body problem for a system consisting of a few electrons confined in the
potential well formed by the conduction band edge within the channel of the
gate-all-around field-effect transistor (GAAFET). We then illustrate the
various physical outputs that can be extracted from the many-body solution.

## 13.3. Setting up the device and loading results from the previous simulation

### 13.3.1. Header

The header for this tutorial is

```
from qtcad.device import constants as ct
from qtcad.device.mesh3d import Mesh
from qtcad.device import io
from qtcad.device import analysis as an
from qtcad.device import materials as mt
from qtcad.device import Device
from qtcad.device.schrodinger import Solver as SingleBodySolver
import pathlib
from qtcad.device.many_body import Solver as ManyBodySolver
from qtcad.device.many_body import SolverParams
```

The only new additions with respect to the Poisson and Schrödinger simulation of a GAAFET tutorial are the
last two lines of the code block, in which the many-body
`Solver`

and
`SolverParams`

classes are imported from the
`many_body`

module.

### 13.3.2. Creating the device

We next load the mesh, instantiate the device, and set up regions and boundary conditions exactly as in the previous tutorial.

```
# Define some device parameters
Vgs = 0.5 # 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.msh"
path_hdf5 = script_dir / "output/" / "gaafet.hdf5"
# Load the mesh
scaling = 1e-9
mesh = Mesh(scaling, str(path_mesh))
# Create device from mesh and set statistics.
# It is important to set statistics before creating boundaries since
# they determine the boundary values, Failing to do so will produce
# inconsistent results.
dvc = Device(mesh, conf_carriers="e")
dvc.statistics = "FD_approx" # Analytic approximation to Fermi-Dirac statistics
# Create regions first.
# Make sure each element is assigned to a region, otherwise default
# (silicon) parameters are used
dvc.new_region('channel_ox',mt.SiO2)
dvc.new_region('source_ox',mt.SiO2)
dvc.new_region('drain_ox',mt.SiO2)
dvc.new_region('channel', mt.Si)
dvc.new_region('source', mt.Si, pdoping=5e20*1e6, ndoping=0)
dvc.new_region('drain', mt.Si, pdoping=5e20*1e6, ndoping=0)
# Then create boundaries
dvc.new_ohmic_bnd('source_bnd')
dvc.new_ohmic_bnd('drain_bnd')
dvc.new_gate_bnd('gate_bnd',Vgs,Ew)
```

The next step is to load the electric potential from the hdf5 output file produced in the previous tutorial.

```
# Load the electic potential
dvc.set_potential(io.load(path_hdf5,"phi"))
```

It is important to use the
`set_potential`

method of the
`Device`

class rather than modifying the
`phi`

attribute directly because the
method automatically recalculates the potential energy
`V`

.

## 13.4. Solving the many-body problem

We first solve the single-electron problem, which will define the single-particle basis set for the many-body problem

```
# Solve the single-electron problem
slv = SingleBodySolver(dvc)
slv.solve()
dvc.print_energies()
```

The single-particle energies will be output as

```
Energy levels (eV)
[0.54798018 0.61320021 0.68496646 0.76317642 0.8478653 0.93920254
1.03757746 1.14369621 1.15198013 1.15240696]
```

We then instantiate and run the many-body
`Solver`

.

```
# Solve the many-body problem
solver_params = SolverParams()
solver_params.num_states = 3
solver_params.n_degen = 2
solver_params.alpha = 0.5
solver_params.num_particles = [0,1,2,3]
slv = ManyBodySolver(dvc, solver_params = solver_params)
slv.solve()
```

Here we limit the number of single-particle orbitals to consider in the basis
set to \(3\), and set the degeneracy of these orbitals to
\(2`\) (e.g. from spin). In addition, we introduce a lever arm
\(\alpha\) of \(0.5\). The lever arm quantifies the extent by which
a potential applied at a gate can shift the energy-level structure of the
device (see the lever arm Theory section). Finally,
we set the `num_particles`

attribute of the many-body
`SolverParams`

to a list of
integers \(N\) labeling the
\(N\)-particle subspaces to consider. Here, we choose to only consider the
0-, 1-, 2-, and 3-particle subspaces even though for
\(n_\mathrm{states} = 3 ` and :math:`n_\mathrm{degen} = 2\),
the maximum number of particles that can be considered in principle is
\(n_\mathrm{states} \times n_\mathrm{degen} = 6\).
This choice is up to the user, but taking a smaller number of particles than
the total accessible amount will typically make calculations faster.

After calling the
`solve`

method,
a progress bar appears to track the calculation of Coulomb integrals

```
Computing Coulomb integrals ████████████████████████████████ 100%
```

## 13.5. Analyzing the results

Once the calculation is done, we produce various outputs that are made
available through quantities stored as new
`Device`

attributes.

### 13.5.1. Many-body subspace properties

We first analyze the many-body subspaces with

```
# Many-body subspace properties
print("Many-body subspaces:", dvc.many_body_subspaces)
print("Number of subspaces:", len(dvc.many_body_subspaces))
print("Electron number in subspace 2:", dvc.many_body_subspaces[2].N)
```

which results in

```
Many-body subspaces: [<device.core.quantum.Subspace object at 0x00000215769AD6A0>, <device.core.quantum.Subspace object at 0x0000021576CC1E80>, <device.core.quantum.Subspace object at 0x0000021576CC1FA0>, <device.core.quantum.Subspace object at 0x0000021576CC1E20>]
Number of subspaces: 4
Electron number in subspace 2: 2
```

In general, for \(n_\mathrm{states}=3\) single-particle orbitals with degeneracy \(n_\mathrm{degen}=2\), we have \(n_\mathrm{states}\times n_\mathrm{degen}+1=7\) subspaces corresponding to \(0\) to \(N_\mathrm{max}=n_\mathrm{states}\times n_\mathrm{degen}=6\) electrons (see Exact diagonalization). However, here, we have kept only the 0, 1, 2, and 3-electron subspaces. With these subspace integer labels starting at 0, the many-body subspace with index 2 is indeed expected to involve 2 electrons, as seen above.

Each subspace contains a many-body basis set. For example, the two-electron subspace basis set is displayed with

```
print("Two-electron many-body basis set (integers): ")
print(dvc.get_N_particle_subspace(2).get_bas_set())
```

which results in

```
Two-electron many-body basis set (integers):
[ 3 5 9 17 33 6 10 18 34 12 20 36 24 40 48]
```

By default, each many-body basis state is labeled by an integer whose binary
representation gives the occupation (0 if unoccupied, 1 if occupied) of each
single-partical orbital (see the
`get_bas_set`

API reference). This is the most compact representation available in QTCAD
in terms of memory usage. The strings for these binaries may be displayed with:

```
print("Two-electron many-body basis set (binary strings): ")
print(dvc.get_N_particle_subspace(2).get_bas_set(dtype="str"))
```

which results in

```
Two-electron many-body basis set (binary strings):
['000011' '000101' '001001' '010001' '100001' '000110' '001010' '010010'
'100010' '001100' '010100' '100100' '011000' '101000' '110000']
```

We remark that these strings must be read from right to left when converting into orbital occupations. The above basis set may be written as an array of integers whose rows correspond to many-body basis states and columns to single-particle orbitals with

```
print("Two-electron many-body basis set (array): ")
print(dvc.get_N_particle_subspace(2).get_bas_set(dtype="array"))
```

which results in

```
Two-electron many-body basis set (array):
[[1 1 0 0 0 0]
[1 0 1 0 0 0]
[1 0 0 1 0 0]
[1 0 0 0 1 0]
[1 0 0 0 0 1]
[0 1 1 0 0 0]
[0 1 0 1 0 0]
[0 1 0 0 1 0]
[0 1 0 0 0 1]
[0 0 1 1 0 0]
[0 0 1 0 1 0]
[0 0 1 0 0 1]
[0 0 0 1 1 0]
[0 0 0 1 0 1]
[0 0 0 0 1 1]]
```

This corresponds to all possible occupancy configurations of the \(n_\mathrm{states}\times n_\mathrm{degen}=6\) single-particle orbitals in which the total number of electrons is \(2\).

We can also display the many-body eigenvalues for a specific subspace, here the two-electron subspace, with

```
print("Two-electron many-body eigenvalues (eV): ")
print(dvc.get_N_particle_subspace(2).eigval/ct.e)
```

which yields

```
Two-electron many-body eigenvalues (eV):
[1.18586872 1.219271 1.219271 1.219271 1.25911222 1.29423484
1.29430848 1.29430848 1.29430848 1.32967697 1.35302216 1.35302216
1.35302216 1.38875253 1.44313977]
```

In addition, we can display, e.g., the first two-electron eigenstates with

```
print("Two-electron many-body ground state:")
print(dvc.get_N_particle_subspace(2).eigvec[0])
print("Two-electron many-body first excited state:")
print(dvc.get_N_particle_subspace(2).eigvec[1])
```

which produces

```
Two-electron many-body ground state:
[ 9.64208035e-01 -7.66382879e-16 1.25976213e-16 8.49194254e-16
1.52173170e-04 -2.56195602e-01 -4.67345620e-15 -3.46790780e-15
2.74378503e-13 -4.74138125e-02 1.35715888e-17 1.92644040e-46
-6.28745796e-17 4.48402977e-05 -4.91791061e-02]
Two-electron many-body first excited state:
[-2.05391260e-15 -7.48120396e-01 -6.51801647e-01 -2.59957417e-02
1.99563673e-15 2.77555756e-15 1.00144701e-06 1.42495485e-05
-4.93700270e-06 2.22044605e-16 1.04975162e-01 1.10441910e-02
6.04396188e-02 -1.15467531e-16 8.32667268e-16]
```

As can be seen from above, the two-body eigenstates are superpositions of the many-body basis states. We remark that because the two-body excited states are degenerate in this example, the two-body first excited states span a Hilbert space. The specific combination of coefficients in the first excited state displayed above is but one of the infinite potential linear combinations that exist within this Hilbert space. This specific outcome is determined by the eigenvalue solver and may thus differ between machines or QTCAD versions.

In the trivial case of the single-body subspace, the basis states are also many-body eigenstates. This can be verified with

```
print("One-electron many-body basis set:")
print(dvc.get_N_particle_subspace(1).get_bas_set(dtype="array"))
print("One-electron many-body eigenstates:")
print(dvc.get_N_particle_subspace(1).eigvec)
```

which produces

```
One-electron many-body basis set:
[[1 0 0 0 0 0]
[0 1 0 0 0 0]
[0 0 1 0 0 0]
[0 0 0 1 0 0]
[0 0 0 0 1 0]
[0 0 0 0 0 1]]
One-electron many-body eigenstates:
[[1. 0. 0. 0. 0. 0.]
[0. 1. 0. 0. 0. 0.]
[0. 0. 1. 0. 0. 0.]
[0. 0. 0. 1. 0. 0.]
[0. 0. 0. 0. 1. 0.]
[0. 0. 0. 0. 0. 1.]]
```

Also, remark that information about the many-body eigenstates and energies can
be equivalently obtained through the
`get_eig_state`

method,
which outputs an
`MbState`

object capturing all the
physical information about a many-body state. For example, we can display the
two-electron ground-state coefficients (of the expansion over the
many-body basis) and energy with

```
# Many-body state objects
print("Two-electron ground-state coefficients and energy (eV)")
print(dvc.get_N_particle_subspace(2).get_eig_state(0).coeff)
print(dvc.get_N_particle_subspace(2).get_eig_state(0).energy/ct.e)
```

which results in

```
Two-electron first ground-state coefficients and energy (eV)
[ 9.64208035e-01 -2.05391260e-15 6.86386230e-05 -1.64798730e-17
1.32048166e-01 -6.86386230e-05 -1.16356984e-29 -1.32048166e-01
1.92592994e-34 -1.86857050e-01 1.96568220e-16 1.23761073e-05
-1.23761073e-05 1.13590657e-30 -2.26685297e-02]
1.1858687241148644
```

### 13.5.2. Chemical potentials and Coulomb peak positions

We next output the chemical potentials and Coulomb peak positions for this device and gate bias configuration (see Chemical potentials and Coulomb peak positions).

This is done with

```
# Chemical potentials and Coulomb peak positions
print("Chemical potentials (eV)")
print(dvc.chem_potentials/ct.e)
print("Coulomb peak positions (V)")
print(dvc.coulomb_peak_pos)
```

which yields

```
Chemical potentials (eV)
[0.54798018 0.63788855 0.74948562]
Coulomb peak positions (V)
[1.09596035 1.27577709 1.49897124]
```

Because we have taken the lever arm \(\alpha\) to be \(0.5\), the position of the Coulomb peaks is twice the chemical potentials.

Note

The position of the Coulomb peaks is given relative to the reference gate potential. The reference gate potential is the potential applied at the boundary condition when solving Poisson’s equation.

For example, at vanishing source–drain bias, to see the Coulomb peak associated with the transition from \(0\) to \(1\) electron inside the potential well in the GAAFET channel, a gate bias of \(V_\mathrm{GS}+\mu(1)/e\alpha = 0.5\;\mathrm V+0.54798\;\mathrm V/0.5 = 1.59596\;\mathrm V\) would need to be applied for a lever arm \(\alpha=0.5\).

### 13.5.3. Many-body particle density

Finally, we may plot the particle densities associated with many-body eigenstates of interest as follows.

```
# Plot electron density
an.plot_slices(mesh, dvc.get_many_body_density(1,0),
title="One-electron ground state particle density (1/m^3)")
an.plot_slices(mesh, dvc.get_many_body_density(1,2),
title="One-electron second-excited state particle density (1/m^3)")
an.plot_slices(mesh, dvc.get_many_body_density(2,0),
title="Two-electron ground state particle density (1/m^3)")
an.plot_slices(mesh, dvc.get_many_body_density(2,2),
title="Two-electron second-excited state particle density (1/m^3)")
an.plot_slices(mesh, dvc.get_many_body_density(3,0),
title="Three-electron ground state particle density (1/m^3)")
an.plot_slices(mesh, dvc.get_many_body_density(3,5),
title="Three-electron fifth-excited state particle density (1/m^3)")
```

The resulting plots are shown below.

Note

The
`Subspace`

objects stored in the
`many_body_subspaces`

of a
`Device`

object account for
single-particle state degeneracy explicitly. As a result, in this example,
the particle density of the first single-electron “excited” state is the
same as for the single-electron ground state, simply because these states
are taken to be degenerate. Therefore, the first single-electron eigenstate
whose particle density differs from the ground state is the one with
eigenstate label \(2\).

## 13.6. Full code

```
__copyright__ = "Copyright 2024, Nanoacademic Technologies Inc."
from qtcad.device import constants as ct
from qtcad.device.mesh3d import Mesh
from qtcad.device import io
from qtcad.device import analysis as an
from qtcad.device import materials as mt
from qtcad.device import Device
from qtcad.device.schrodinger import Solver as SingleBodySolver
import pathlib
from qtcad.device.many_body import Solver as ManyBodySolver
from qtcad.device.many_body import SolverParams
# Define some device parameters
Vgs = 0.5 # 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.msh"
path_hdf5 = script_dir / "output/" / "gaafet.hdf5"
# Load the mesh
scaling = 1e-9
mesh = Mesh(scaling, str(path_mesh))
# Create device from mesh and set statistics.
# It is important to set statistics before creating boundaries since
# they determine the boundary values, Failing to do so will produce
# inconsistent results.
dvc = Device(mesh, conf_carriers="e")
dvc.statistics = "FD_approx" # Analytic approximation to Fermi-Dirac statistics
# Create regions first.
# Make sure each element is assigned to a region, otherwise default
# (silicon) parameters are used
dvc.new_region('channel_ox',mt.SiO2)
dvc.new_region('source_ox',mt.SiO2)
dvc.new_region('drain_ox',mt.SiO2)
dvc.new_region('channel', mt.Si)
dvc.new_region('source', mt.Si, pdoping=5e20*1e6, ndoping=0)
dvc.new_region('drain', mt.Si, pdoping=5e20*1e6, ndoping=0)
# Then create boundaries
dvc.new_ohmic_bnd('source_bnd')
dvc.new_ohmic_bnd('drain_bnd')
dvc.new_gate_bnd('gate_bnd',Vgs,Ew)
# Load the electic potential
dvc.set_potential(io.load(path_hdf5,"phi"))
# Solve the single-electron problem
slv = SingleBodySolver(dvc)
slv.solve()
dvc.print_energies()
# Solve the many-body problem
solver_params = SolverParams()
solver_params.num_states = 3
solver_params.n_degen = 2
solver_params.alpha = 0.5
solver_params.num_particles = [0,1,2,3]
slv = ManyBodySolver(dvc, solver_params = solver_params)
slv.solve()
# Many-body subspace properties
print("Many-body subspaces:", dvc.many_body_subspaces)
print("Number of subspaces:", len(dvc.many_body_subspaces))
print("Electron number in subspace 2:", dvc.many_body_subspaces[2].N)
print("Two-electron many-body basis set (integers): ")
print(dvc.get_N_particle_subspace(2).get_bas_set())
print("Two-electron many-body basis set (binary strings): ")
print(dvc.get_N_particle_subspace(2).get_bas_set(dtype="str"))
print("Two-electron many-body basis set (array): ")
print(dvc.get_N_particle_subspace(2).get_bas_set(dtype="array"))
print("Two-electron many-body eigenvalues (eV): ")
print(dvc.get_N_particle_subspace(2).eigval/ct.e)
print("Two-electron many-body ground state:")
print(dvc.get_N_particle_subspace(2).eigvec[0])
print("Two-electron many-body first excited state:")
print(dvc.get_N_particle_subspace(2).eigvec[1])
print("One-electron many-body basis set:")
print(dvc.get_N_particle_subspace(1).get_bas_set(dtype="array"))
print("One-electron many-body eigenstates:")
print(dvc.get_N_particle_subspace(1).eigvec)
# Many-body state objects
print("Two-electron ground-state coefficients and energy (eV)")
print(dvc.get_N_particle_subspace(2).get_eig_state(0).coeff)
print(dvc.get_N_particle_subspace(2).get_eig_state(0).energy/ct.e)
# Chemical potentials and Coulomb peak positions
print("Chemical potentials (eV)")
print(dvc.chem_potentials/ct.e)
print("Coulomb peak positions (V)")
print(dvc.coulomb_peak_pos)
# Plot electron density
an.plot_slices(mesh, dvc.get_many_body_density(1,0),
title="One-electron ground state particle density (1/m^3)")
an.plot_slices(mesh, dvc.get_many_body_density(1,2),
title="One-electron second-excited state particle density (1/m^3)")
an.plot_slices(mesh, dvc.get_many_body_density(2,0),
title="Two-electron ground state particle density (1/m^3)")
an.plot_slices(mesh, dvc.get_many_body_density(2,2),
title="Two-electron second-excited state particle density (1/m^3)")
an.plot_slices(mesh, dvc.get_many_body_density(3,0),
title="Three-electron ground state particle density (1/m^3)")
an.plot_slices(mesh, dvc.get_many_body_density(3,5),
title="Three-electron fifth-excited state particle density (1/m^3)")
```