12. Valley splitting (MVEMT)

12.1. Requirements

12.1.1. Software components

12.1.1.1. Basic Python modules

  • pathlib

  • numpy

  • matplotlib

12.1.1.2. Specific modules

12.1.2. Geometry file

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

12.1.3. Python script

  • qtcad/examples/tutorials/valleysplitting_1D_MVEMT.py

12.1.4. References

12.2. Briefing

In this tutorial we explain how to calculate the valley splitting in a silicon quantum well using the multi-valley effective-mass theory (MVEMT) solver implemented in QTCAD. In this example, we consider an electron confined to a \(\mathrm{Si/SiO}_2\) interface with an electric field. We take this confinement to be along the \(x\) direction and perform a one-dimensional calculation. This example aims to use QTCAD’s finite-element MVEMT solver to reproduce the results of [GHCJ+16] (where MVEMT is discretized over a Gaussian basis), to understand the valley splitting in a silicon quantum well. Specifically, we will use QTCAD to reproduce Fig. S1 (a) of the supplemental material of [GHCJ+16] (reproduced as Fig. 12.2.1 of this tutorial).

Reproduction of Fig. S1 (a) from the supplemental material of :cite:`Gamble2016valley`

Fig. 12.2.1 Valley splitting in a silicon quantum well as a function of electric field. This figure is a reproduction of Fig. S1 (a) from the supplemental material of [GHCJ+16].

It is assumed that the user is familiar with basic Python and has some familiarity with the functionality of the QTCAD software.

The full code may be found at the bottom of this page, or in valleysplitting_1D_MVEMT.py.

12.4. Setup

Now that the necessary modules have been imported, we define paths to relevant files and/or directories and the constants that will be used throughout the example. We start with paths to the mesh file and the output directory, to two files which will contain the plots this tutorial will generate, and to three files which will contain the data generated in this tutorial.

# Setup -----------------------------------------------------------------------
# File paths
script_dir  = pathlib.Path(__file__).parent.resolve()
path_mesh   = script_dir / "meshes" / "Si_well.msh"
path_out    = script_dir / "output"

# Outputs
## Plots
path_valley = path_out / f"valley_wavefunctions.png"
path_vs_E   = path_out / f"valley_splitting_vs_E.png"
## Data
path_triv   = path_out / f"valley_splitting_qw_trivial.txt"
path_ff     = path_out / f"valley_splitting_qw_ff.txt"
path_None   = path_out / f"valley_splitting_qw_None.txt"

We then specify the path to files containing information about the Bloch amplitudes of silicon. Specifically, these files contain a description of the Bloch amplitudes in the plane-wave basis. The Bloch amplitudes contain a description of the periodic part of the wavefunction of the crystal Hamiltonian [see Multivalley effective-mass theory solver for more details].

# Inputs
## Bloch amplitudes
qtcad_dir   = pathlib.Path(os.path.dirname(qtcad.__file__))
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"
]

Finally we specify constants that will be used in the calculation

## Constants
scale       = 1e-9                      # nm
a           = 0.543 * scale             # Lattice constant of Si
k0          = 0.84 * 2 * np.pi / a      # Wavevector
valleys     = mt.Si.valleys[0:6]        # Valley wavevectors
me_inv      = mt.Si.Me_inv_valley[0:6]  # Inverse effective mass tensors

These are mainly material parameters for silicon. They include the lattice constant, the magnitude of the wavevectors where the conduction-band minima are located in the first Brillouin zone, the location of each individual valley in the first Brillouin zone, and the inverse effective mass tensors for each valley.

12.5. Setting up the device

Next we create a helper function which takes in a value of the electic field and outputs a device defined over the mesh located at path_mesh.

# Create the device ------------------------------------------------------------
def create_well(E) -> Device:
    """ Creates a quantum-well device with a given electric field.

    Args:
        E (float): Electric field in V/m.

    Returns:
        Device: Device object.
    """

    # Mesh
    mesh = Mesh(scale, path_mesh)

    # Device
    d = Device(mesh, conf_carriers="e")
    d.new_region('Domain', mt.Si)
    d.new_region('Barrier', mt.SiO2)

    # Set potential
    def V(z):
        v = -ct.e * E * z

        return v

    d.set_V(V)
    d.add_to_V(3 * ct.e, "Barrier") # Add barrier offset

    return d

The create_well function creates a Device composed of silicon from \(x=-13 \, \mathrm{nm}\) to \(x=0\) and silicon dioxide from \(x=0\) to \(x=3 \, \mathrm{nm}\). These dimensions are specified in the geometry file qtcad/examples/tutorials/meshes/qw_1d_Si.geo. In addition, there is an electric field E applied in the \(x\) direction. The potential is set to be linear in the \(x\) direction with a slope proportional to the electric field E. In addition, the potential is raised by \(3 \, \mathrm{eV}\) in the barrier region which represents the band offset between silicon and silicon dioxide. The device is chosen to have the above structure to mimic the quantum well studied in [GHCJ+16].

12.6. MVEMT calculation

Now that we have a helper function to set up the device for a given electric field, E, we create an additional helper function to solve the MVEMT Schrödinger equation (see Multivalley effective-mass theory solver).

# Solving MVEMT Schrödinger equation -------------------------------------------
def solve_MVEMT(d: Device, approx: str=None) -> None:
    """ Solve the multivalley effective mass theory for the device.

    Args:
        d (Device): Device for which we wish to solve the MVEMT.
        approx (None or str): Level of approximation to use in the MVEMT.
            - "trivial" for trivial approximation.
            - "ff" for form-factor approximation.
            - None for full calculation.
    """

    # Compute valley coupling (using form-factor approximation)
    params               = SolverParams()
    params.num_states    = 20
    params.approx        = approx
    params.val           = valleys
    params.bloch_files   = bloch_paths
    params.m_inv_tensors = me_inv
    params.lattice_const = a
    params.method        = "fast"
    params.tol           = 1e-7
    params.maxiter       = 5000

    mvemt = Solver(d, solver_params=params)
    mvemt.solve()
    d.print_energies()

The solve_MVEMT function takes in a device d and an approximation approx and solves the MVEMT equations over the device under the approximation approx. The possible values for approx that QTCAD can currently handle are

  • "trivial" for the trivial approximation where the Bloch amplitudes are set to unity, \(u_{\nu}(\mathbf{r}) = 1\),

  • "ff" for the form-factor approximation where products of Bloch amplitudes are approximated by constants, \(u_{\nu}^{\ast}(\mathbf{r})u_{\nu'}(\mathbf{r}) = a_0^{\nu\nu'}\) where \(a_0^{\nu\nu'}\) are constants that depend only on the valleys \(\nu\) and \(\nu'\),

  • None where the Bloch amplitudes are treated exactly.

The function first sets up the solver parameters using a SolverParams object. The number of states to solve for is set to 20, and the approximation to use is set to approx. We also pass the material properties defined earlier, set the solving method to "fast", set the tolerance to 1e-7, and the maximum number of iterations of the "fast" iterative eigensolver to 5000. The solver parameters are then passed to the Solver object which is used to solve the MVEMT equations.

12.7. Computing the valley splitting

As mentioned above, we aim to reproduce the results of [GHCJ+16], where the authors compute the valley splitting of an electron bound to a silicon quantum well formed at the Si/SiO2 interface of a heterostructure using the MVEMT. This calculation is performed as a function an applied electric field and under the various approximation schemes listed above ("trivial", "ff", and None).

# Calculation -----------------------------------------------------------------

# Levels of approximation
approximations = ['trivial', 'ff', None]
for approx in approximations:

    # Electric field values
    E_fields = np.linspace(5, 50, 10)
    for E in E_fields:

        d = create_well(E * 1e6)
        solve_MVEMT(d, approx)

        save_file = path_out / f"valley_splitting_qw_{approx}.txt"
        with open(save_file, 'a') as file:
            out = np.zeros((1, 2))
            out[0, 0] = E
            out[0, 1] = (d.energies[1] - d.energies[0]) / ct.e
            np.savetxt(file, out)

Here we loop over the various approximation schemes and electric field values and for each we perform a MVEMT calculation. We expect the valley splitting to be much smaller than the orbital level spacing generated by the quantum-well confinement. Accordingly, we can identify the valley splitting as the energy difference between the ground and first excited eigenenergies. We save this splitting, along with the electric field value to a file for each approximation scheme.

12.8. Visualizing the results

As discussed in Multivalley effective-mass theory solver, the eigenstates of the MVEMT Hamiltonian can be thought of as vectors whose components span the different valleys [see in particular Eq. (5.6)].

Note

The MVEMT Solver in QTCAD computes eigenstates as 3 dimensional array-like objects stored in the eigenfunctions attribute of a Device. The indices of the array correspond to the global nodes of the mesh, the different eigenstates, and the valley contributions, respectively. This ordering is similar to that used for holes, where we have replaced the band index by a valley index.

Due to the confinement along the \(x\) direction and the larger effective mass of the \(\pm x\) valleys along the \(x\) direction relative to the \(\pm y\) and \(\pm z\) valleys, we expect the lowest energy states to possess predominantly \(\pm x\) valley character. Here we plot the valley contributions to the ground state and first-excited state wavefunctions.

# Visualizing the results ------------------------------------------------------

def plot_env(wf: np.ndarray, mesh: Mesh, path_fig: pathlib.Path) -> None:
    """Plot the envelope contributions to the wavefunctions from each valley.

    Args:
        wf (np.ndarray): Wavefunctions computed from a MVEMT solver.
        mesh (Mesh): Mesh object over which the wavefunctions are defined.
        path_fig (pathlib.Path): Path to save the plot.
    """
    # Sort the wavefunctions over the global nodes of the mesh
    ix  = np.argsort(mesh.glob_nodes, axis=0)[:, 0]       # index
    x   = mesh.glob_nodes[ix, 0] / scale                  # sorted nodes
    env = wf[ix, :, :]                                    # sorted wf

    # 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]
        ax2 = ax.twinx()
        for j in range(env.shape[1]):
            ax.plot(x, np.real(env[:, j, i]), label=f'Re[$F_{"{"+labels[i]+"}"}^{j}$]')
            ax.plot(x, np.imag(env[:, j, i]), label=f'Im[$F_{"{"+labels[i]+"}"}^{j}$]')
            ax2.plot(x, np.abs(env[:, j, i])**2, '--', label=f'$|F_{"{"+labels[i]+"}"}^{j}|^2$')
        ax.set_title(f'{labels[i]} valley')
        ax.set_xlabel('$x$ [nm]')
        ax.set_ylabel(f'$F_{"{"+labels[i]+"}"}$ [1/m$^{{3/2}}$]')
        ax2.set_ylabel(f'$|F_{"{"+labels[i]+"}"}|^2$ [1/m$^{{3}}$]')
        ax.legend(loc='upper left')
        ax2.legend(loc='lower left')
        ax.grid()

    fig.tight_layout()
    plt.savefig(path_fig)

plot_env(d.eigenfunctions[:, 0:2, :], d.mesh, path_valley)
Ground and first-excited state.

Fig. 12.8.1 Contributions to the ground and first-excited state wavefunctions from each valley in a silicon quantum well formed at a Si/SiO2 interface. The wavefunctions are displayed using solid lines while the densities are displayed using dashed lines. Due to the difference in effective mass along the x direction, the ground and first-excited states are predominantly composed of the \(\pm x\) valleys. The other valleys have contributions that are negligible compared to the \(\pm x\) valleys.

In Fig. 12.8.1 we plot the projection of the ground and first-excited states (\(\vec F^0(\mathbf r)\) and \(\vec F^1(\mathbf r)\), respectively) onto each valley states using solid lines. We also plot the density associated with these projections with dashed lines. As expected, the ground and first-excited states are predominantly composed of the \(\pm x\) valleys due to the larger effective mass of the \(\pm x\) valleys along the \(x\) direction relative to the \(\pm y\) and \(\pm z\) valleys. We also see that the ground and first-excited states have identical densities associated with the \(+ x\) valley and \(-x\) valleys (the dashed lines are on top of each other), as expected since the sharp interface breaks the envelope function approximation and provides a coupling between the valley states.

Finally, we also plot the valley splitting as a function of the electric field for the various approximation schemes.

def plot_vs(data_None: np.ndarray, data_ff: np.ndarray, data_triv: np.ndarray,
        save_file: pathlib.Path) -> None:
    """ Plot the valley splitting for the different approximations.
    Args:
        data_None (np.ndarray): Valley splitting data for full calculation.
        data_ff (np.ndarray): Valley splitting data for form-factor
            approximation.
        data_triv (np.ndarray): Valley splitting data for trivial
            approximation.
        save_file (pathlib.Path): Path to save the plot.
    """

    # Create Figure
    fig, ax = plt.subplots(nrows=1, ncols=1, layout="tight")
    fig.set_size_inches(8, 4)

    # Plot data in [meV]
    ax.plot(data_ff[:, 0], data_ff[:, 1] * 1000, '-r', linewidth=2, label="form factors")
    ax.plot(data_triv[:, 0], data_triv[:, 1] * 1000, '-k', linewidth=2, label="trivial")
    ax.plot(data_None[:, 0], data_None[:, 1] * 1000, '-g', linewidth=2, label="full")

    # Settings
    ax.set_title(f"Fig. (S1) [Gamble et al. APL 109 (2016)] using QTCAD", fontsize=18)
    ax.set_xlabel("Electric field [MV/m]", fontsize=16)
    ax.set_ylabel("valley splitting [meV]", fontsize=16)

    ax.tick_params(axis='x', labelsize=14)
    ax.tick_params(axis='y', labelsize=14)

    ax.legend(fontsize=16)
    ax.grid()

    # Save
    plt.savefig(save_file, dpi=300)

# Reload Data
data_None = np.loadtxt(path_None)
data_ff = np.loadtxt(path_ff)
data_triv = np.loadtxt(path_triv)

# Plot valley splitting vs. E
plot_vs(data_None, data_ff, data_triv, path_vs_E)

This part of the script outputs the equivalent of Fig. S1 (a) from [GHCJ+16].

Valley splitting as a function of electric field.

Fig. 12.8.2 Valley splitting as a function of electric field for an electron confined to a quantum well formed at a Si/SiO2 interface. The calculation is performed using the MVEMT under various approximation schemes – 'trivial', 'ff', and None. The corresponding results are in good agreement with those of [GHCJ+16]. Note that the green line overlaps with the red line.

Comparing Fig. 12.8.2 with Fig. S1 (a) of [GHCJ+16], we see that the results obtained with QTCAD are in good agreement with those of [GHCJ+16].

12.9. 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 qtcad.device.mesh1d import Mesh
from qtcad.device import Device
from qtcad.device import materials as mt
from qtcad.device import constants as ct

from qtcad.device.multivalley_EMT import Solver, SolverParams

# Setup -----------------------------------------------------------------------
# File paths
script_dir  = pathlib.Path(__file__).parent.resolve()
path_mesh   = script_dir / "meshes" / "Si_well.msh"
path_out    = script_dir / "output"

# Outputs
## Plots
path_valley = path_out / f"valley_wavefunctions.png"
path_vs_E   = path_out / f"valley_splitting_vs_E.png"
## Data
path_triv   = path_out / f"valley_splitting_qw_trivial.txt"
path_ff     = path_out / f"valley_splitting_qw_ff.txt"
path_None   = path_out / f"valley_splitting_qw_None.txt"

# Inputs
## Bloch amplitudes
qtcad_dir   = pathlib.Path(os.path.dirname(qtcad.__file__))
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"
]

## Constants
scale       = 1e-9                      # nm
a           = 0.543 * scale             # Lattice constant of Si
k0          = 0.84 * 2 * np.pi / a      # Wavevector
valleys     = mt.Si.valleys[0:6]        # Valley wavevectors
me_inv      = mt.Si.Me_inv_valley[0:6]  # Inverse effective mass tensors 

# Create the device ------------------------------------------------------------
def create_well(E) -> Device:
    """ Creates a quantum-well device with a given electric field.

    Args:
        E (float): Electric field in V/m.

    Returns:
        Device: Device object.
    """

    # Mesh
    mesh = Mesh(scale, path_mesh)

    # Device
    d = Device(mesh, conf_carriers="e")
    d.new_region('Domain', mt.Si)
    d.new_region('Barrier', mt.SiO2)

    # Set potential
    def V(z):
        v = -ct.e * E * z

        return v

    d.set_V(V)
    d.add_to_V(3 * ct.e, "Barrier") # Add barrier offset

    return d

# Solving MVEMT Schrödinger equation -------------------------------------------
def solve_MVEMT(d: Device, approx: str=None) -> None:
    """ Solve the multivalley effective mass theory for the device.

    Args:
        d (Device): Device for which we wish to solve the MVEMT.
        approx (None or str): Level of approximation to use in the MVEMT.
            - "trivial" for trivial approximation.
            - "ff" for form-factor approximation.
            - None for full calculation.
    """

    # Compute valley coupling (using form-factor approximation)
    params               = SolverParams()
    params.num_states    = 20
    params.approx        = approx
    params.val           = valleys
    params.bloch_files   = bloch_paths
    params.m_inv_tensors = me_inv
    params.lattice_const = a
    params.method        = "fast"
    params.tol           = 1e-7
    params.maxiter       = 5000

    mvemt = Solver(d, solver_params=params)
    mvemt.solve()
    d.print_energies()

# Calculation -----------------------------------------------------------------

# Levels of approximation
approximations = ['trivial', 'ff', None]
for approx in approximations:

    # Electric field values
    E_fields = np.linspace(5, 50, 10)
    for E in E_fields:

        d = create_well(E * 1e6)
        solve_MVEMT(d, approx)

        save_file = path_out / f"valley_splitting_qw_{approx}.txt"
        with open(save_file, 'a') as file:
            out = np.zeros((1, 2))
            out[0, 0] = E 
            out[0, 1] = (d.energies[1] - d.energies[0]) / ct.e
            np.savetxt(file, out)

# Visualizing the results ------------------------------------------------------

def plot_env(wf: np.ndarray, mesh: Mesh, path_fig: pathlib.Path) -> None:
    """Plot the envelope contributions to the wavefunctions from each valley.

    Args:
        wf (np.ndarray): Wavefunctions computed from a MVEMT solver.
        mesh (Mesh): Mesh object over which the wavefunctions are defined.
        path_fig (pathlib.Path): Path to save the plot.
    """
    # Sort the wavefunctions over the global nodes of the mesh
    ix  = np.argsort(mesh.glob_nodes, axis=0)[:, 0]       # index
    x   = mesh.glob_nodes[ix, 0] / scale                  # sorted nodes
    env = wf[ix, :, :]                                    # sorted wf

    # 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]
        ax2 = ax.twinx()
        for j in range(env.shape[1]):
            ax.plot(x, np.real(env[:, j, i]), label=f'Re[$F_{"{"+labels[i]+"}"}^{j}$]')
            ax.plot(x, np.imag(env[:, j, i]), label=f'Im[$F_{"{"+labels[i]+"}"}^{j}$]')
            ax2.plot(x, np.abs(env[:, j, i])**2, '--', label=f'$|F_{"{"+labels[i]+"}"}^{j}|^2$')
        ax.set_title(f'{labels[i]} valley')
        ax.set_xlabel('$x$ [nm]')
        ax.set_ylabel(f'$F_{"{"+labels[i]+"}"}$ [1/m$^{{3/2}}$]')
        ax2.set_ylabel(f'$|F_{"{"+labels[i]+"}"}|^2$ [1/m$^{{3}}$]')
        ax.legend(loc='upper left')
        ax2.legend(loc='lower left')
        ax.grid()

    fig.tight_layout()  
    plt.savefig(path_fig)

plot_env(d.eigenfunctions[:, 0:2, :], d.mesh, path_valley)

def plot_vs(data_None: np.ndarray, data_ff: np.ndarray, data_triv: np.ndarray,
        save_file: pathlib.Path) -> None:
    """ Plot the valley splitting for the different approximations.
    Args:
        data_None (np.ndarray): Valley splitting data for full calculation.
        data_ff (np.ndarray): Valley splitting data for form-factor 
            approximation.
        data_triv (np.ndarray): Valley splitting data for trivial 
            approximation.
        save_file (pathlib.Path): Path to save the plot.
    """ 

    # Create Figure
    fig, ax = plt.subplots(nrows=1, ncols=1, layout="tight")
    fig.set_size_inches(8, 4)

    # Plot data in [meV] 
    ax.plot(data_ff[:, 0], data_ff[:, 1] * 1000, '-r', linewidth=2, label="form factors")
    ax.plot(data_triv[:, 0], data_triv[:, 1] * 1000, '-k', linewidth=2, label="trivial")
    ax.plot(data_None[:, 0], data_None[:, 1] * 1000, '-g', linewidth=2, label="full")

    # Settings
    ax.set_title(f"Fig. (S1) [Gamble et al. APL 109 (2016)] using QTCAD", fontsize=18)
    ax.set_xlabel("Electric field [MV/m]", fontsize=16)
    ax.set_ylabel("valley splitting [meV]", fontsize=16)

    ax.tick_params(axis='x', labelsize=14)
    ax.tick_params(axis='y', labelsize=14)

    ax.legend(fontsize=16)
    ax.grid()

    # Save
    plt.savefig(save_file, dpi=300)

# Reload Data
data_None = np.loadtxt(path_None)
data_ff = np.loadtxt(path_ff)
data_triv = np.loadtxt(path_triv)

# Plot valley splitting vs. E
plot_vs(data_None, data_ff, data_triv, path_vs_E)