Lecture 24.4

Long-range strain — boundary conditions and dot configurations

In the lecture

Four strain boundary conditions for an In${0.6}$Ga${0.4}$As dome (15 nm radius, 5.4 nm height) as the buffer grows: free surfaces underestimate the strain and the electron energy, fixed surfaces overestimate them with the opposite trend; growing only the vertical buffer reproduces the free result, a hint that the vertical strain dominates. The strain reaches about 16 nm.
Two realistic scenarios for a flat InAs dot on a wetting layer: a deeply buried dot (all faces fixed to GaAs) and a lateral array of dots under a finite cap (bottom fixed, sides periodic, top free).
With separate strain (up to 64 million atoms) and electronic (2 million atoms) domains, the mixed boundary conditions converge for much smaller systems: the free top lets the strain relax vertically.
A strain domain extended vertically converges faster than one extended laterally.
In dot arrays, the cap thickness and the distance between dots shift the levels by a few meV, not always monotonically; strain creates the dots, reaches far, further vertically, which is why dots stack, while the wavefunctions stay confined.

In this notebook

You compute, with open-nemo's finite-element elasticity solver and the effective-mass model on the strained local band edges of the previous lectures: the electron level of the alloyed dome under free, vertically extended free, periodic and fixed boundaries; the S and P electron and the hole levels of a flat InAs dot under fixed and mixed boundaries for strain domains up to 160 nm wide; laterally against vertically extended domains; and the ground state of a dot array as a function of the cap thickness and the distance between dots.

Run the notebookEvery number and figure on this page comes from it. Change it, rerun it, take it further.

Download .ipynb View on GitHub

Run it on your own machine

The notebook needs open-nemo (Python ≥ 3.10 and a C++20 compiler; BLAS/LAPACK).

git clone https://github.com/denphi/open-nemo.git && cd open-nemo
python3 -m venv .venv && source .venv/bin/activate
pip install -e ".[notebooks]"
jupyter lab From_Atoms_to_Devices/10_long_range_strain/24_4_boundary_conditions_and_dot_configurations.ipynb

Code13 lines · cell 1

import time

import numpy as np
import matplotlib.pyplot as plt
import opennemo as nemo
from opennemo.models import kp
from opennemo.solvers import elasticity as el
from opennemo.solvers.poisson import Grid
from opennemo.structures.shapes import Box, Dome
from opennemo.visualization import style

style.use(notebook=True)
T_START = time.time()

Code24 lines · cell 2

GaAs = kp.material("GaAs", temperature=300)            # Vurgaftman et al. (2001), 300 K
InAs = kp.material("InAs", temperature=300)
E0 = GaAs.valence_band_offset                           # energy zero: unstrained GaAs valence-band top


def stiffness(m):
    """Cubic elastic stiffness (GPa, Voigt) of a k·p material."""
    return el.cubic_stiffness(m.c11, m.c12, m.c44)


def misfit(m):
    """Eigenstrain of m on the GaAs reference lattice: a_m / a_GaAs - 1."""
    return m.lattice_constant / GaAs.lattice_constant - 1


for m in (GaAs, InAs):
    print(f"{m.name}: a = {m.lattice_constant:.5f} nm, E_g = {m.band_gap:.3f} eV, "
          f"E_c = {m.valence_band_offset + m.band_gap - E0:+.3f} eV, "
          f"E_v = {m.valence_band_offset - E0:+.3f} eV, "
          f"a_c = {m.ac:+.2f} eV, m_e = {m.electron_mass:.3f} m0")
print(f"InAs misfit on GaAs: {100 * misfit(InAs):.2f} %")

InGaAs = kp.material("InGaAs", x=0.6, temperature=300)
print(f"{InGaAs.name}: a = {InGaAs.lattice_constant:.5f} nm, misfit on GaAs {100 * misfit(InGaAs):.2f} %")

GaAs: a = 0.56532 nm, E_g = 1.422 eV, E_c = +1.422 eV, E_v = +0.000 eV, a_c = -7.17 eV, m_e = 0.067 m0
InAs: a = 0.60583 nm, E_g = 0.354 eV, E_c = +0.564 eV, E_v = +0.210 eV, a_c = -5.08 eV, m_e = 0.026 m0
InAs misfit on GaAs: 7.16 %
InGaAs(x=0.6): a = 0.58963 nm, misfit on GaAs 4.30 %

Code44 lines · cell 3

VOIGT = [(0, 0), (1, 1), (2, 2), (1, 2), (0, 2), (0, 1)]
_LINEAR = {}


def local_edges(m, eps):
    """Local band edges (eV, from E0) of material m under the strains eps (n, 3, 3): the zone-centre
    eigenvalues of its 8-band k·p Hamiltonian (Bir-Pikus strain terms). Returns E_c, E_HH, E_LH; the
    heavy holes are the valence pair with the least |Z> character (z: growth axis)."""
    if m.name not in _LINEAR:                            # H(eps) = H0 + sum eps_ij D_ij: linear in eps
        h0 = kp.bulk_hamiltonian(m, np.zeros(3), 8)
        d = []
        for i, j in VOIGT:
            u = np.zeros((3, 3))
            u[i, j] = u[j, i] = 1.0
            d.append(kp.bulk_hamiltonian(m, np.zeros(3), 8, strain=u) - h0)
        _LINEAR[m.name] = h0, np.array(d)
    h0, d = _LINEAR[m.name]
    v = np.stack([eps[:, i, j] for i, j in VOIGT], -1)
    w, u = np.linalg.eigh(h0[None] + np.einsum("na,aij->nij", v, d))
    z = np.abs(u[:, 3, 2:6]) ** 2 + np.abs(u[:, 7, 2:6]) ** 2   # basis S, X, Y, Z per spin
    order = np.argsort(z, axis=1)
    top = w[:, 2:6]
    e_hh = np.take_along_axis(top, order[:, :2], 1).max(1)
    e_lh = np.take_along_axis(top, order[:, 2:], 1).max(1)
    return w[:, 6] - E0, e_hh - E0, e_lh - E0


def region_at(result, points):
    """Region index (0: the GaAs matrix, >0: the dot) of the finite-element cell containing each
    point: the cell whose material local_strain_at measures the strain against."""
    p = result.problem
    ids = []
    for d, a in enumerate(p.grid.axes):
        x = points[:, d]
        if p.periodic[d]:
            x = a[0] + np.mod(x - a[0], a[-1] - a[0])
        ids.append(np.clip(np.searchsorted(a, x, side="right") - 1, 0, len(a) - 2))
    return p.material[np.ravel_multi_index(tuple(ids), p.cell_shape)]


eps_test = np.array([[-0.05, 0.01, 0.0], [0.01, -0.05, 0.0], [0.0, 0.0, 0.04]])
direct = np.linalg.eigvalsh(kp.bulk_hamiltonian(InAs, np.zeros(3), 8, strain=eps_test)) - E0
print("direct 8-band:  E_c, top valence =", np.round(direct[[6, 5]], 5))
print("local_edges:    E_c, E_HH, E_LH  =", np.round(np.ravel(local_edges(InAs, eps_test[None])), 5))

direct 8-band:  E_c, top valence = [0.86859 0.32936]
local_edges:    E_c, E_HH, E_LH  = [0.86859 0.32936 0.1134 ]

Code82 lines · cell 4

HB2M = 0.0380998                                         # hbar^2 / 2 m0 (eV nm^2)


def bdd_hamiltonian(shape, h, inv_mass, potential):
    """Single-band BenDaniel-Duke operator -div (hbar^2 / 2m) grad + V on a box grid of spacings h
    (nm) with hard walls; inv_mass (n, 3): m0/m along x, y, z at each point."""
    import scipy.sparse as sp
    n = int(np.prod(shape))
    idx = np.arange(n).reshape(shape)
    diag = np.asarray(potential, dtype=float).copy()
    rows, cols, vals = [], [], []
    for ax in range(3):
        a = np.take(idx, range(shape[ax] - 1), axis=ax).ravel()
        b = np.take(idx, range(1, shape[ax]), axis=ax).ravel()
        t = HB2M * 0.5 * (inv_mass[a, ax] + inv_mass[b, ax]) / h[ax] ** 2
        rows += [a, b]
        cols += [b, a]
        vals += [-t, -t]
        np.add.at(diag, a, t)
        np.add.at(diag, b, t)
        for wall in (np.take(idx, [0], axis=ax).ravel(), np.take(idx, [shape[ax] - 1], axis=ax).ravel()):
            diag[wall] += HB2M * inv_mass[wall, ax] / h[ax] ** 2   # psi = 0 one step outside
    off = sp.coo_matrix((np.concatenate(vals), (np.concatenate(rows), np.concatenate(cols))), shape=(n, n))
    return off.tocsr() + sp.diags(diag)


def hole_inv_mass(m):
    """Heavy-hole m0/m in the plane (gamma1 + gamma2) and along z (gamma1 - 2 gamma2)."""
    return np.array([m.gamma1 + m.gamma2, m.gamma1 + m.gamma2, m.gamma1 - 2 * m.gamma2])


def dot_levels(result, dot, lo, hi, h=(1.0, 1.0, 0.5), num=1, holes=True):
    """Lowest electron and highest hole levels (eV, from E0) of a dot of material `dot` in GaAs, in
    the hard-wall box lo..hi (nm) on a grid of spacings h, with the strain and the regions of the
    finite-element `result`. Returns a dict with the levels, the grid and the band edges on it."""
    axes = [np.arange(lo[i] + h[i] / 2, hi[i], h[i]) for i in range(3)]
    shape = tuple(len(a) for a in axes)
    pts = np.stack(np.meshgrid(*axes, indexing="ij"), -1).reshape(-1, 3)
    ind = region_at(result, pts) > 0
    eps = result.local_strain_at(pts)
    ec, ev = np.empty(len(pts)), np.empty(len(pts))
    for m, sel in ((GaAs, ~ind), (dot, ind)):
        c, hh, lh = local_edges(m, eps[sel])
        ec[sel], ev[sel] = c, np.maximum(hh, lh)
    inv_me = 1.0 / np.where(ind, dot.electron_mass, GaAs.electron_mass)
    re = nemo.eigenstates(bdd_hamiltonian(shape, h, np.repeat(inv_me[:, None], 3, 1), ec), num=num)
    out = {"e": np.sort(re["energies"]), "e_states": re, "points": pts, "shape": shape, "axes": axes,
           "ec": ec, "ev": ev, "inside": ind}
    if holes:                                            # holes: lowest states of (kinetic - E_v)
        inv_mh = np.where(ind[:, None], hole_inv_mass(dot)[None], hole_inv_mass(GaAs)[None])
        rh = nemo.eigenstates(bdd_hamiltonian(shape, h, inv_mh, -ev), num=num)
        out["h"], out["h_states"] = -np.sort(rh["energies"]), rh
    return out


def graded_axis(lo, hi, core_lo, core_hi, h, h_max=2.0, growth=1.25):
    """Grid nodes from lo to hi (nm): spacing h on [core_lo, core_hi] around the dot, growing by
    `growth` per cell outside it up to h_max. Elasticity has no length scale, so the far buffer
    needs no fine grid."""
    def grow(start, end, sign):
        pts, x, s = [], start, h
        while True:
            s = min(s * growth, h_max)
            x += sign * s
            if (x >= end - 0.3 * s) if sign > 0 else (x <= end + 0.3 * s):
                return pts + [end]
            pts.append(x)
    core = list(np.arange(core_lo, core_hi + 1e-9, h))
    return np.array(grow(core_lo, lo, -1)[::-1] + core + grow(core_hi, hi, +1))


def boundary_nodes(prob, faces):
    """Nodes on the named faces ('x', 'y': both lateral faces, 'bottom', 'top') of the grid."""
    nodes, axes = prob.grid.nodes, prob.grid.axes
    mask = np.zeros(len(nodes), dtype=bool)
    for f in faces:
        if f in ("x", "y"):
            d = "xy".index(f)
            mask |= np.isclose(nodes[:, d], axes[d][0]) | np.isclose(nodes[:, d], axes[d][-1])
        else:
            mask |= np.isclose(nodes[:, 2], axes[2][0 if f == "bottom" else -1])
    return mask

1. Four boundary conditions for one dot#

The free outer surfaces of Lectures 24.2 and 24.3 let the dot expand the whole crystal. What happens with other choices? The lecture compares four for a dome of In${0.6}$Ga${0.4}$As (15 nm radius, 5.4 nm height) with a growing GaAs buffer:

free: no constraint on the outer surfaces;
free*: free surfaces, but only the vertical buffer grows; the lateral buffer stays at 4 nm;
periodic: the box repeats in all three directions (a three-dimensional array of dots);
fixed: the outer surfaces are pinned to the GaAs lattice.

Slide 3: electron ground state against the buffer for free, free\ (vertical growth only), periodic and fixed boundary conditions: free underestimates and fixed overestimates the strain and the energy, with opposite trends; free* stays close to free.*

open-nemo's periodic axes keep the period of the GaAs reference lattice (the solver's documentation: "the period stays fixed, pseudomorphic to the reference"), so the periodic case here is an array of dots whose supercell cannot expand. The electronic box reaches 4 nm beyond the dot in every direction, enough by Lecture 24.3.

Code35 lines · cell 6

R, H = 15.0, 5.4
dome = Dome((0, 0, 0), (R, R, H))
lo, hi = (-R - 4, -R - 4, -4), (R + 4, R + 4, H + 4)


def solve_dome(bc, b):
    """The alloyed dome in GaAs with buffer b (nm) under boundary condition bc."""
    bl = 4.0 if bc == "free*" else b                                 # lateral buffer
    xa = graded_axis(-R - bl, R + bl, -R - 1, R + 1, 1.0)
    za = graded_axis(-b, H + b, -1, H + 1, 0.5)
    prob = el.Elasticity.from_regions(Grid(xa, xa, za), [(None, stiffness(GaAs), 0.0),
                                                         (dome, stiffness(InGaAs), misfit(InGaAs))],
                                      periodic=(bc == "periodic",) * 3)
    if bc == "fixed":
        prob.fix(boundary_nodes(prob, ("x", "y", "bottom", "top")))
    return prob.solve()


t0 = time.time()
BUF = [4, 8, 12, 16]
BCS = ["free", "free*", "periodic", "fixed"]
e_bc = {bc: np.array([dot_levels(solve_dome(bc, b), InGaAs, lo, hi, holes=False)["e"][0] for b in BUF])
        for bc in BCS}
for bc in BCS:
    print(f"{bc:9s}: e1 = " + ", ".join(f"{e:.4f}" for e in e_bc[bc]) + f" eV for b = {BUF} nm; "
          f"last step {1e3 * (e_bc[bc][-1] - e_bc[bc][-2]):+.1f} meV")
print(f"{len(BCS) * len(BUF)} strain fields and levels in {time.time() - t0:.1f} s")

fig, ax = plt.subplots(figsize=(6.5, 3.8))
for k, bc in enumerate(BCS):
    ax.plot(BUF, e_bc[bc], marker="osd^"[k], color=style.sweep(4)[k], label=bc)
ax.set(xlabel="GaAs buffer b (nm)", ylabel="electron ground state e1 (eV)",
       title="In$_{0.6}$Ga$_{0.4}$As dome, four boundary conditions")
ax.legend(frameon=False)
fig.tight_layout()

free     : e1 = 1.2349, 1.2889, 1.3117, 1.3222 eV for b = [4, 8, 12, 16] nm; last step +10.5 meV
free*    : e1 = 1.2349, 1.2852, 1.3057, 1.3122 eV for b = [4, 8, 12, 16] nm; last step +6.5 meV
periodic : e1 = 1.3918, 1.3559, 1.3457, 1.3417 eV for b = [4, 8, 12, 16] nm; last step -4.0 meV
fixed    : e1 = 1.4451, 1.3822, 1.3589, 1.3488 eV for b = [4, 8, 12, 16] nm; last step -10.1 meV
16 strain fields and levels in 20.6 s

Figure computed by the notebook with open-nemo

The four curves split into two families, as on the slide:

Free surfaces give the lowest energies, rising by 87 meV from the 4 nm to the 16 nm buffer and still by 10 meV in the last step: the dot can push the small block apart, so it is too little compressed, and the strain and the energy are underestimated.
Fixed surfaces give the highest energies, falling by 96 meV, again by 10 meV in the last step: the GaAs frame around the dot cannot give way, and the strain is overestimated. The two families approach the answer from opposite sides and are still 27 meV apart at a 16 nm buffer; the strain reaches at least that far.
Free* grows only the vertical buffer and keeps 4 nm laterally, yet it reproduces 89 % of the free curve's rise (1.312 against 1.322 eV at 16 nm): the vertical extension of the buffer matters most.
Periodic with a fixed GaAs period behaves like the fixed case, between it and the free one: the supercell of the dot array cannot expand on average, so close neighbours compress each other. The lecture's periodic curve instead rises like the free one at a higher energy; that needs a supercell that can relax its period, which open-nemo's elasticity solver does not offer.

2. Two realistic scenarios: a buried dot and a dot array#

Real dots are not surrounded by free space. The lecture's second half uses the flatter dots seen in experiments, an InAs dome of 18.1 nm diameter and 1.7 nm height on a 0.6 nm InAs wetting layer, and two realistic sets of boundary conditions:

fixed: a deeply buried single dot, all six faces of the domain fixed to the GaAs lattice;
mixed: a lateral array of dots under a finite cap, the bottom (substrate) fixed to GaAs, the sides periodic, the top surface free.

Slide 4: fixed and mixed boundary conditions with strain and electronic domains Slide 4: the deeply buried single dot (fixed boundary conditions) and the lateral array of dots with a finite cap (mixed boundary conditions); the electronic domain (inner box) is a part of the strain domain (outer box) of width d and height h.

The strain domain is $d$ wide and $h = d/2$ high, the dot at its centre; the electronic domain stays a fixed box around the dot (4 nm beyond it on every side) with a grid of 1 nm laterally and 0.25 nm vertically. The finite-element grid is 1 nm laterally and 0.2 nm vertically near the dot, growing to 6 nm cells in the buffer.

Code31 lines · cell 9

D, HD, WL = 18.1, 1.7, 0.6
flat = Dome((0, 0, 0), (D / 2, D / 2, HD)) | Box((-1e3, -1e3, -WL), (1e3, 1e3, 0.0))   # dot + wetting layer
flat_lo, flat_hi = (-D / 2 - 4, -D / 2 - 4, -WL - 4), (D / 2 + 4, D / 2 + 4, HD + 4)


def solve_flat(bc, d, z_lo, z_hi):
    """Strain of the flat InAs dot in a GaAs domain d wide, from z_lo to z_hi (nm): 'fixed' (all faces
    on the GaAs lattice) or 'mixed' (bottom fixed, sides periodic, top free)."""
    xa = graded_axis(-d / 2, d / 2, -D / 2 - 1, D / 2 + 1, 1.0, h_max=6.0, growth=1.4)
    za = graded_axis(z_lo, z_hi, -WL - 1, HD + 1, 0.2, h_max=6.0, growth=1.4)
    prob = el.Elasticity.from_regions(Grid(xa, xa, za), [(None, stiffness(GaAs), 0.0),
                                                         (flat, stiffness(InAs), misfit(InAs))],
                                      periodic=(bc == "mixed",) * 2 + (False,))
    prob.fix(boundary_nodes(prob, ("x", "y", "bottom", "top") if bc == "fixed" else ("bottom",)))
    return prob.solve()


def flat_levels(res, holes=True):
    return dot_levels(res, InAs, flat_lo, flat_hi, h=(1.0, 1.0, 0.25), num=3, holes=holes)


t0 = time.time()
DS = [40, 80, 120]
scen = {(bc, d): flat_levels(solve_flat(bc, d, -d / 4, d / 4)) for bc in ("fixed", "mixed") for d in DS}
t_scen = time.time() - t0
n_el = 8 * np.prod(np.subtract(flat_hi, flat_lo)) / 0.56532**3
print(f"electronic domain {np.subtract(flat_hi, flat_lo).round(1)} nm: {n_el / 1e6:.2f} million atoms")
for (bc, d), lv in scen.items():
    print(f"{bc:5s} d = {d:3d} nm ({8 * d * d * d / 2 / 0.56532**3 / 1e6:5.1f} M atoms): "
          f"e = {np.round(lv['e'], 4)}, h = {np.round(lv['h'], 4)} eV")
print(f"{len(scen)} cases in {t_scen:.1f} s")

electronic domain [26.1 26.1 10.3] nm: 0.31 million atoms
fixed d =  40 nm (  1.4 M atoms): e = [1.3459 1.4202 1.4207], h = [0.1881 0.1248 0.1246] eV
fixed d =  80 nm ( 11.3 M atoms): e = [1.32   1.3947 1.3952], h = [0.1995 0.1344 0.1342] eV
fixed d = 120 nm ( 38.3 M atoms): e = [1.3135 1.3882 1.3887], h = [0.2024 0.137  0.1368] eV
mixed d =  40 nm (  1.4 M atoms): e = [1.2961 1.3715 1.372 ], h = [0.2113 0.1445 0.1443] eV
mixed d =  80 nm ( 11.3 M atoms): e = [1.3016 1.3761 1.3766], h = [0.2078 0.142  0.1418] eV
mixed d = 120 nm ( 38.3 M atoms): e = [1.3022 1.3767 1.3772], h = [0.2073 0.1416 0.1414] eV
6 cases in 16.7 s

Code13 lines · cell 10

fig, axes = plt.subplots(1, 2, figsize=(11, 3.8))
for k, bc in enumerate(("fixed", "mixed")):
    col = style.sweep(2)[k]
    for ax, key in ((axes[0], "e"), (axes[1], "h")):
        v = np.array([scen[bc, d][key] for d in DS])
        ax.plot(DS, v[:, 0], marker="o", ls="-", color=col, label=f"S, {bc}")
        ax.plot(DS, v[:, 1:].mean(1), marker="s", ls="--", color=col, label=f"P, {bc}")
axes[0].set(xlabel="strain domain width d (nm), h = d/2", ylabel="electron energy (eV)",
            title="electron states")
axes[1].set(xlabel="strain domain width d (nm), h = d/2", ylabel="hole energy (eV)", title="hole states")
for ax in axes:
    ax.legend(frameon=False, fontsize=8, ncol=2)
fig.tight_layout()

Under fixed boundaries the S electron level falls from 1.346 eV (d = 40 nm, 1.4 million atoms in the strain domain) to 1.314 eV (120 nm, 38 million atoms) and still moves by 6.5 meV in the last step; the hole S level rises from 0.188 to 0.202 eV. Under mixed boundaries the S electron is 1.296, 1.302 and 1.302 eV: converged to 0.6 meV at d = 80 nm, and the holes likewise within 0.5 meV. The P electron states, a nearly degenerate pair, follow the S state 75 meV higher. The S levels lie in the ranges of the slide (electrons 1.30–1.35 eV, holes 0.19–0.21 eV), from a much smaller electronic domain (0.31 million atoms instead of 2 million); the single-band model puts the P levels further from S than tight binding does, the hole P pair most of all (63 meV). The free top of the mixed domain lets the strain relax vertically, which is easier than pushing a fixed frame, so the mixed domain converges for much smaller systems.

Slide 5: S and P electron and hole states against the domain width for fixed and mixed boundaries Slide 5: electron and hole S and P states of the flat InAs dot against the strain domain width d (d/h = 2, electronic domain 2 million atoms, strain domain up to 64 million atoms): the mixed boundary conditions converge for much smaller systems than the fixed ones.

3. Lateral or vertical: which extension matters?#

Lecture 24.2 suggested that the strain reaches further along the growth direction. The lecture tests it with two shapes of the strain domain under mixed boundary conditions: A, laterally dominated ($h = d/2$, the domains of §2), and B, vertically dominated ($h = 2d$).

Slide 6: laterally and vertically dominated strain domains Slide 6: the S and P electron states in a laterally dominated (A, h = d/2) and a vertically dominated (B, h = 2d) strain domain: B converges faster than A.

Code15 lines · cell 12

t0 = time.time()
DB = [40, 80, 120]
vert = {d: flat_levels(solve_flat("mixed", d, -d, d), holes=False) for d in DB}
print(f"B (h = 2d)  d = {DB} nm: S = " + ", ".join(f"{vert[d]['e'][0]:.4f}" for d in DB) + " eV")
print(f"A (h = d/2) d = {DS} nm: S = " + ", ".join(f"{scen['mixed', d]['e'][0]:.4f}" for d in DS) + " eV")
print(f"({time.time() - t0:.1f} s)")
fig, ax = plt.subplots(figsize=(6.5, 3.8))
for k, (lab, ds, lv) in enumerate((("A: h = d/2", DS, [scen["mixed", d] for d in DS]),
                                   ("B: h = 2d", DB, [vert[d] for d in DB]))):
    ax.plot(ds, [x["e"][0] for x in lv], marker="o", ls="-", color=style.sweep(2)[k], label=f"S, {lab}")
    ax.plot(ds, [x["e"][1:].mean() for x in lv], marker="s", ls="--", color=style.sweep(2)[k],
            label=f"P, {lab}")
ax.set(xlabel="strain domain width d (nm)", ylabel="electron energy (eV)", title="mixed boundary conditions")
ax.legend(frameon=False, fontsize=8, ncol=2)
fig.tight_layout()

B (h = 2d)  d = [40, 80, 120] nm: S = 1.3007, 1.3023, 1.3024 eV
A (h = d/2) d = [40, 80, 120] nm: S = 1.2961, 1.3016, 1.3022 eV
(9.9 s)

At the same width $d$ the vertically dominated domain B is closer to convergence: 1.7 meV below its d = 120 nm value at d = 40 nm, where A is 6.1 meV below; at d = 80 nm, B is within 0.1 meV and A within 0.6 meV. Counting atoms instead of width the advantage shrinks (B at d = 40 nm already holds 5.7 million atoms, A at 80 nm 11.3 million), but the conclusion is the lecture's: these flat dots have more vertical than lateral strain influence, and the strain domain should be extended along the growth direction first.

4. Arrays of dots: capping layer and distance between dots#

With mixed boundary conditions the periodic sides describe a lateral array of identical dots, and the free top a cap of finite thickness. The lecture varies both: the substrate below the wetting layer is large (here 40 nm), the cap $c$ is 11.87, 17.52 or 23.18 nm, and the edge-to-edge distance between neighbouring dots $b$ goes from 10 to 30 nm (the period is $D + b$; the lecture goes to 35 nm).

Slide 7: ground state of a dot array against cap thickness and dot distance Slide 7: electron ground state of a lateral dot array for three cap thicknesses c against the edge-to-edge distance b between the dots: a few meV of variation, not monotonic for the thinnest cap.

Code16 lines · cell 14

t0 = time.time()
CAPS, GAPS = [11.87, 17.52, 23.18], [10, 15, 20, 30]
array_e = {c: np.array([flat_levels(solve_flat("mixed", D + b, -WL - 40.0, c), holes=False)["e"][0]
                        for b in GAPS]) for c in CAPS}
for c in CAPS:
    print(f"cap {c:5.2f} nm: e1 = " + ", ".join(f"{e:.5f}" for e in array_e[c]) + f" eV for b = {GAPS} nm "
          f"(range {1e3 * np.ptp(array_e[c]):.1f} meV, maximum at b = {GAPS[int(np.argmax(array_e[c]))]} nm)")
print(f"({time.time() - t0:.1f} s)")
fig, ax = plt.subplots(figsize=(6.5, 3.8))
for k, c in enumerate(CAPS):
    ax.plot(GAPS, array_e[c], marker="os^"[k], color=style.sweep(3)[k], label=f"cap c = {c} nm")
ax.set(xlabel="edge-to-edge distance between dots b (nm)", ylabel="electron ground state (eV)",
       title="lateral dot array, mixed boundary conditions")
ax.legend(frameon=False, fontsize=8)
fig.tight_layout()
print(f"notebook run time: {time.time() - T_START:.0f} s")

cap 11.87 nm: e1 = 1.29609, 1.29749, 1.29771, 1.29743 eV for b = [10, 15, 20, 30] nm (range 1.6 meV, maximum at b = 20 nm)
cap 17.52 nm: e1 = 1.29687, 1.29913, 1.30016, 1.30082 eV for b = [10, 15, 20, 30] nm (range 4.0 meV, maximum at b = 30 nm)
cap 23.18 nm: e1 = 1.29688, 1.29918, 1.30035, 1.30139 eV for b = [10, 15, 20, 30] nm (range 4.5 meV, maximum at b = 30 nm)
(21.6 s)
notebook run time: 69 s

The levels of the array move by only a few meV (1.6 to 4.5 meV) over these geometries, but in an instructive way. For the two thicker caps the electron level rises monotonically with the distance between the dots: close neighbours stretch the GaAs between them, the dots are less compressed, and their conduction band is lower. For the thinnest cap (11.87 nm) the level rises up to b = 20 nm and then falls again by 0.3 meV: the strain also relaxes through the nearby free top surface, and the two relaxation paths compete. The lecture finds the same non-monotonic curve for its thinnest cap at 1.310–1.316 eV; the model's levels are lower by about 15 meV, its trends the same.

Summary#

The strain boundary conditions decide the trend: free surfaces underestimate the strain and the electron energy and approach the result from below, fixed surfaces overestimate them and approach it from above. Extending only the vertical buffer gives nearly the free result.
For a flat InAs dot on a wetting layer, mixed boundary conditions (substrate fixed, periodic sides, free top) converge in much smaller domains than fixed ones: the free top lets the strain relax.
A vertically extended strain domain converges faster than a laterally extended one.
In a dot array the cap thickness and the distance between dots shift the ground state by a few meV, non-monotonically for a thin cap, where strain relaxes through the surface.
The lecture's summary: strain creates the dots, it is long-ranged and reaches further vertically, so dots stack on top of each other, while the wavefunctions stay confined to a few nanometres around the dot.

Exercises#

Add hole levels to §1 (holes=True). Do the four boundary conditions order the hole levels the same way as the electron levels?
Stack two flat dots vertically, 10 nm apart, in the mixed domain of §2 (flat shifted by 10 nm along $z$, combined with |). How does the strain under the upper dot differ from that under the lower one?
Repeat §4 with a 4 nm cap. Is the dependence on the distance between dots stronger, and why?