4. Computing \(g\)-tensors in an FD-SOI device

4.1. Requirements

4.1.1. Software components

  • QTCAD

  • Paraview

4.1.2. Python script

  • qtcad/examples/tutorials/g_tensor.py

4.1.3. References

4.2. Briefing

In this tutorial, we will return to the fully-depleted silicon-on-insulator (FD-SOI) device studied in A double quantum dot device in a fully-depleted silicon-on-insulator transistor and explore the magnetic properties of confined electrons using QTCAD’s tight-binding (TB) approach. In particular, we will demonstrate how to compute the \(g\)-tensor of a confined electron and show how the spin splitting can depend on the direction of an applied in-plane magnetic field. The analysis presented is inspired by the work of Cifuentes et al. from Ref. [CTG+24].

Note

This tutorial involves a TB simulation of \(\sim 250,000\) atoms including spin effects. We expect the simulations presented here to run in less than (but on the order of) an hour on a relatively advanced desktop workstation (as described in Software requirements).

4.4. Constructing the atomic structure and the TB Hamiltonian

We start by setting up a Mesh object, which defines a mesh over a geometric model of the FD-SOI device of interest (see A double quantum dot device in a fully-depleted silicon-on-insulator transistor for more details). This Mesh is used to create a Device. We load into this Device the potential obtained in Exchange coupling in a double quantum dot in FD-SOI—Part 1: Perturbation theory which led to bound states localized at each quantum dot.

# Setup
script_dir          = pathlib.Path(__file__).parent.resolve()
out_dir             = script_dir / "output"

## Load the mesh
path_mesh           = script_dir / "meshes" / "refined_dqdfdsoi.msh"
scaling             = 1e-9
mesh                = Mesh(scaling, path_mesh)

## Create device and set potential in device to result of FEM simulation
d                   = Device(mesh, conf_carriers="e")
path_phi_hdf5       = out_dir / "detuned_potential.hdf5"
d.set_potential(io.load(path_phi_hdf5))

Next, we construct the atomic structure of a relevant region of the device. As explained in Unit cell tiling, this is done by tiling a unit cell so as to fill a specified region of space. By default, the unit cell being tiled is the conventional unit cell of bulk silicon. We choose to tile this unit cell over the region below the left plunger gate of the FD-SOI device of A double quantum dot device in a fully-depleted silicon-on-insulator transistor. The specific dimensions of the chosen region are justified a posteriori by looking at the confinement of the eigenfunctions of the system (see Fig. 4.6.1, below). Because the eigenfunctions do not extend to the boundaries of the atomistically-modeled region, we can be sure that the chosen region is large enough to capture the relevant physics of the system. In particular, in this configuration, the buried oxide layer (see A double quantum dot device in a fully-depleted silicon-on-insulator transistor) is not necessary for electron confinement. Because the system will not be relaxed, this buried oxide layer cannot affect the atomic structure near the locations of the eigenfunctions of interest. However, in calculations involving relaxation, we expect that the inclusion of the buried oxide layer can lead to important modifications of the presented results.

Model of the FD-SOI device studied in the tutorial.

Fig. 4.4.1 Model of the FD-SOI device studied in this tutorial. The gates B1, B2, and B3 are barrier gates, while the gates P1 and P2 are the left and right plunger gates, respectively. The region below the left plunger gate (indicated with the red dashed line) is modeled atomistically. From the results of Exchange coupling in a double quantum dot in FD-SOI—Part 1: Perturbation theory, this region is expected to contain bound states.

# Dimensions - Box below the left plunger gate
A                   = 1e-10 # Angstroms
x                   = np.array([-200.0, 200.0]) * A
y                   = np.array([-250.0, -100.0]) * A
z                   = np.array([-80.0, 10.0]) * A

# Construct the atomic structure of the heterostructure
atoms = Atoms(
    x = x,
    y = y,
    z = z,
    verbose=True
)

Next, we specify a RoughSurface object to model the interface between the two layers of the heterostructure (see Multiscale simulations of a SiGe–Si–SiGe quantum dot—Part 2: Tight-binding model for valley splitting for details on how to create a rough surface). We specify the \(\mathrm{SiO}_2\) layer using the new_layer method of the atoms object.

Note

In QTCAD, silicon oxide is modeled within the virtual crystal approximation, wherein fictitious \(\mathrm{SiO}_2\) “atoms” are arranged on a diamond lattice with the same lattice constant as that of bulk silicon. See Unit cell tiling for more details.

# Define the material stack of the heterostructure
## Create rough surface describing interface between Si and SiO2 layers
rng_seed = 1000
interface = RoughSurface(hurst=0.3, mean=0.0, rms=2.0*A, rng_seed=rng_seed)
## Create the heterostructure layer
SiO2 = np.array([["SiO2", 1.0]])
atoms.new_layer(z_top=10.0*A, z_bot=interface, atom_species=SiO2)

Below, we show the atomic structure being simulated with color coding for the different atoms: silicon atoms are shown as red spheres, and the fictitious \(\mathrm{SiO}_2\) “atoms” are shown as blue spheres.

Atomic structure of the FD-SOI device.

Fig. 4.4.2 Atomic structure of the FD-SOI device visualized with Paraview. The silicon atoms are shown as red spheres, while the fictitious \(\mathrm{SiO}_2\) “atoms” are shown as blue spheres. The file containing this data is created below via the save_vtu function.

Finally, we set additional attributes of atoms which specify how the structure is modeled. We set the potential over the atomic structure using that of the Device object d using the set_potential_from_device method. Because we aim to understand the magnetic/spin properties of the electrons confined to the modeled device, we include spin–orbit coupling in the model by setting the soc attribute via the set_soc method. We also use the B attribute to specify the magnetic field, which is set to a small value of \(10^{-6}\;\mathrm{T}\) in the \(z\) direction, and is used, for our purposes, to break the time-reversal symmetry of the system and lift the spin degeneracy. When magnetic fields are applied to the system, by default, the only magnetic effect that is included in the computation is the coupling of the long-range (device-scale) orbital magnetic moments of the electrons to the applied field ; this coupling is handled using the Peierls substitution, as detailed in Magnetic orbital effects—Peierls substitution, and affects the offsite blocks of the TB Hamiltonian. The set_Zeeman method is used to specify whether the Zeeman effect at the atomic scale is to be included in the computation; this effect alters the onsite block of the TB Hamiltonian, as detailed in Zeeman effect. Finally, because we are interested in the spin properties of the confined electrons, our total basis size can be quite large; it is two times larger than a basis in the absence of spin. To reduce the computational cost of the simulation, we use an \(sp3s^{\ast}\) basis (which has ten basis states per atom when spin is included). This basis is a good compromise between accuracy and computational cost. It is enforced by setting the params_name attribute using the set_tb_params method to Klimeck_cb and Kim.

atoms.set_potential_from_device(d)
atoms.set_soc(True)
atoms.set_Bfield(np.array([0.0, 0.0, 1e-6]))
atoms.set_Zeeman(True)
atoms.set_tb_params(("Klimeck_cb", "Kim"))

4.5. Adding strain

In the previous tutorial, the structure was relaxed using the Keating valence force-field model so that strain effects were accounted for in the model. Here, instead, we will specify the strain tensor applied to the atomic structure by defining a function strain_func which outputs the \(3 \times 3\) strain tensor for every value of its (three) inputs, the \(x\), \(y\), and \(z\) coordinates of a point in space, and passing this function to the strain_atomic_structure method of the Atoms object.

To model the strain tensor, we will use a simple formulation that accounts for the difference in thermal expansion of the \(\mathrm{SiO}_2\) layer compared to the silicon channel. We will assume that the strain is in the plane of the heterostructure and is due entirely to the cooling of the device from the growth temperature (assumed to be \(T_i \sim 1000 \, \mathrm{K}\), the rough temperature of a furnace used to grow \(\text{SiO}_2\)) to the operating temperature of the device (\(T_f \sim 0 \, \mathrm{K}\)). Moreover, we will assume that the system experiences no shear strain and that the normal (in-plane) strain is isotropic. In the oxide layer we will take the magnitude of the strain to be

\[\epsilon_{\parallel} = \left(\alpha_{\mathrm{SiO}_2} - \alpha_{\mathrm{Si}} \right) \Delta T,\]

where \(\alpha_{m}\) is the coefficient of thermal expansion of material \(m\), and \(\Delta T = T_f - T_i\). In the silicon layer, the strain is in the opposite direction. Indeed, if the \(\mathrm{SiO}_2\) layer contracts relative to the silicon channel, the silicon channel expands relative to the \(\mathrm{SiO}_2\) layer and vice-versa. Moreover, we will assume that the strain decays exponentially with depth into the silicon layer as we expect the strain to be less and less important the further we get from the interface. We will assume a decay length of \(\lambda = 2 \, \mathrm{nm}\). All-in-all, the non-vanishing component of the strain tensor are taken to be

\[\begin{split}\epsilon_{xx}(z) = \epsilon_{yy}(z) &=\left\{ \begin{array}{ll} \epsilon_{\parallel} &\textrm{if}\;z > 0,\\ -\epsilon_{\parallel} e^{z/\lambda}&\textrm{if}\;z \le 0. \end{array} \right.\end{split}\]

The mathematical formalism used to strain the atomic structure from an input strain tensor is described in Displacing atoms according to a strain tensor field.

# Define strain function
def strain_func(x, y, z):
    alpha_SiO2  = 0.55e-6
    alpha_Si    = 2.6e-6
    delta_T     = 1000.0 - 0.0 # K
    epsilon_parallel = (alpha_SiO2 - alpha_Si) * delta_T
    lambda_ = 2.0 * 1e-9 # 2 nm
    if z > 0:
        epsilon = epsilon_parallel
    else:
        epsilon = -epsilon_parallel * np.exp(z / lambda_)

    return np.array([
        [epsilon, 0.0, 0.0],
        [0.0, epsilon, 0.0],
        [0.0, 0.0, 0.0]
    ])

# Set the strain
atoms.strain_atomic_structure(strain_func)

4.6. Solving the Schrödinger equation

Once we have set up the Atoms object, we set up the TB Schrödinger solver. To do so, we first define the SolverParams object, which contains the parameters of the solver. We specify the number of eigenpairs to be computed (num_states) and whether messages pertaining to the solver are printed to the screen (verbose).

# Solve the Schrödinger equation

## Parametrize the Schrödinger solver
solver_params               = SolverParams()
solver_params.num_states    = 10
solver_params.verbose       = True
solver = Solver(atoms, solver_params=solver_params)

# Solve the Schrödinger equation
energy_target = 0.04 * ct.e
solver.solve(energy_target)
atoms.print_energies()

After parametrizing the solver, we solve the Schrödinger equation by calling the solve method of the Solver object, specifying that we are looking for solutions near energy_target. The value of this parameter is chosen by looking at the energy levels from the finite-element simulations of Exchange coupling in a double quantum dot in FD-SOI—Part 1: Perturbation theory. There, the ground-state energy was found to be \(\sim 0.04 \, \mathrm{eV}\), i.e. near the Fermi energy \(E_F=0\). Accordingly, we set energy_target = 0.04 * ct.e in the code above. Once the solver has finished running, we can print the eigenenergies nearest energy_target using the print_energies method and find relatively good agreement with the energies of Exchange coupling in a double quantum dot in FD-SOI—Part 1: Perturbation theory.

Atomistic tight-binding eigenenergies (eV):
[0.04836001 0.04836001 0.04843321 0.04843321 0.05010970 0.05010970 0.05022183 0.05022183 0.05576065 0.05576065]

We also save the probability density of each state over the atoms of the system to a .vtu file using the save_vtu function.

## Save modulus-square of projection of states on atoms to a .vtu file
psi = atoms.eigenfunctions
out_dict = {f"state {i}": psi[i].get_prob_on_atoms() for i in range(len(psi))}
path_vtu = out_dir / "FD-SOI_psi.vtu"
save_vtu(atoms=atoms, out_dict=out_dict, path=path_vtu, verbose=True)

Note

The save_vtu function is used to save the probability density of the states over the atoms of the system to a .vtu file, which can be visualized using Paraview. This function also saves the chemical species of each atom in the system (e.g. \(\mathrm{Si}\), \(\mathrm{SiO}_2\), etc.) to the .vtu file.

The ground state of the system is shown in the figure below

Ground state of the FD-SOI device.

Fig. 4.6.1 Ground state of the FD-SOI device. The color coding indicates the probability density of the ground state over the atoms of the system. The file containing this data is created above via the save_vtu function.

4.7. \(g\)-tensor calculation

QTCAD is also equipped with a \(g\)-tensor Solver which can be used to compute the effective \(g\)-tensor for a user-specified two-level system. To summarize what this solver does: it projects the terms in the TB Hamiltonian linear in the applied magnetic field onto the two-level system of interest. From there, it parametrizes the resulting \(2 \times 2\) Hamiltonian matrix in terms of a pseudospin operator \(\mathbf{s}\) (a vector of four \(2\times 2\) matrices, Eq. (6.9), whose \(0^{\mathrm{th}}\) component is \(I/2\), where \(I\) is the identity matrix) and an effective \(g\)-tensor (see Magnetic effects for more details). The \(g\)-tensor contains information about how the pseudospin is influenced by a magnetic field (to first order in the magnetic field). This solver is applied using

# Compute the g-tensor
states          = (psi[1], psi[0])
g_tensor_solver = GTensorSolver(atoms, states)
g               = g_tensor_solver.solve()

print("g-tensor")
print(f"g = {g}")

Once the \(g\)-tensor is computed, we can visualize its features by computing the associated Zeeman splittings as a function of the direction of the applied magnetic field. The effective \(g\)-tensor is used to parametrize the effective Zeeman Hamiltonian in a two-level subspace. As such, it depends only on the basis chosen for this subspace and is magnetic-field independent. However, the projected Zeeman Hamiltonian does depend on the (direction and magnitude of the) applied magnetic field. Consequently, the splitting \(\Delta\) between the eigenenergies of the Zeeman Hamiltonian will also depend on the (direction and magnitude of the) applied magnetic field. We investigate how \(\Delta\) depends on the direction of the applied in-plane magnetic field using the get_Zeeman_splitting method of the Solver class:

## Zeeman splitting for in-plane magnetic field
Theta = np.linspace(0, np.pi, 41) # angles in radians
Delta = []
for theta in Theta:
    B0 = 1 # units: T
    # Vary direction of in-plane field
    B = B0*np.array([np.cos(theta), np.sin(theta), 0])
    Delta.append(g_tensor_solver.get_Zeeman_splitting(B))

Delta = np.array(Delta) / (ct.h * 1e6) # Convert to MHz / T
Delta0 = np.mean(Delta)
print("Average spin-splitting: " + str(Delta0) + " MHz / T")

The resulting array Delta contains the Zeeman splitting for each angle theta. The values of Delta are computed by diagonalizing the projected Zeeman Hamiltonian, which is a \(2 \times 2\) matrix (Eq. (6.8)), for different directions of the applied in-plane magnetic field, whose magnitude is fixed to \(B=1 \, \mathrm{T}\). Because the Zeeman Hamiltonian is linear in the applied magnetic field, we can convert the units of the splittings to MHz/T by multiplying by \(\frac{e}{10^6 h}\), where \(-e\) is the electron charge and \(h\) is Planck’s constant. We also compute the average Zeeman splitting Delta0.

We now plot the Zeeman splitting relative to the average Zeeman splitting as a function of the in-plane angle of the magnetic field \(\theta\):

## Plotting the spin splitting
fig, ax = plt.subplots()
ax.plot(Theta / np.pi, (Delta - Delta0))
ax.set_xlabel(r"$\theta$ [$\pi$]")
ax.set_ylabel(r"$\Delta - \Delta_0$ [MHz / T]")
ax.grid()
fig.savefig(out_dir / "in_plane_g.png", dpi=300)

## Save spin splitting
with open(out_dir / "in_plane_g.txt", "a") as f:
    np.savetxt(f, Theta[np.newaxis, :] / np.pi, fmt="%.6f")
    np.savetxt(f, Delta[np.newaxis, :], fmt="%.6f")

The resulting plot is shown below:

Zeeman splitting as a function of the in-plane angle of the magnetic field.

Fig. 4.7.4 Zeeman splitting as a function of the in-plane angle of the magnetic field. The average Zeeman splitting is subtracted from the Zeeman splitting to emphasize magnitude of the variations on the scale of the axes.

This plot can be compared to that of Fig. 3 (e) of Ref. [CTG+24]. In [CTG+24] dots with quasi-cylindrical symmetry are investigated: in the absence of strain and interface effects, the in-plane confinement potential is isotropic. As such, the Zeeman splitting for a magnetic field along the \(x\) or \(y\) direction is the same. Here, the cylindrical symmetry is broken by the geometry of the device which is elongated in the \(y\) direction (compared to the \(x\) direction). Accordingly, we do not expect to have the same Zeeman splitting for a magnetic field along the \(x\) or \(y\) direction. Moreover, the amplitude of variations in the Zeeman splitting is on the order of \(\sim 2 \, \mathrm{MHz / T}\), roughly an order of magnitude smaller than that presented in Fig. 3 (e) of Ref. [CTG+24]. However, the quantum dots studied in Ref. [CTG+24] are much smaller (diameter of \(\sim 5-10 \, \mathrm{nm}\)) than those analyzed here (\(\sim 15 \, \mathrm{nm} \times 40 \, \mathrm{nm}\)) which could explain the difference in amplitude of the Zeeman splitting.

Note

While the analysis of the \(g\)-tensor presented here is performed using the \(g\)-tensor Solver, users can also compute the \(g\)-tensor and other magnetic/spin properties using the full TB Schrödinger Solver. Concretely, we can solve the full TB Hamiltonian for different values of an applied magnetic field. The downside of this approach is that it requires diagonalizing the full TB Hamiltonian for each value of the applied magnetic field, which can be computationally expensive.

4.8. Full code

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

import pathlib
import numpy as np
import matplotlib.pyplot as plt
import qtcad.device.constants as ct
from qtcad.device.mesh3d import Mesh
from qtcad.device import Device
from qtcad.device import io
from qtcad.atoms.atoms import Atoms
from qtcad.atoms.rough_surface import RoughSurface
from qtcad.atoms.schrodinger import Solver, SolverParams
from qtcad.atoms.analysis import save_vtu
from qtcad.atoms.g_tensor import Solver as GTensorSolver

# Setup
script_dir          = pathlib.Path(__file__).parent.resolve()
out_dir             = script_dir / "output"

## Load the mesh
path_mesh           = script_dir / "meshes" / "refined_dqdfdsoi.msh"
scaling             = 1e-9
mesh                = Mesh(scaling, path_mesh)

## Create device and set potential in device to result of FEM simulation
d                   = Device(mesh, conf_carriers="e")
path_phi_hdf5       = out_dir / "detuned_potential.hdf5"
d.set_potential(io.load(path_phi_hdf5))

# Dimensions - Box below the left plunger gate
A                   = 1e-10 # Angstroms 
x                   = np.array([-200.0, 200.0]) * A
y                   = np.array([-250.0, -100.0]) * A
z                   = np.array([-80.0, 10.0]) * A 

# Construct the atomic structure of the heterostructure
atoms = Atoms(
    x = x,
    y = y,
    z = z,
    verbose=True
)

# Define the material stack of the heterostructure
## Create rough surface describing interface between Si and SiO2 layers
rng_seed = 1000
interface = RoughSurface(hurst=0.3, mean=0.0, rms=2.0*A, rng_seed=rng_seed)
## Create the heterostructure layer
SiO2 = np.array([["SiO2", 1.0]])
atoms.new_layer(z_top=10.0*A, z_bot=interface, atom_species=SiO2) 

atoms.set_potential_from_device(d)
atoms.set_soc(True)
atoms.set_Bfield(np.array([0.0, 0.0, 1e-6]))
atoms.set_Zeeman(True)
atoms.set_tb_params(("Klimeck_cb", "Kim"))


# Define strain function
def strain_func(x, y, z):
    alpha_SiO2  = 0.55e-6
    alpha_Si    = 2.6e-6
    delta_T     = 1000.0 - 0.0 # K
    epsilon_parallel = (alpha_SiO2 - alpha_Si) * delta_T
    lambda_ = 2.0 * 1e-9 # 2 nm
    if z > 0:
        epsilon = epsilon_parallel
    else:
        epsilon = -epsilon_parallel * np.exp(z / lambda_)

    return np.array([
        [epsilon, 0.0, 0.0],
        [0.0, epsilon, 0.0],
        [0.0, 0.0, 0.0]
    ])

# Set the strain
atoms.strain_atomic_structure(strain_func)

# Solve the Schrödinger equation

## Parametrize the Schrödinger solver
solver_params               = SolverParams()
solver_params.num_states    = 10
solver_params.verbose       = True
solver = Solver(atoms, solver_params=solver_params)

# Solve the Schrödinger equation
energy_target = 0.04 * ct.e
solver.solve(energy_target)
atoms.print_energies()

## Save modulus-square of projection of states on atoms to a .vtu file
psi = atoms.eigenfunctions
out_dict = {f"state {i}": psi[i].get_prob_on_atoms() for i in range(len(psi))}
path_vtu = out_dir / "FD-SOI_psi.vtu"
save_vtu(atoms=atoms, out_dict=out_dict, path=path_vtu, verbose=True)

# Compute the g-tensor
states          = (psi[1], psi[0])
g_tensor_solver = GTensorSolver(atoms, states)
g               = g_tensor_solver.solve()

print("g-tensor")
print(f"g = {g}")

## Zeeman splitting for in-plane magnetic field
Theta = np.linspace(0, np.pi, 41) # angles in radians
Delta = []
for theta in Theta:
    B0 = 1 # units: T
    # Vary direction of in-plane field
    B = B0*np.array([np.cos(theta), np.sin(theta), 0]) 
    Delta.append(g_tensor_solver.get_Zeeman_splitting(B))

Delta = np.array(Delta) / (ct.h * 1e6) # Convert to MHz / T
Delta0 = np.mean(Delta)
print("Average spin-splitting: " + str(Delta0) + " MHz / T")

## Plotting the spin splitting
fig, ax = plt.subplots()
ax.plot(Theta / np.pi, (Delta - Delta0))
ax.set_xlabel(r"$\theta$ [$\pi$]")
ax.set_ylabel(r"$\Delta - \Delta_0$ [MHz / T]")
ax.grid()
fig.savefig(out_dir / "in_plane_g.png", dpi=300)

## Save spin splitting
with open(out_dir / "in_plane_g.txt", "a") as f:
    np.savetxt(f, Theta[np.newaxis, :] / np.pi, fmt="%.6f")
    np.savetxt(f, Delta[np.newaxis, :], fmt="%.6f")