9. Schrödinger simulation of a quantum well

9.1. Requirements

9.1.1. Software components

  • QTCAD

  • Gmsh

9.1.2. Geometry file

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

9.1.3. Python script

  • qtcad/examples/tutorials/quantum_well_holes.py

9.1.4. References

9.2. Briefing

In this tutorial, we explain how QTCAD’s Schrödinger solver can be used to simulate structures in lower than three dimensions. An important example is a quantum well where we have confinement along a single direction (\(z\) in our case) and translational invariance along the other two directions (\(x\) and \(y\)). For electrons under the effective mass approximation (where the effective mass is isotropic), this calculation is straightforward. For holes, when there is coupling between bands, it is more complicated and requires projecting the Hamiltonian onto a lower dimension accounting for the wavefunction along the confined dimension. In this tutorial we will calculate the band structure for a hole described by a six-band model [see [For93] for details], bound to a quantum well (1d confinement). We note that the six-band model is only available (for now) for GaAs.

We assume that the user is familiar with basic Python. We also assumed that the user has access to the Gmsh geometry file describing the quantum well. Geometry files for all the QTCAD tutorials are stored in the folder examples/tutorials/meshes.

The full code may be found at the bottom of this page, or in qtcad/examples/tutorials/quantum_well_holes.py.

Do not hesitate to look into the API source code for detailed descriptions of the inputs and outputs of each class, method, or function used in the tutorials. The source code also contains documentation for many other useful features not described here.

9.3. Mesh generation

The mesh is generated from the geometry file qtcad/examples/tutorials/meshes/qw_1d.geo by running Gmsh. This can be done in the terminal with

gmsh qw_1d.geo

This will produce a .msh2 file in the same directory as the geometry file (quantum_well.msh2). See Getting Started for more information about interfacing with Gmsh. Currently, QTCAD only accepts 1D Gmsh meshes along the \(x\) axis.

9.4. Setting up the device

9.4.2. Defining the mesh

The mesh is defined by loading the mesh file produced by Gmsh and setting units of length.

# Paths
path = pathlib.Path(__file__).parent.resolve()
path_mesh = str(path/"meshes"/"quantum_well.msh2")

# Load the mesh
scale = 1e-9
mesh = Mesh(scale, path_mesh)

For more details, see e.g. Poisson and Schrödinger simulation of a GAAFET. Because we have confinement only along 1 direction, we only need a 1d mesh to solve our problem. In the other directions, the system is assumed to be translationally invariant. Thus, the eigenstates are plane waves along these translationally invariant directions and do not need to be solved using QTCAD. We focus here only on the confined part of the states (envelope functions along \(z\)).

9.4.3. Creating the device

The device to simulate is then initialized on the mesh defined above:

# Create the device
dvc = Device(mesh, conf_carriers='h')
# In this example we will consider a 6 valence band model
dvc.hole_kp_model = "luttinger_kohn_foreman_6band"

# Set the material to be GaAs
mat = mt.GaAs
dvc.new_region("well", mat)

# Then create insulator boundaries
dvc.new_insulator("left_barrier")
dvc.new_insulator("right_barrier")

The first line calls the constructor of the device class with the mesh object as an argument. We also specify that we wish to consider holes confined to the quantum well with the optional argument conf_carriers. The second line of code indicates which \(\mathbf{k} \cdot \mathbf{p}\) model we will consider for the confined holes. In this case we use the 6 band model [For93]. We could also consider the four band model by using dvc.hole_kp_model = "luttinger_kohn_foreman". See k.p module for more details. After the model is specified, we specify the material that hosts our well, GaAs and give boundary conditions. The new_insulator method is used to enforce the boundary conditions of the Schrödinger solver. When solving the Schrödinger eqaution, the solver will look for solutions that vanish at the 'left_barrier' and the 'right_barrier' of the mesh. Since there is no potential defined over the device, we are essential modelling a particle-in-a-box-like quantum well.

9.5. Solving the Schrödinger solver

9.5.1. Setting up the solver

The parameters of our Schrödinger solver are specified through an instance of a SolverParams object. The num_states attribute indicates we wish to solve for 20 holes states in ascending order, starting from the ground state. We also specify ti_directions which is an attribute that indicates the directions along which we have translational invariance. Here we specify the 'x' and 'y' directions in a list (the device dvc describes the problem along 'z'). We note that while 1D meshes must be defined along the 'x' axis using Gmsh, the ti_directions attribute allows users to specify the confinement (which can be different than 'x') and translationally-invariant directions. In this specific example, the Gmsh mesh defined along the 'x' is reinterpreted by QTCAD to represent a mesh defined along the 'z' axis.

# Solve Schrodinger's equation
# Parameters
params = SolverParams()
params.num_states = 20
params.ti_directions=['x', 'y']
# Solver
s = Solver(dvc, solver_params=params)

9.5.2. Solving the equation

Next, we use the solve method to solve the Schrödinger equation. This method will store the eigenfunctions and eigenenergies in attributes of dvc. As mentioned previously, we are assuming that we have translational invariance along \(x\) and \(y\). Thus, the momentum along each of these two directions (or equivalently, the wavevector) is a good quantum number. Therefore, the Schrödinger equation can, in principle, be solved for any in-plane wavevector \(\mathbf{k} = (k_x, k_y)\). By default, \(\mathbf{k} = (0, 0)\).

# Energies at (0, 0)
s.solve()
print("Energies at (0,0)")
dvc.print_energies()

The print_energies method prints the eigenenergies. The convention used in QTCAD is one where the holes have positive mass and charge (see QTCAD convention for holes for more details).

The solve method can also be used to solve the Schrödinger equation for a different in-plane wavevector. This is done using the keyword argument k.

# Energies at (1, 1)/aB
s.solve(k=np.array([1, 1])/ct.aB, mixing=True)
print("Energies at (1, 1)/aB")
dvc.print_energies()

Here we also solve the Schrödinger equation but for an in-plane wavevector \(\mathbf{k} = a_B^{-1}(1, 1)\), where \(a_B\) is the Bohr radius.

9.5.3. Band Structure

Given that we can solve the Schrödinger equation for different values of \(\mathbf{k}\), we can also compute a band structure for our quantum well.

# Band structure
K = np.linspace(0, 0.3/ct.aB, 30)
basis_vec = np.array([1, 0])

# With mixing
band_structure = s.band_structure(K, basis_vec)

We first specify an array, K, and a basis vector, basis_vec. Together these 2 quantities specify a series of \(\mathbf{k} = k\mathbf{v}\), where \(k\) is an element of K, and \(\mathbf{v}\) represents basis_vec. We pass these two objects into the band_structure method of the solver s. The output is a 2D array that contains the energies of the hole states for every \(\mathbf{k} = k\mathbf{v}\) specified. This output is \(30 \times 20\) in this case: we solve for 20 eigenstates for each of the 30 vectors \(\mathbf{k} = k\mathbf{v}\) specified.

The band_structure can be plotted using basic python plotting functionality

fig = plt.figure(figsize=(7,4))
ax1 = fig.add_subplot(1,1,1)
ax1.plot(K * ct.aB, band_structure/ct.e, '-')
ax1.set_xlabel("$k_x$ ($a_B^{-1}$)")
ax1.set_ylabel("Energy (eV)")
ax1.set_ylim([0, 5])
ax1.set_title("Band structure including band mixing")
plt.show()

Resulting in the following plot:

Band Structure

We see from this plot that the bands have positive curvature, consistent with the QTCAD convention for holes.

We can also compute the band structure in the absence of band mixing, i.e. by ignoring all coupling (off diagonal elements) between the valence bands of the model.

# Without mixing
# Compute band structure using QTCAD
band_structure = s.band_structure(K, basis_vec, mixing=False)

This is achieved using the mixing keyword argument of the band_structure method and setting it to False. Before plotting this band structure, we also compute it analytically. This is a straightforward calculation since we are computing particle in a box eigenstates for the confinement along \(z\) and free particle states along \(x\) and \(y\). Since the potential is separable, we can write the total energy as the sum of the energies which solve each of the 1d problems (along \(x\), \(y\), and \(z\) separately).

# Analytic solutions
# Particle in a box energies
def particle_in_box_E(m, L, n):
    return n**2*np.pi**2*ct.hbar**2 / ( 2 * m * L**2 )

g1, g2, g3, Delta = tuple(mat.hole_kp_params)
# Out of plane effective masses
mHH = ct.me / (g1-2*g2) # HH effective mass GaAs
mLH = ct.me / (g1+2*g2) # LH effective mass
mSO = ct.me / g1 # SO effective mass
# In plane effective masses
mHHip = ct.me / (g1 + g2) # HH in-plane effective mass GaAs
mLHip = ct.me / (g1 - g2) # LH in-plane effective mass GaAs
mSOip = ct.me / (g1) # SO in-plane effective mass GaAs

# Analytic energies
PIB_HH = np.array([particle_in_box_E(mHH, 2e-9, n) for n in range(1,5)])
free_HH = ct.hbar**2 / 2 / mHHip * K**2
E_HH = (PIB_HH[np.newaxis, :] + free_HH[:, np.newaxis])

PIB_LH = np.array([particle_in_box_E(mLH, 2e-9, n) for n in range(1,5)])
free_LH = ct.hbar**2 / 2 / mLHip * K**2
E_LH = (PIB_LH[np.newaxis, :] + free_LH[:, np.newaxis])

PIB_SO = np.array([particle_in_box_E(mSO, 2e-9, n) for n in range(1,5)])
free_SO = ct.hbar**2 / 2 / mSOip * K**2
E_SO = (PIB_SO[np.newaxis, :] + free_SO[:, np.newaxis]) + Delta

Here we have defined a function that computes the particle-in-a-box eigenstates. We then get the Luttinger parameters from the material file (in this case mat=mt.GaAs). After computing the in-plane and out-of-plane effective masses for each band, we sum the particle-in-a-box energies and the free electron energies for each band. (We also include the spin-orbit splitting for the split-off band).

Finally, we plot our results.

# Plot
fig = plt.figure(figsize=(7,4))
ax1 = fig.add_subplot(1,1,1)
ax1.plot(K * ct.aB, band_structure/ct.e, '-')
ax1.plot(K * ct.aB, E_HH/ct.e, '--b') # Plot analytic HH energies in blue
ax1.plot(K * ct.aB, E_LH/ct.e, '--r') # Plot analytic LH energies in blue
ax1.plot(K * ct.aB, E_SO/ct.e, '--g') # Plot analytic SO energies in blue
ax1.set_xlabel("$k_x$ ($a_B^{-1}$)")
ax1.set_ylabel("Energy (eV)")
ax1.set_ylim([0, 5])
ax1.set_title("Band structure in the absence of band mixing")
plt.show()

We find that the QTCAD solution (solid lines) matches the analytic solution (dashed lines):

Band Structure without mixing

9.6. Full code

__copyright__ = "Copyright 2023, Nanoacademic Technologies Inc."

from device import constants as ct
from device import materials as mt
from device.mesh1d import Mesh
from device.device import Device
from device.schrodinger import Solver, SolverParams
import numpy as np
from matplotlib import pyplot as plt
import pathlib

# Paths
path = pathlib.Path(__file__).parent.resolve()
path_mesh = str(path/"meshes"/"quantum_well.msh2")

# Load the mesh
scale = 1e-9
mesh = Mesh(scale, path_mesh)

# Create the device
dvc = Device(mesh, conf_carriers='h')
# In this example we will consider a 6 valence band model
dvc.hole_kp_model = "luttinger_kohn_foreman_6band"

# Set the material to be GaAs
mat = mt.GaAs
dvc.new_region("well", mat)

# Then create insulator boundaries
dvc.new_insulator("left_barrier")
dvc.new_insulator("right_barrier")

# Solve Schrodinger's equation
# Parameters
params = SolverParams()
params.num_states = 20
params.ti_directions=['x', 'y']
# Solver
s = Solver(dvc, solver_params=params)

# Energies at (0, 0)
s.solve()
print("Energies at (0,0)")
dvc.print_energies()


# Energies at (1, 1)/aB
s.solve(k=np.array([1, 1])/ct.aB, mixing=True)
print("Energies at (1, 1)/aB")
dvc.print_energies()

# Band structure
K = np.linspace(0, 0.3/ct.aB, 30)
basis_vec = np.array([1, 0])

# With mixing
band_structure = s.band_structure(K, basis_vec)

fig = plt.figure(figsize=(7,4))
ax1 = fig.add_subplot(1,1,1)
ax1.plot(K * ct.aB, band_structure/ct.e, '-')
ax1.set_xlabel("$k_x$ ($a_B^{-1}$)")
ax1.set_ylabel("Energy (eV)")
ax1.set_ylim([0, 5])
ax1.set_title("Band structure including band mixing")
plt.show()

# Without mixing
# Compute band structure using QTCAD
band_structure = s.band_structure(K, basis_vec, mixing=False)

# Analytic solutions
# Particle in a box energies
def particle_in_box_E(m, L, n):
    return n**2*np.pi**2*ct.hbar**2 / ( 2 * m * L**2 )

g1, g2, g3, Delta = tuple(mat.hole_kp_params)
# Out of plane effective masses
mHH = ct.me / (g1-2*g2) # HH effective mass GaAs
mLH = ct.me / (g1+2*g2) # LH effective mass
mSO = ct.me / g1 # SO effective mass
# In plane effective masses
mHHip = ct.me / (g1 + g2) # HH in-plane effective mass GaAs
mLHip = ct.me / (g1 - g2) # LH in-plane effective mass GaAs
mSOip = ct.me / (g1) # SO in-plane effective mass GaAs

# Analytic energies
PIB_HH = np.array([particle_in_box_E(mHH, 2e-9, n) for n in range(1,5)])
free_HH = ct.hbar**2 / 2 / mHHip * K**2
E_HH = (PIB_HH[np.newaxis, :] + free_HH[:, np.newaxis]) 

PIB_LH = np.array([particle_in_box_E(mLH, 2e-9, n) for n in range(1,5)])
free_LH = ct.hbar**2 / 2 / mLHip * K**2
E_LH = (PIB_LH[np.newaxis, :] + free_LH[:, np.newaxis]) 

PIB_SO = np.array([particle_in_box_E(mSO, 2e-9, n) for n in range(1,5)])
free_SO = ct.hbar**2 / 2 / mSOip * K**2
E_SO = (PIB_SO[np.newaxis, :] + free_SO[:, np.newaxis]) + Delta

# Plot
fig = plt.figure(figsize=(7,4))
ax1 = fig.add_subplot(1,1,1)
ax1.plot(K * ct.aB, band_structure/ct.e, '-')
ax1.plot(K * ct.aB, E_HH/ct.e, '--b') # Plot analytic HH energies in blue
ax1.plot(K * ct.aB, E_LH/ct.e, '--r') # Plot analytic LH energies in blue
ax1.plot(K * ct.aB, E_SO/ct.e, '--g') # Plot analytic SO energies in blue
ax1.set_xlabel("$k_x$ ($a_B^{-1}$)")
ax1.set_ylabel("Energy (eV)")
ax1.set_ylim([0, 5])
ax1.set_title("Band structure in the absence of band mixing")
plt.show()