4. Computing the chemical potential of an SET device

4.1. Requirements

4.1.1. Software components

• QTCAD

4.1.2. Python script

• qtcad/examples/practical_application/FDSOI/chemical_potential.py

4.1.3. References

• Mesh generation using QTCAD Builder

• Creating the FD-SOI Device

• Simulating electrostatics using the linear Poisson solver

• Transport in the Coulomb blockade regime

4.2. Briefing

In this section of the tutorial, we will demonstrate how to compute the chemical potential, \(\mu\), of the single-electron transistor (SET) as a function of the gate voltage applied to the associated gate, in this case, "plunger_gate_1_bnd".

In principle, to properly compute the chemical potential, one would need to solve the full many-body Schrödinger equation. However, this is computationally very expensive and not feasible for large systems. Instead, we will use a simplified approach to account for many-body effects by solving the Schrödinger–Poisson equation, which self-consistently couples the single-particle Schrödinger equation with the Poisson equation. This approach captures the essential physics of electron-electron interactions in a mean-field manner, making it computationally more tractable.

Note

While the Schrödinger–Poisson approach (outlined in Schrödinger–Poisson solver) incorporates many-body effects in a mean-field manner, it has important limitations. In particular, it does not account for exchange and correlation contributions to the energy. Consequently, results obtained with this method should be interpreted with caution in regimes where strong correlations are expected to be significant.

Furthermore, unlike approaches such as Hartree–Fock [BSA00], the Schrödinger–Poisson method described in Schrödinger–Poisson solver does not eliminate electron self-interaction, which can introduce inaccuracies in the computed chemical potentials. We expect the relative impact of self-interaction errors to diminish as the number of electrons increases, making this approach more reliable for systems with a large number of electrons.

4.3. Header

We start by importing the necessary Python modules.

import pathlib
import numpy as np
from numpy.typing import NDArray

from qtcad.device import constants as ct
from qtcad.device import materials as mt
from qtcad.device import Device, SubDevice
from qtcad.device.schrodinger_poisson import Solver
from qtcad.device.schrodinger_poisson import SolverParams
from qtcad.device.schrodinger import SolverParams as SchrodSolverParams
from qtcad.device.poisson import SolverParams as PoissonSolverParams
from qtcad.device import io

from double_dot_fdsoi import save_slice, state_probability_density

These modules should be familiar from previous tutorials. The state_probability_density helper will be used later to convert each valley-resolved eigenstate into the corresponding probability density before saving it for visualization.

We also define the path to the current directory

# Script directory
script_dir = pathlib.Path(__file__).parent.resolve()

and the degeneracy factor for silicon which specifies the degeneracy of each quantum state in our simulations.

# Degeneracy factor for silicon
# Start from the default confined-electron degeneracy and remove the
# valley factor, since valley splitting is treated explicitly.
g = int(mt.Si.gqc / 2)

Here we assume that the growth axis of the device (\(z\)) is oriented along the \([001]\) crystallographic direction. Because the device is more strongly confined in the \(z\) direction than in \(x\) or \(y\), the valley states associated with the \(z\) axis respond differently to confinement. In silicon, electrons in the \(z\) valleys have a larger effective mass along the \(z\) axis than electrons in the \(x\) and \(y\) valleys. Since the quantization energy in a confined direction decreases as the effective mass increases, states with a heavier effective mass end up at lower energies for the same confinement geometry. Consequently, the \(z\)-valley states lie lowest in energy, leading to a twofold valley degeneracy associated with the two equivalent \(z\) valleys. If the valleys were not accounted for explicitly, each confined level would therefore carry a twofold spin degeneracy and a twofold valley degeneracy.

In QTCAD, this default confined-electron degeneracy is stored in the silicon Material attribute mt.Si.gqc. For silicon under quantum confinement, its value is 4. In this practical application, however, the FD-SOI device is constructed with a nonzero valley splitting through set_valley_splitting (see Creating the device). This means that the two low-lying valley states are accounted for explicitly in the single-particle spectrum rather than through the degeneracy factor. We therefore do not want to keep the valley contribution inside g. Instead, we start from the default value mt.Si.gqc = 4, divide it by two, and use g=2 so that g accounts only for the residual spin degeneracy.

Note

In reality there is typically a coupling and an associated splitting between the two z valleys due to device-specific features such as interfaces, confinement, electric fields, and disorder. In a realistic device, the resulting valley splitting would therefore depend on the details of the structure being simulated. In this tutorial, however, we simply choose a fixed representative value of \(0.5\,\mathrm{meV}\) and include it explicitly in the device definition (see Creating the device). This is a typical order of magnitude for silicon MOS/SiO₂ quantum dots [YRR+13], so the residual degeneracy factor g accounts only for spin. It is also possible to compute the valley splitting direvtly via tight binding using QTCAD Atoms (see, e.g., Multiscale simulations of a SiGe–Si–SiGe quantum dot—Part 2: Tight-binding model for valley splitting).

4.4. Solver parameters

The main goal of this tutorial series is to simulate the behavior of an SET using the Schrödinger–Poisson solver. Because most of the steps of the tutorial involve solving the Schrödinger–Poisson equations, we define here the default SolverParams parameters that will be used in the calculations.

# Solver parameters
# (Poisson)
p_params                            = PoissonSolverParams()
p_params.tol                        = 1e-8
# (Schrödinger)
schrod_params                       = SchrodSolverParams()
schrod_params.tol                   = 1e-8
schrod_params.num_states            = 10
# (Schrödinger-Poisson)
sp_params                           = SolverParams()
sp_params.tol                       = 1e-5
sp_params.maxiter                   = 2000
sp_params.bound_state_charges_only  = True
sp_params.sc_method                 = 'underrelax'
sp_params.initialization            = False
sp_params.mixing_algo               = "adaptive_linear"
sp_params.mixing_param              = 0.1
sp_params.poisson_solver_params     = p_params
sp_params.schrod_solver_params      = schrod_params

We start by setting up the parameters for the Poisson solver (PoissonSolverParams) and the Schrödinger solver (SchrodSolverParams). Then we create an instance of the Schrödinger–Poisson SolverParams class that contains the parameters for the Schrödinger–Poisson Solver. These parameters include:

• tol — Convergence tolerance for the self-consistent loop. Here the tolerance is set to 1e-5 eV, meaning that the solver will consider the solution converged when the maximal relative change in the potential between successive iterations is less than this value. The reason we require such a strict tolerance is that we will later compute many-body energies, chemical potentials, and addition energies which are expressed as sums or difference of single-particle eigenenergies. The typical values we will be interested in are on the order of one to ten meV. Because we will be interested in energies of \(N \sim 12\) electrons, the error on the quantities of interest can be on the order of 10 times larger than the we error on the single-particle eigenenergies. Therefore, we need an accuracy of each single-particle eigenenergy to be significantly better than 1 meV to ensure that the computed energies are accurate.

• maxiter — Maximum number of iterations for the self-consistent loop.

• bound_state_charges_only — Whether or not to consider only bound-state charges. Here it is set to True, specifying that the only charges contributing to the Poisson equation are those from bound states (electrons occupying eigenstates of the Schrödinger equation). If False, classical charges from conduction and valence bands would also be included. When bound_state_charges_only is True, the sc_method parameter (see below) is restricted to 'underrelax'.

• sc_method — Self-consistent method to be used. Here it is set to 'underrelax'. The other option is 'p-c' for predictor–corrector. While predictor–corrector has been shown to be more efficient when simulating systems containing electron gases [TGPR97] originating from the filling of parabolic bands, here, all such charge is neglected (see description of bound_state_charges_only parameter, above).

• initialization — Whether to generate an initial guess for the electric potential phi. Here, it is set to False since we will provide our own initial guess, either by first solving the linear Poisson equation (see Simulating electrostatics using the linear Poisson solver) or by using a previously computed solution.

• mixing_algo — Mixing algorithm to be used. Here, we use the "adaptive_linear" algorithm which adaptively adjusts the mixing parameter for a linear-based mixing scheme, based on the convergence behavior.

• mixing_param — Mixing parameter. In the case of the "adaptive_linear" algorithm, this parameter serves as the initial mixing parameter.

• poisson_solver_params — The Poisson solver parameters defined earlier.

• schrod_solver_params — The Schrödinger solver parameters defined earlier.

4.5. Function to compute the chemical potential

Next, we define the main content of the corresponding Python file: a function, compute_chemical_potential, which computes the chemical potential as a function of the gate voltage applied to the SET gate, "plunger_gate_1_bnd".

4.5.1. Function description

The first step is to define the function signature and docstrings.

def compute_chemical_potential(
    V_SET: NDArray,
    Nmin: int,
    Nmax: int,
    d: Device,
    set_device: SubDevice,
    out_path: pathlib.Path=script_dir
) -> NDArray:
    """Compute the chemical potential for different number of particles and
    different voltage configurations. This function also saves the potential
    to HDF5 files, eigenstate probability densities to VTU files, and the
    energies and population factors to text files. All output files are saved
    in the specified output directory.

    Args:
        V_SET: Array of SET gate voltages.
        Nmin: Minimum number of electrons.
        Nmax: Maximum number of electrons.
        d: Device object representing the entire FD-SOI device.
        set_device: SubDevice representing the SET region of the FD-SOI device.
        out_path: Output directory path.

    Returns:
        2D array of chemical potentials for each voltage configuration.
        The first index runs over the applied voltages to the SET gate, and the
        second index runs over the number of particles.

    Note:
        The chemical potential is calculated as the difference in many-body
        ground-state energies when adding an electron to the system. The
        chemical potentials are computed for the electron numbers in the range
        [Nmin+1, Nmax].
    """

As is explained in the docstring, the function takes as input an array of gate voltages, the minimum and maximum number of electrons to consider, the Device object representing the entire FD-SOI device, the SubDevice object representing the SET region, and the output directory path. It returns a 2D array of chemical potentials for each voltage configuration and number of particles.

As mentioned above, in this tutorial, we will compute the chemical potential when the number of electrons in the SET is around \(N=12\), for different values of the "plunger_gate_1_bnd" voltage. The chemical potential \(\mu(N)\) is defined as the energy required to add the \(N\)-th electron to the system, i.e., \(\mu(N) = U(N) - U(N-1)\), where \(U(N)\) is the ground-state energy of the system with \(N\) electrons [HKP+07]. Therefore, for each value of the gate voltage, we will compute the total energy of the system by converging a Schrödinger–Poisson calculation for \(N-1 = 12\), \(N = 13\), and \(N+1=14\) electrons to obtain the chemical potentials \(\mu(N=13)\) and \(\mu(N=14)\). Users are free to modify the range of electron numbers and gate voltages as needed.

4.5.2. Preliminaries

Inside the function, the first step is to define the paths to the output files that will be used to save the results of the simulations

# output files
# (results)
phi_out_file         = str(out_path / "phi_vset{vset:.2f}_N{n}.hdf5")
mu_file              = out_path / "mu.txt"
pop_file             = out_path / "population.txt"
E_file               = out_path / "energies.txt"

# (visualization)
wf_vtu_file          = str(out_path / "wf_vset{vset:.2f}_N{n}.vtu")
rho_slice_file       = str(out_path / "rho_vset{vset:.2f}_N{n}_slice.png")
wf_slice_file        = str(out_path / "wf_vset{vset:.2f}_N{n}_slice.png")

Then, the outputs of intermediary steps and of the entire calculation are initialized.

# Range of number of particles
N = range(Nmin, Nmax + 1)
# Many-body energies
mb_E = np.zeros((len(V_SET), len(N)))
# Chemical potential
Mu = np.zeros((len(V_SET), len(N) - 1))

We start by defining the range of number of particles to consider, from Nmin to Nmax. We then initialize a 2D array, mb_E, to store the many-body ground-state energies for each value of the gate voltage and number of particles. Finally, we initialize a 2D array, Mu, to store the chemical potentials for each value of the gate voltage and number of particles.

4.5.3. Solving the Schrödinger–Poisson equation

Next, we loop over the values of the gate voltage, print some information about which calculation is currently being performed, and set the potential on the SET gate, "plunger_gate_1_bnd", accordingly

for j, vset in enumerate(V_SET):
    # Print info
    print("\n")
    print("==================================================")
    print(f"Setting V_SET = {vset} V")
    print("==================================================")
    print("\n")

    # Set the potential in the SET region
    d.set_applied_potential("plunger_gate_1_bnd", vset)

For each of these applied voltages, we then loop over the number of particles and perform a Schrödinger–Poisson calculation to obtain the ground-state energy of the system

for i, n in enumerate(N):
    # Print info
    print("--------------------------------------------------")
    print(f"Solving for N = {n} electrons...")
    print("--------------------------------------------------")
    print("\n")

    # Solve Schrödinger-Poisson in the SET region
    sp_solver = Solver(d, subdevice=set_device, solver_params=sp_params)
    sp_solver.solve(N=n, g=g)
    print("Energies and population factors of SET SubDevice:")
    set_device.print_energies()

The argument N passed to solve specifies the total number of electrons included in the Schrödinger–Poisson calculation. To enforce this electron number, the solver fills the single-particle energy levels according to the Fermi–Dirac distribution. The Fermi level is adjusted self-consistently so that the total occupation equals the requested number of electrons N.

You may also provide the degeneracy factor g to account for any residual degeneracy of the single-particle levels. This factor multiplies the number of electrons allowed in each energy level (i.e., each level can accommodate g electrons instead of one). In this practical application, the residual degeneracy is due to spin since valleys are accounted for explicitly. The g argument is optional. If it is not provided, the Solver automatically determines an effective degeneracy factor by averaging the Device attribute gqc over the mesh, weighted by the ground-state electron density.

4.5.3.1. Valley-Resolved Indexing

Because the single-particle spectrum is printed explicitly at this stage, it is helpful to introduce the corresponding valley-resolved notation here. We denote by \(\psi_{i,\nu}(\mathbf{r})\) the \(\nu\)-th valley component of the \(i\)-th valley-resolved single-particle state, where \(i\) labels the entries of energies and \(\nu\) labels the valley index of eigenfunctions.

In the present electron-only FD-SOI example without explicit spin, the eigenfunction array has shape (N_nodes, N_levels, 2). Therefore, set_device.eigenfunctions[:, i, 0] and set_device.eigenfunctions[:, i, 1] are the two explicit valley components associated with the level set_device.energies[i]. In this simplified example, each eigenstate has weight in only one valley component. For instance, eigenstate 0 has a nonzero contribution only in component 0, so set_device.eigenfunctions[:, 0, 0] is nonzero while set_device.eigenfunctions[:, 0, 1] vanishes. Likewise, eigenstate 1 has a nonzero contribution only in component 1, so set_device.eigenfunctions[:, 1, 1] is nonzero while set_device.eigenfunctions[:, 1, 0] vanishes. The corresponding scalar probability density of state \(i\) is then \(\sum_\nu |\psi_{i,\nu}(\mathbf{r})|^2\).

Because the FD-SOI device was constructed with an explicit valley splitting, the energy list printed by print_energies already shows each original orbital as a pair of consecutive valley-partner levels. In other words, when two adjacent printed levels belong to the same orbital, their separation is the prescribed valley splitting \(\Delta_v\).

4.5.4. Saving results

After solving the Schrödinger and Poisson equations self-consistently, we save the potential to an .hdf5 file for later use and write the probability densities associated with the Schrödinger eigenstates to a .vtu file for visualization. We also save slices of the ground-state probability density and total electron density for quick inspection.

# Save phi, eigenstate densities, and charge density
io.save(phi_out_file.format(vset=vset, n=n), {"phi": d.phi})
wf_dict = {
    f"WF_{i}": state_probability_density(set_device, i)
    for i in range(set_device.eigenfunctions.shape[1])
}
io.save(wf_vtu_file.format(vset=vset, n=n), wf_dict, set_device.mesh)

save_slice(
    set_device,
    state_probability_density(set_device, 0),
    wf_slice_file.format(vset=vset, n=n),
    label="Ground-state probability density [1 / m$^3$]",
)
save_slice(
    d,
    d.rho / ct.e,
    rho_slice_file.format(vset=vset, n=n),
    label="$\\rho$ [$e$ / m$^3$]"
)

Because the valley splitting is accounted for explicitly in this FD-SOI workflow, the Schrödinger eigenfunctions carry an additional trailing axis that indexes the valley components. When the solver applies a valley splitting \(\Delta_v\), each original single-particle level \(E\) is replaced by two consecutive levels, \(E - \Delta_v/2\) and \(E + \Delta_v/2\), and the corresponding wavefunction is copied into the two valley sectors. In other words, state 1 is the valley partner of state 0, state 3 is the valley partner of state 2, and so on. The state_probability_density helper acts on one state index at a time and then computes \(\sum_\nu |\psi_{i,\nu}(\mathbf{r})|^2\) from the band/valley axis, yielding the scalar probability density associated with that single valley-resolved state. This is the quantity we save in the WF_i datasets and use for the ground-state slice.

Next, we save the population factors and the single-particle energies to output files.

# Save energies and population factors
with open(pop_file, "a") as f:
    out = [vset, n] + list(set_device.population_factors)
    np.savetxt(f, np.array(out)[np.newaxis, :], fmt="%.6f")
with open(E_file, "a") as f:
    out = [vset, n] + list(set_device.energies / ct.e)
    np.savetxt(f, np.array(out)[np.newaxis, :], fmt="%.8f")

The population factors indicate the occupation of each single-particle energy level in the SET region. Each of the factors is a real number between 0 and 1, where a value of 1 indicates that the corresponding energy level is fully occupied, while a value of 0 indicates that the level is empty. It is important to note that the population factors do not include the degeneracy. In other words, each population factor should be multiplied by the degeneracy of the corresponding energy level to obtain the actual number of electrons occupying that level. In the present FD-SOI example, the valley splitting is included explicitly in the device, so the only residual degeneracy is the twofold spin degeneracy. Therefore each population factor should be multiplied by g=2 to obtain the actual occupancy of the corresponding level.

4.5.5. Computing chemical potentials

The final step is to compute the many-body ground-state energy and the chemical potential. The many-body ground-state energy is computed as the sum of the single-particle energies weighted by their respective population factors.

# Compute many-body ground-state energy
mb_E[j, i] = np.sum(g * set_device.population_factors * set_device.energies)
print("  Many-body ground-state energy (eV): "f"{mb_E[j, i] / ct.e:.10f}")

The chemical potential is then computed as the difference in many-body ground-state energies when adding an electron to the system [HKP+07].

        if i > 0:
            # Compute chemical potential
            Mu[j, i-1] = (mb_E[j, i] - mb_E[j, i-1])
            print(f"  Chemical potentials (eV): {Mu[j, i-1] / ct.e}")
            with open(mu_file, "a") as f:
                np.savetxt(f, np.array([[vset, n, Mu[j, i-1] / ct.e]]), fmt="%.10f")

return Mu

4.6. Full code

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

import pathlib
import numpy as np
from numpy.typing import NDArray
from qtcad.device import constants as ct
from qtcad.device import materials as mt
from qtcad.device import Device, SubDevice
from qtcad.device.schrodinger_poisson import Solver 
from qtcad.device.schrodinger_poisson import SolverParams
from qtcad.device.schrodinger import SolverParams as SchrodSolverParams
from qtcad.device.poisson import SolverParams as PoissonSolverParams
from qtcad.device import io
from double_dot_fdsoi import save_slice, state_probability_density

# Script directory
script_dir = pathlib.Path(__file__).parent.resolve()

# Degeneracy factor for silicon
# Start from the default confined-electron degeneracy and remove the
# valley factor, since valley splitting is treated explicitly.
g = int(mt.Si.gqc / 2)

# Solver parameters
# (Poisson)
p_params                            = PoissonSolverParams()
p_params.tol                        = 1e-8
# (Schrödinger)
schrod_params                       = SchrodSolverParams()
schrod_params.tol                   = 1e-8
schrod_params.num_states            = 10
# (Schrödinger-Poisson)
sp_params                           = SolverParams()
sp_params.tol                       = 1e-5
sp_params.maxiter                   = 2000
sp_params.bound_state_charges_only  = True
sp_params.sc_method                 = 'underrelax'
sp_params.initialization            = False
sp_params.mixing_algo               = "adaptive_linear"
sp_params.mixing_param              = 0.1
sp_params.poisson_solver_params     = p_params
sp_params.schrod_solver_params      = schrod_params

def compute_chemical_potential(
    V_SET: NDArray, 
    Nmin: int, 
    Nmax: int, 
    d: Device, 
    set_device: SubDevice, 
    out_path: pathlib.Path=script_dir
) -> NDArray:
    """Compute the chemical potential for different number of particles and
    different voltage configurations. This function also saves the potential
    to HDF5 files, eigenstate probability densities to VTU files, and the
    energies and population factors to text files. All output files are saved
    in the specified output directory.

    Args:
        V_SET: Array of SET gate voltages.
        Nmin: Minimum number of electrons.
        Nmax: Maximum number of electrons.
        d: Device object representing the entire FD-SOI device.
        set_device: SubDevice representing the SET region of the FD-SOI device.
        out_path: Output directory path.

    Returns:
        2D array of chemical potentials for each voltage configuration. 
        The first index runs over the applied voltages to the SET gate, and the 
        second index runs over the number of particles.

    Note:
        The chemical potential is calculated as the difference in many-body 
        ground-state energies when adding an electron to the system. The 
        chemical potentials are computed for the electron numbers in the range
        [Nmin+1, Nmax].
    """

    # output files
    # (results)
    phi_out_file         = str(out_path / "phi_vset{vset:.2f}_N{n}.hdf5")
    mu_file              = out_path / "mu.txt"
    pop_file             = out_path / "population.txt"
    E_file               = out_path / "energies.txt"

    # (visualization)
    wf_vtu_file          = str(out_path / "wf_vset{vset:.2f}_N{n}.vtu")
    rho_slice_file       = str(out_path / "rho_vset{vset:.2f}_N{n}_slice.png")
    wf_slice_file        = str(out_path / "wf_vset{vset:.2f}_N{n}_slice.png")

    # Range of number of particles
    N = range(Nmin, Nmax + 1)
    # Many-body energies
    mb_E = np.zeros((len(V_SET), len(N)))
    # Chemical potential
    Mu = np.zeros((len(V_SET), len(N) - 1))

    for j, vset in enumerate(V_SET):
        # Print info
        print("\n")
        print("==================================================")
        print(f"Setting V_SET = {vset} V")
        print("==================================================")
        print("\n")

        # Set the potential in the SET region
        d.set_applied_potential("plunger_gate_1_bnd", vset)

        for i, n in enumerate(N):
            # Print info
            print("--------------------------------------------------")
            print(f"Solving for N = {n} electrons...")
            print("--------------------------------------------------")
            print("\n")

            # Solve Schrödinger-Poisson in the SET region
            sp_solver = Solver(d, subdevice=set_device, solver_params=sp_params)
            sp_solver.solve(N=n, g=g)
            print("Energies and population factors of SET SubDevice:")
            set_device.print_energies()

            # Save phi, eigenstate densities, and charge density
            io.save(phi_out_file.format(vset=vset, n=n), {"phi": d.phi})
            wf_dict = {
                f"WF_{i}": state_probability_density(set_device, i)
                for i in range(set_device.eigenfunctions.shape[1])
            }
            io.save(wf_vtu_file.format(vset=vset, n=n), wf_dict, set_device.mesh)

            save_slice(
                set_device,
                state_probability_density(set_device, 0),
                wf_slice_file.format(vset=vset, n=n),
                label="Ground-state probability density [1 / m$^3$]",
            )
            save_slice(
                d, 
                d.rho / ct.e, 
                rho_slice_file.format(vset=vset, n=n), 
                label="$\\rho$ [$e$ / m$^3$]"
            )

            # Save energies and population factors
            with open(pop_file, "a") as f:
                out = [vset, n] + list(set_device.population_factors)
                np.savetxt(f, np.array(out)[np.newaxis, :], fmt="%.6f")
            with open(E_file, "a") as f:
                out = [vset, n] + list(set_device.energies / ct.e)
                np.savetxt(f, np.array(out)[np.newaxis, :], fmt="%.8f")

            # Compute many-body ground-state energy
            mb_E[j, i] = np.sum(g * set_device.population_factors * set_device.energies)
            print("  Many-body ground-state energy (eV): "f"{mb_E[j, i] / ct.e:.10f}")

            if i > 0:
                # Compute chemical potential
                Mu[j, i-1] = (mb_E[j, i] - mb_E[j, i-1])
                print(f"  Chemical potentials (eV): {Mu[j, i-1] / ct.e}")
                with open(mu_file, "a") as f:
                    np.savetxt(f, np.array([[vset, n, Mu[j, i-1] / ct.e]]), fmt="%.10f")

    return Mu