13. MVEMT applied to a donor in silicon
13.1. Requirements
13.1.1. Software components
13.1.1.1. Basic Python modules
pathlib
numpy
matplotlib
13.1.1.2. Specific modules
13.1.2. Geometry file
qtcad/examples/tutorials/meshes/sphere.geo
13.1.3. Python script
qtcad/examples/tutorials/donor_MVEMT.py
13.1.4. References
13.2. Briefing
In this tutorial we explain how to use the multi-valley effective-mass theory (MVEMT) solver implemented in QTCAD to understand the eigenenergies and eigenstates of a phosphorus donor in silicon. In contrast to the Valley splitting (MVEMT) tutorial which applied the MVEMT solver to a 1D problem, here we will apply it to a 3D problem. To investigate the phosphorus donor, we consider a hydrogenic donor potential with central-cell corrections [PLBL14]
where \(r\) is the distance from the donor, \(\epsilon_{\mathrm {Si}}\) is the dielectric constant of silicon, and \(b\) and \(B\) are central-cell correction parameters. The benefit of using such a potential is that it avoids the singularity at the origin that arises from the Coulomb potential:
It has also been shown, using a variational procedure, that this potential can lead to eigenenergies that are in good agreement with those obtained from experimental measurements [PLBL14]. Qualitatively, the central cell correction potential can be thought of as a potential that reproduces the screened Coulomb interation at long distances, while modifying the short-range behaviour of the potential to account for corrections beyond the point-charge model.
The simulation presented below will consist of solving the MVEMT Schrödinger equation for an electron subject to the donor potential described above in a sphere of silicon. Because the MVEMT solver requires a mesh with atomic resolution [since the valley-orbit coupling potential, \(V^{\mathrm{VO}}(\mathbf r)\) depends on the atomic wavefunctions, see Multivalley effective-mass theory solver], very fine meshes are typically required to converge MVEMT calculations. In order for the script associated with this tutorial to run in a reasonable amount of time, we will use a suboptimal mesh. However, we will also list the mesh parameters that we have used to give reasonable convergence of the eigenenergies and eigenstates.
It is assumed that the user is familiar with basic Python and has some familiarity with the functionality of the QTCAD code.
The full code may be found at the bottom of this page,
or in donor_MVEMT.py.
13.3. Header
In the header, we import relevant modules.
import pathlib
import os
import qtcad
import numpy as np
import matplotlib.pyplot as plt
from scipy.interpolate import LinearNDInterpolator
from qtcad.device.mesh3d import Mesh
from qtcad.device import Device
from qtcad.device import materials as mt
from qtcad.device import constants as ct
from qtcad.device.analysis import linecut
These have been used in previous tutorials.
We also import two additional functions: LinearNDInterpolator from
scipy.interpolate and dump() from joblib.
The first function lets us create a callable function from a linear interpolation of
discrete data, and the second lets us save the generated function to a file
for future use.
from scipy.interpolate import LinearNDInterpolator
from joblib import dump
Just like in the Valley splitting (MVEMT) tutorial, we also import the MVEMT
Solver class which can be used
to solve the Schrödinger equation derived from the multi-valley effective-mass
theory implemented in QTCAD (and described in Multivalley effective-mass theory solver) and the
SolverParams class which
allows us to parametrize the solver.
from qtcad.device.multivalley_EMT import Solver, SolverParams
13.4. Setup
Next we set up the calculation. We start by defining certain constants that will be used throughout the tutorial.
# Constants --------------------------------------------------------------------
h = 0.050 # characteristic length
R = 4 # Radius of sphere [a]
nm = 1 / 10**9 # nanometer
a = 5.431e-10 # lattice constant Si
k0 = 0.84 * 2 * np.pi / a # wavevector
valleys = mt.Si.valleys # Valley wavevectors
me_inv = mt.Si.Me_inv_valley # Inverse effective mass tensors
B = 250 # Potential parameters [units = 1/nm]
b = 30
Here we have defined the nanometer in meters, the lattice constant of silicon, the distance of the conduction-band minima from the \(\Gamma\) point in the Brillouin zone, the position of the valleys in the Brillouin zone, the inverse effective-mass tensors for each valley and the two parameters in units of inverse nanometers that define the central-cell correction potential. The central-cell correction parameters chosen here are those that lead to eigenenergies that we have found to be in good agreement with experimental values. They differ slightly from those used in [PLBL14] (i.e. \(B = 246.1 \mathrm{nm}^{-1}\) and \(b = 19.96 \mathrm{nm}^{-1}\)). This is because, in contrast with this study, we have optimized the parameters of the central-cell correction using the numerical solution of the MVEMT problem rather than a variational procedure using analytic trial wavefunctions.
We also define paths to relevant files and directories that will be used throughout the example.
# Setup -----------------------------------------------------------------------
# File paths
script_dir = pathlib.Path(__file__).parent.resolve()
path_mesh = script_dir / "meshes"
mesh_file = path_mesh / f"sphere_{h:.3f}_{R:.1f}.msh"
path_out = script_dir / "output"
gs_file = path_out / "donor_ground_wavefunction.png"
es_file = path_out / "donor_excited_wavefunction.png"
E_file = path_out / "donor_energies.txt"
dens_file = path_out / f"pc_density_profile_{h:.3f}_{R:.1f}.joblib"
qtcad_dir = qtcad_dir = pathlib.Path(os.path.dirname(qtcad.__file__))
# Bloch amplitudes
vc_path = qtcad_dir / "device" / "valleycoupling"
bloch_path = vc_path / "bloch" / "silicon"
bloch_paths = [
bloch_path/"BA_px_Si.data",
bloch_path/"BA_mx_Si.data",
bloch_path/"BA_py_Si.data",
bloch_path/"BA_my_Si.data",
bloch_path/"BA_pz_Si.data",
bloch_path/"BA_mz_Si.data"
]
13.5. Setting up the device
Next we create a helper function which takes a
Mesh and central-cell correction
parameters, B and b and returns a
Device with a central-cell corrected
impurity potential centered at the origin.
# Create the device ------------------------------------------------------------
def create_device(mesh: Mesh, B: float, b: float) -> Device:
""" Create a device representing a sphere of silicon with a phosphorus donor
at the center. The phosphorus donor is modelled by a hydrogenic potential
with central-cell corrections.
Args:
mesh (Mesh): Mesh object discretizing the silicon sphere
B (float): First central-cell correction potential parameter in units
of inverse nanometers.
b (float): Second central-cell correction potential parameter in units
of inverse nanometers.
Returns:
Device: Device object.
"""
# Create Device
d = Device(mesh, conf_carriers="e")
d.new_region('QD', mt.Si) # Physical volume
B = B / nm # Convert parameter to SI units
b = b / nm # Convert parameter to SI units
# Set impurity potential
def impurity_potential(x, y, z):
r = np.sqrt(x**2+y**2+z**2)
eps_Si = mt.Si.eps
if r > 1e-16:
factor = (1 - np.exp(-b * r) + B * r * np.exp(-b * r)) / r
else:
factor = (B + b) - (b/2 + B) * b * r
U = -ct.e**2/(4 * np.pi * eps_Si) * factor
return U
d.set_V(impurity_potential)
return d
To apply the potential to the device, we define a function of space (x,
y, z) called impurity_potential, and pass this function to the
set_V method of the
Device, d.
The impurity_potential function implements the central-cell corrected
potential described above for all distances r > 1e-16.
For small distances, we use a Taylor expansion of the potential to avoid the
numerical singularity associated with \(0/0\).
13.6. MVEMT calculation
Now that we have a helper function to set up the device, we create an additional function to solve the MVEMT Schrödinger equation (see Multivalley effective-mass theory solver).
# Solving MVEMT Schrödinger equation -------------------------------------------
def solve_MVEMT(d: Device, params: SolverParams=None) -> None:
""" Solve the multivalley effective mass theory for the device.
Args:
d (Device): Device for which we wish to solve the MVEMT.
params (SolverParams): Parameters for the solver.
"""
if params is None:
params = SolverParams()
params.num_states = 10
params.approx = None
params.val = valleys
params.bloch_files = bloch_paths
params.m_inv_tensors = me_inv
params.lattice_const = a
params.method = "fast"
params.tol = 1e-8
params.maxiter = 5000
mvemt = Solver(d, solver_params=params)
mvemt.solve()
# Energy output
print("Valley-coupled eigenenergies [eV]:")
energies_mvemt = d.energies/ct.e
print(energies_mvemt)
E_out = np.append(np.array([h, R, B, b]), energies_mvemt)
E_header = '' # set empty header
if not os.path.isfile(E_file): # checks if the file exists
# if it doesn't then add the header
E_header = "h, R [a], B [1/nm], b [1/nm], energies [eV]"
with open(E_file, 'a') as file:
np.savetxt(file, np.array(E_out)[np.newaxis, :], header=E_header)
The solve_MVEMT function takes in a device d and, optionally, a
SolverParams object and
applies the MVEMT Solver to the
device.
By default (if no
SolverParams object is
passed to solve_MVEMT), the function sets up the solver parameters such
that the solver will solve for 10 states, using the Bloch amplitudes from the
files in bloch_paths (without any approximations, see e.g.
Valley splitting (MVEMT)), with the material properties (lattice constant,
valley positions, inverse effective mass tensors for the different valleys)
defined earlier.
The default SolverParams
also sets the solving method to "fast", with a tolerance of 1e-8, and
sets the maximum number of iterations of the "fast" iterative eigensolver
to 5000.
After solving the MVEMT equations, the function prints the eigenenergies and saves them to a file.
13.7. Visualizing the results
Next we create a function, plot_psi, to plot the contributions to the
wavefunctions from each valley.
# Visualizing the results ------------------------------------------------------
def plot_psi(wf: np.ndarray, mesh: Mesh, path_fig: pathlib.Path, id: str) -> None:
"""Plot the contributions to the wavefunctions from each valley.
Args:
wf (np.ndarray): Wavefunction computed from a MVEMT solver.
mesh (Mesh): Mesh object over which the wavefunctions are defined.
path_fig (pathlib.Path): Path to save the plot.
id (str): Identifier for the state.
"""
# Linecut through the center of the sphere
R0 = np.max(mesh.glob_nodes) # Radius of the domain
begin = np.array([-1, -1, -1]) * R0 / np.sqrt(3)
end = np.array([1, 1, 1]) * R0 / np.sqrt(3)
# Linecut through the center of the sphere
for i in range(6): # loop over all valley contributions
x, psi_i = linecut(mesh, wf[:, i], begin, end)
try:
psi_lc = np.vstack((psi_lc, psi_i))
except NameError:
psi_lc = psi_i
# Labels for the valley states
labels = ['+x', '-x', '+y', '-y', '+z', '-z']
# Create the figure
fig, axs = plt.subplots(3, 2, figsize=(10, 10))
# Plotting the valley contributions
for i in range(6):
ax = axs[i//2, i%2]
label = labels[i]
ax.plot(x / nm, np.real(psi_lc[i]), label=f'Re[$F^{id}_{{{label}}}$]')
ax.plot(x / nm, np.imag(psi_lc[i]), label=f'Im[$F^{id}_{{{label}}}$]')
ax.set_title(f'{label} valley')
ax.set_xlabel('$x$ [nm]')
ax.set_ylabel(f'$F_{{{label}}}$' + ' [1/m$^{3/2}$]')
ax.legend()
ax.grid()
fig.suptitle(f'Wavefunction contributions for state {id}')
fig.tight_layout()
plt.savefig(path_fig)
This function takes the wavefunctions, wf, and mesh over which the
wavefunction is defined, mesh, and saves a plot of a linecut of the
projection of the wavefunction onto each valley.
13.8. Calculation
We now use the functions defined above to perform a donor calculation.
# Calculation ------------------------------------------------------------------
mesh = Mesh(nm, mesh_file)
d = create_device(mesh, B, b)
solve_MVEMT(d)
We start by creating the Mesh and the
Device (using the create_device
function).
Then we use the solve_MVEMT function to solve the MVEMT Schrödinger
equation.
This function prints the following eigenenergies to the console:
Valley-coupled eigenenergies [eV]:
[-0.10373137 -0.09200582 -0.091813 -0.09160682 -0.08987396 -0.08959993
-0.03366751 -0.03351823 -0.03343959 -0.03322084]
We also plot linecuts of the projection of the ground and first excited states onto each valley component and save them to files using
plot_psi(d.eigenfunctions[:, 0, :], mesh, gs_file, '0')
plot_psi(d.eigenfunctions[:, 1, :], mesh, es_file, '1')
The ground state envelope functions are shown in the figure below.
In addition, the first excited state envelope functions are shown in the figure below.
To get an idea of the qualitative validity of the calculation, we consider the Hamiltonian for a donor under the envelope function approximation [see Eq. (5.2) with \(V_{\mathrm{conf}}(\mathbf r)\) replaced by the central-cell corrected donor potential \(V(\mathbf r)\) defined above [PLBL14]]. Because the potential is spherically symmetric and the six valleys in silicon are equivalent, we expect the ground subspace to be sixfold degenerate. We write these six degenerate states in vector notation [see (5.1)]
where the \(F_{\nu}(\mathbf r)\) are the envelope functions for the valleys and are related to \(F_{\nu'}(\mathbf r)\) (for \(\nu \ne \nu'\)) by inversion, if \(\nu = -\nu'\), or by rotations of \(\pi/2\) about the \(x\), \(y\), and \(z\) axes, if \(\nu \ne -\nu'\). Once we turn on valley coupling, the spherical symmetry gets reduced to tetrahedral as the information about the crystal is included in the Hamiltonian [via the Bloch functions, see Eq. (5.4)]. The sixfold degeneracy gets lifted and is replaced by a singly-degenerate ground state labelled \(A_1\), triply-degenerate excited states labelled \(T_x\), \(T_y\), and \(T_z\), and doubly-degenerate second excited states labelled \(E_z\) and \(E_{xy}\) [SBCalderonK15]. The energies of these states have been measured experimentally [JGR81]:
The computed energy levels are also grouped as a singly-degenerate ground state, (nearly) triply-degenerate excited states, and (nearly) doubly-degenerate second excited states, as expected from group theory. Moreover, while the computed energies are quite far from the experimental values, the energy differences between the states and the degeneracies are in good agreement.
Note
While the calculation presented here runs in a reasonable amount of time
(less than 10 minutes) using limited computational resources (a laptop), it
is important to note that to achieve this performance, the calculation was
run over an unconverged mesh.
Ideally, the outputs of the calculation should be converged with respect to
the characteristic length h and the size of the geometry R.
Looking at the figures above, we see that the wavefunctions do not vanish
at the edges of the domain indicating that we may be far from convergence.
We have observed rough convergence (using a workstation over hours) with
h = 0.075 and R = 12 where we find
\(E_{A_1} = -0.0446 \, \mathrm{eV}\),
\(E_{T} = -0.0365 \, \mathrm{eV}\), and
\(E_{E} = -0.0356 \, \mathrm{eV}\) with
\(B = 250 \, \mathrm{nm}^{-1}\) and \(b = 30 \, \mathrm{nm}^{-1}\).
Group theory also tells us what the eigenstates in this subspace have the form:
With these states in mind, we can get some qualitative understanding of the plots above. We start with the ground state which we identify with \(A_1\) since it is singly-degenerate. We expect this state to be a symmetric linear combination of all six valley states, \(\vec{F}^{\nu}(\mathbf{r})\), \(\nu \in \{\pm x, \pm y, \pm z\}\). Because of the properties described above (inversion and rotation symmetries), we expect the contributions from the \(+x\), \(+y\), and \(+z\) i.e. \(F_{+x}(\mathbf r)\), \(F_{+y}(\mathbf r)\), and \(F_{+z}(\mathbf r)\) to be identical along a linecut from \((-1, -1, -1)\) to \((1, 1, 1)\) (which treats the \(x\), \(y\), and \(z\) directions on equal footing) and to be related to the contributions from the \(-x\), \(-y\), and \(-z\), i.e \(F_{-x}(\mathbf r)\), \(F_{-y}(\mathbf r)\), and \(F_{-z}(\mathbf r)\) valleys by inversion symmetry along the linecut. This behavior is seen in the plot of the ground state wavefunction above.
We next consider the excited state. Because the excited states are triply-degenerate, the result of the MVEMT calculation is a linear combination of the \(T_x\), \(T_y\), and \(T_z\) states. Any linear combination of these states is expected to have a contribution from a valley \(\nu\) that is related to the contribution from the valley \(-\nu\) by a minus sign. This behaviour is also observed in the plots above.
13.9. Saving the result
The calculation described above can also be used to understand how point-charge impurities influence a system. If the impurities are phosphorus dopants, we expect an electron bound (or not) to the dopant to add (or remove) the charge density associated with ground state computed above. In this section, we create an interpolation function from the computed ground state density and save it to a file. This density can then be loaded as a point-charge density profile to model phosphorus-donor charge impurities in a device.
# Save -------------------------------------------------------------------------
# A1 state
psi_A1 = d.eigenfunctions[:, 0, :]
dens_A1 = np.sum(np.abs(psi_A1)**2, axis=1)
coords = mesh.glob_nodes
# Create interpolation function
interpolator = LinearNDInterpolator(coords, dens_A1, fill_value=0)
# Save interpolation
dump(interpolator, dens_file)
In the above code, we have used the LinearNDInterpolator from
scipy.interpolate to create a callable object that can be used to
interpolate the ground-state density over all points in space (linear
interpolation over the region of space where the mesh is defined and zero
elsewhere).
We also use the dump() function from joblib to save the
interpolation function to a file for later use
(see Point charges in a double quantum dot in FD-SOI).
13.10. Full code
__copyright__ = "Copyright 2022-2024, Nanoacademic Technologies Inc."
import pathlib
import os
import qtcad
import numpy as np
import matplotlib.pyplot as plt
from scipy.interpolate import LinearNDInterpolator
from qtcad.device.mesh3d import Mesh
from qtcad.device import Device
from qtcad.device import materials as mt
from qtcad.device import constants as ct
from qtcad.device.analysis import linecut
from scipy.interpolate import LinearNDInterpolator
from joblib import dump
from qtcad.device.multivalley_EMT import Solver, SolverParams
# Constants --------------------------------------------------------------------
h = 0.050 # characteristic length
R = 4 # Radius of sphere [a]
nm = 1 / 10**9 # nanometer
a = 5.431e-10 # lattice constant Si
k0 = 0.84 * 2 * np.pi / a # wavevector
valleys = mt.Si.valleys # Valley wavevectors
me_inv = mt.Si.Me_inv_valley # Inverse effective mass tensors
B = 250 # Potential parameters [units = 1/nm]
b = 30
# Setup -----------------------------------------------------------------------
# File paths
script_dir = pathlib.Path(__file__).parent.resolve()
path_mesh = script_dir / "meshes"
mesh_file = path_mesh / f"sphere_{h:.3f}_{R:.1f}.msh"
path_out = script_dir / "output"
gs_file = path_out / "donor_ground_wavefunction.png"
es_file = path_out / "donor_excited_wavefunction.png"
E_file = path_out / "donor_energies.txt"
dens_file = path_out / f"pc_density_profile_{h:.3f}_{R:.1f}.joblib"
qtcad_dir = qtcad_dir = pathlib.Path(os.path.dirname(qtcad.__file__))
# Bloch amplitudes
vc_path = qtcad_dir / "device" / "valleycoupling"
bloch_path = vc_path / "bloch" / "silicon"
bloch_paths = [
bloch_path/"BA_px_Si.data",
bloch_path/"BA_mx_Si.data",
bloch_path/"BA_py_Si.data",
bloch_path/"BA_my_Si.data",
bloch_path/"BA_pz_Si.data",
bloch_path/"BA_mz_Si.data"
]
# Create the device ------------------------------------------------------------
def create_device(mesh: Mesh, B: float, b: float) -> Device:
""" Create a device representing a sphere of silicon with a phosphorus donor
at the center. The phosphorus donor is modelled by a hydrogenic potential
with central-cell corrections.
Args:
mesh (Mesh): Mesh object discretizing the silicon sphere
B (float): First central-cell correction potential parameter in units
of inverse nanometers.
b (float): Second central-cell correction potential parameter in units
of inverse nanometers.
Returns:
Device: Device object.
"""
# Create Device
d = Device(mesh, conf_carriers="e")
d.new_region('QD', mt.Si) # Physical volume
B = B / nm # Convert parameter to SI units
b = b / nm # Convert parameter to SI units
# Set impurity potential
def impurity_potential(x, y, z):
r = np.sqrt(x**2+y**2+z**2)
eps_Si = mt.Si.eps
if r > 1e-16:
factor = (1 - np.exp(-b * r) + B * r * np.exp(-b * r)) / r
else:
factor = (B + b) - (b/2 + B) * b * r
U = -ct.e**2/(4 * np.pi * eps_Si) * factor
return U
d.set_V(impurity_potential)
return d
# Solving MVEMT Schrödinger equation -------------------------------------------
def solve_MVEMT(d: Device, params: SolverParams=None) -> None:
""" Solve the multivalley effective mass theory for the device.
Args:
d (Device): Device for which we wish to solve the MVEMT.
params (SolverParams): Parameters for the solver.
"""
if params is None:
params = SolverParams()
params.num_states = 10
params.approx = None
params.val = valleys
params.bloch_files = bloch_paths
params.m_inv_tensors = me_inv
params.lattice_const = a
params.method = "fast"
params.tol = 1e-8
params.maxiter = 5000
mvemt = Solver(d, solver_params=params)
mvemt.solve()
# Energy output
print("Valley-coupled eigenenergies [eV]:")
energies_mvemt = d.energies/ct.e
print(energies_mvemt)
E_out = np.append(np.array([h, R, B, b]), energies_mvemt)
E_header = '' # set empty header
if not os.path.isfile(E_file): # checks if the file exists
# if it doesn't then add the header
E_header = "h, R [a], B [1/nm], b [1/nm], energies [eV]"
with open(E_file, 'a') as file:
np.savetxt(file, np.array(E_out)[np.newaxis, :], header=E_header)
# Visualizing the results ------------------------------------------------------
def plot_psi(wf: np.ndarray, mesh: Mesh, path_fig: pathlib.Path, id: str) -> None:
"""Plot the contributions to the wavefunctions from each valley.
Args:
wf (np.ndarray): Wavefunction computed from a MVEMT solver.
mesh (Mesh): Mesh object over which the wavefunctions are defined.
path_fig (pathlib.Path): Path to save the plot.
id (str): Identifier for the state.
"""
# Linecut through the center of the sphere
R0 = np.max(mesh.glob_nodes) # Radius of the domain
begin = np.array([-1, -1, -1]) * R0 / np.sqrt(3)
end = np.array([1, 1, 1]) * R0 / np.sqrt(3)
# Linecut through the center of the sphere
for i in range(6): # loop over all valley contributions
x, psi_i = linecut(mesh, wf[:, i], begin, end)
try:
psi_lc = np.vstack((psi_lc, psi_i))
except NameError:
psi_lc = psi_i
# Labels for the valley states
labels = ['+x', '-x', '+y', '-y', '+z', '-z']
# Create the figure
fig, axs = plt.subplots(3, 2, figsize=(10, 10))
# Plotting the valley contributions
for i in range(6):
ax = axs[i//2, i%2]
label = labels[i]
ax.plot(x / nm, np.real(psi_lc[i]), label=f'Re[$F^{id}_{{{label}}}$]')
ax.plot(x / nm, np.imag(psi_lc[i]), label=f'Im[$F^{id}_{{{label}}}$]')
ax.set_title(f'{label} valley')
ax.set_xlabel('$x$ [nm]')
ax.set_ylabel(f'$F_{{{label}}}$' + ' [1/m$^{3/2}$]')
ax.legend()
ax.grid()
fig.suptitle(f'Wavefunction contributions for state {id}')
fig.tight_layout()
plt.savefig(path_fig)
# Calculation ------------------------------------------------------------------
mesh = Mesh(nm, mesh_file)
d = create_device(mesh, B, b)
solve_MVEMT(d)
plot_psi(d.eigenfunctions[:, 0, :], mesh, gs_file, '0')
plot_psi(d.eigenfunctions[:, 1, :], mesh, es_file, '1')
# Save -------------------------------------------------------------------------
# A1 state
psi_A1 = d.eigenfunctions[:, 0, :]
dens_A1 = np.sum(np.abs(psi_A1)**2, axis=1)
coords = mesh.glob_nodes
# Create interpolation function
interpolator = LinearNDInterpolator(coords, dens_A1, fill_value=0)
# Save interpolation
dump(interpolator, dens_file)