3. Position-Dependent Two-Valley \(\mathbf{k} \cdot \mathbf{p}\) Model
3.1. Requirements
3.1.1. Software components
QTCAD
NumPy
SciPy
Matplotlib
3.1.2. Python script
qtcad/examples/practical_application/Tunnel_Falls_DAPS/valley_kp_model.py
3.1.3. References
3.2. Briefing
This file defines the position-dependent two-valley electron
\(\mathbf{k} \cdot \mathbf{p}\) model used later in the detuning-spectrum
calculation. It is not intended to be run directly. Instead, it is imported by
Calculating the Energy Spectrum as a Function of Detuning, where the returned
ElectronKPModel is attached
to the double-dot subdevice before solving the single-electron Schrödinger
equation.
The two components of the envelope-function basis represent the \(+z\) and \(-z\) conduction-band valleys of silicon. These are the low-energy valleys selected by confinement along the heterostructure growth direction in silicon quantum wells [FCTC07, ZDM+13]. Sharp interfaces, alloy disorder, and atomic-scale roughness can mix the two valley states through a complex intervalley-coupling field \(\Delta(\mathbf r)\) [SCalderonC+11, SCalderonH+09]. Here, \(\Delta(\mathbf r)\) denotes the local valley coupling. The corresponding local valley splitting is denoted \(E_\mathrm{V}(\mathbf r)\).
In the \((+z,-z)\) valley basis, the local valley Hamiltonian is modeled as
The off-diagonal quantity \(\Delta(\mathbf r)\) is the coupling that mixes the two valley basis states. Diagonalizing this two-by-two matrix gives local eigenvalues \(\pm|\Delta(\mathbf r)|\), so the energy separation between the two valley eigenstates is
Thus, the valley coupling \(\Delta(\mathbf r)\) is the complex matrix element used by the electron \(\mathbf{k} \cdot \mathbf{p}\) model, whereas the valley splitting \(E_\mathrm{V}(\mathbf r)\) is the resulting energy separation. The splitting fixes the magnitude of the coupling through
but it does not by itself determine the complex phase of \(\Delta(\mathbf r)\).
The Tunnel Falls DAPS experiment motivates treating the valley coupling as a
position-dependent correlated random field across the dot array
[MEB+25]. The experimentally relevant input is the statistical
distribution of valley splittings, \(E_\mathrm{V}\). The functions defined
below use these splitting statistics to generate a spatially varying coupling
field \(\Delta(\mathbf r)\), whose magnitude is consistent with
\(E_\mathrm{V}(\mathbf r)=2|\Delta(\mathbf r)|\). This sampled coupling
field is then passed to
ElectronKPModel so that
valley mixing is included in the later Schrödinger calculations.
3.3. Header
We begin with the imports needed to sample a correlated valley-coupling field, interpolate it onto finite-element nodes, and define the two-valley \(\mathbf{k} \cdot \mathbf{p}\) model.
from pathlib import Path
import matplotlib.pyplot as plt
import numpy as np
from mpl_toolkits.axes_grid1 import make_axes_locatable
from scipy.interpolate import RegularGridInterpolator
from qtcad.device import ElectronKPModel
from qtcad.device import ElectronKPParameter
from qtcad.device import constants as ct
The RegularGridInterpolator
is used to map regularly sampled fields onto the unstructured finite-element nodes.
The QTCAD imports provide the objects and constants needed to build the
ElectronKPModel.
3.4. Constructing a valley-coupling map
The goal of this section is to generate one spatially correlated realization of the intervalley-coupling field \(\Delta(\mathbf r)\) over the relevant region of the Tunnel Falls device. The model parameters can be chosen from, or fit to, experimentally measured valley-splitting statistics [LEJ+23, MEB+25].
The important distinction is that the random field sampled by the code is the complex valley coupling \(\Delta(\mathbf r)\), whereas the local valley splitting is the derived quantity \(E_\mathrm{V}(\mathbf r)=2|\Delta(\mathbf r)|\).
Thus, the code samples \(\Delta\) first and then computes \(E_\mathrm{V}\) from its magnitude.
3.4.1. Gaussian covariance
The valley-coupling field is sampled at in-plane positions \(\boldsymbol{\rho}_j=(x_j, y_j)\). The random map has no explicit \(z\) dependence: the same value \(\Delta(x, y)\) is applied to all mesh nodes with the same \((x, y)\) coordinates. This is an effective description of in-plane variations, rather than a fully atomistic model of the interfaces.
Following the Tunnel Falls analysis of [MEB+25], the intervalley coupling is modeled as a zero-mean complex Gaussian random field. Writing the valley coupling at every position as
the real and imaginary parts, \(X_j\) and \(Y_j\), are sampled as independent Gaussian random fields with identical covariance matrices. The valley splitting is then obtained afterward from
For the spatial covariance, we use the Gaussian form motivated by Appendix D.1 of [LEJ+23]. For an \(N\)-point map, the covariance matrix used for either the real or imaginary components is
The real and imaginary components are taken to be uncorrelated, so that the real-imaginary cross covariances vanish. The parameters \(\ell_x\) and \(\ell_y\) are correlation lengths along the transport and transverse directions, respectively.
The function gaussian_covariance defined below computes this covariance matrix.
def gaussian_covariance(
positions: np.ndarray,
*,
sigma_delta_complex: float,
correlation_length_x: float,
correlation_length_y: float,
) -> np.ndarray:
"""Return the covariance for one component of the complex field Delta.
Args:
positions: Array of ``(x, y)`` positions where the field is sampled.
sigma_delta_complex: Square root of the variance of the complex
valley-coupling field, defined so that
``<|Delta|**2> = sigma_delta_complex**2``.
correlation_length_x: Correlation length along x.
correlation_length_y: Correlation length along y.
Returns:
Gaussian covariance matrix for either the real or imaginary part of
the complex valley-coupling field.
"""
if correlation_length_x <= 0.0:
raise ValueError("correlation_length_x must be positive.")
if correlation_length_y <= 0.0:
raise ValueError("correlation_length_y must be positive.")
dx = positions[:, None, 0] - positions[None, :, 0]
dy = positions[:, None, 1] - positions[None, :, 1]
r2 = (
dx * dx / correlation_length_x**2
+ dy * dy / correlation_length_y**2
)
return 0.5 * sigma_delta_complex**2 * np.exp(-r2)
The arrays dx and dy contain all pairwise coordinate differences
between sampled positions. The returned covariance matrix is dense because
each sampled point is correlated with every other sampled point, with the
correlation decaying as a Gaussian function of their in-plane separation.
3.4.2. Sampling the valley map
The function below draws one correlated realization of the complex intervalley-coupling field \(\Delta(\mathbf r)\). The experimental input is the mean valley splitting, \(\langle E_\mathrm{V}\rangle\), while the Gaussian covariance model is parameterized by the complex-coupling scale \(\sigma_\Delta\). Using \(\mathrm{Var}[\mathrm{Re}\,\Delta] =\mathrm{Var}[\mathrm{Im}\,\Delta]=\sigma_\Delta^2/2\) [LEJ+23], these quantities are related by [LEJ+23, MEB+25]
The function below generates a valley-coupling map and associated valley-splitting map
using the above relation and the covariance matrix output by the gaussian_covariance
function.
def sample_valley_coupling_map(
positions: np.ndarray,
*,
mean_valley_splitting: float,
correlation_length_x: float,
correlation_length_y: float,
seed: int,
) -> tuple[np.ndarray, np.ndarray]:
"""Sample a correlated valley-coupling field and its valley splitting.
Args:
positions: Array of ``(x, y)`` positions where the field is sampled.
mean_valley_splitting: Target mean valley splitting ``<E_V>``.
correlation_length_x: Correlation length along x.
correlation_length_y: Correlation length along y.
seed: Random seed used for reproducibility.
Returns:
Tuple containing the derived valley splitting ``E_V = 2*abs(Delta)``
and the sampled complex intervalley coupling ``Delta``.
"""
if mean_valley_splitting <= 0.0:
raise ValueError("mean_valley_splitting must be positive.")
# See Marcks et al., Nat. Commun. 16, 11381 (2025).
sigma_delta_complex = mean_valley_splitting / np.sqrt(np.pi)
# Get covariance matrix assuming Gaussian distribution.
covariance = gaussian_covariance(
positions,
sigma_delta_complex=sigma_delta_complex,
correlation_length_x=correlation_length_x,
correlation_length_y=correlation_length_y,
)
# Sample multivariate Gaussian to account for correlation lengths.
rng = np.random.default_rng(seed)
mean = np.zeros(covariance.shape[0])
delta_real = rng.multivariate_normal(mean=mean, cov=covariance)
delta_imag = rng.multivariate_normal(mean=mean, cov=covariance)
delta = delta_real + 1j * delta_imag
# Compute valley splitting from valley coupling.
valley_splitting = 2.0 * np.abs(delta)
return valley_splitting, delta
The random seed is exposed so that later detuning-spectrum calculations can
use a fixed disorder realization while varying only the plunger-gate detuning.
The function returns both quantities: the real-valued valley splitting
\(E_\mathrm{V}\), and the complex coupling field \(\Delta\), which is used to
define the two-valley
ElectronKPModel.
3.5. Interpolating to mesh nodes
The sampled valley-coupling map is defined on a regular two-dimensional grid.
To use it in the Schrödinger calculation, however, the coupling must be
evaluated at the nodes of a potentially unstructured finite-element
SubMesh. The next helper performs this
interpolation for the double-dot SubDevice.
def interpolate_map_to_nodes(
x_coords: np.ndarray,
y_coords: np.ndarray,
map_values: np.ndarray,
xy_nodes: np.ndarray,
) -> np.ndarray:
"""Interpolate a regular x-y map onto submesh nodes.
Args:
x_coords: x coordinates of the regular map.
y_coords: y coordinates of the regular map.
map_values: Values sampled on the regular ``(x, y)`` grid.
xy_nodes: Mesh-node coordinates where the map should be evaluated.
Returns:
Interpolated map values on the mesh nodes.
"""
interpolator = RegularGridInterpolator(
(x_coords, y_coords),
map_values,
method="linear",
bounds_error=False,
)
clipped_points = np.column_stack(
(
np.clip(xy_nodes[:, 0], x_coords[0], x_coords[-1]),
np.clip(xy_nodes[:, 1], y_coords[0], y_coords[-1]),
)
)
return interpolator(clipped_points)
The coordinates are clipped with np.clip before interpolation so that
points lying just outside the sampled grid due to numerical roundoff are projected back
onto the nearest grid boundary. This avoids out-of-bounds interpolation errors.
3.6. Plotting the valley map
Before constructing the final \(\mathbf{k} \cdot \mathbf{p}\) model, we add a helper function which saves a color plot of the sampled valley-splitting map.
def save_valley_map(
file_name: str | Path,
x_coords: np.ndarray,
y_coords: np.ndarray,
valley_splitting_grid: np.ndarray,
) -> None:
"""Save an image of the sampled valley-splitting map.
Args:
file_name: Output image path.
x_coords: x coordinates of the regular map.
y_coords: y coordinates of the regular map.
valley_splitting_grid: Valley splitting sampled on the regular
``(x, y)`` grid.
"""
figure, axis = plt.subplots(figsize=(5.5, 3.5))
color_mesh = axis.pcolormesh(
x_coords / 1e-9,
y_coords / 1e-9,
valley_splitting_grid.T / ct.e * 1e6,
shading="auto",
)
axis.set_aspect("equal")
axis.set_xlabel("x [nm]")
axis.set_ylabel("y [nm]")
divider = make_axes_locatable(axis)
colorbar_axis = divider.append_axes("right", size="5%", pad=0.08)
figure.colorbar(
color_mesh,
cax=colorbar_axis,
label=r"Valley splitting [$\mu$eV]",
)
figure.tight_layout()
figure.savefig(file_name, dpi=200)
plt.close(figure)
3.7. Generating the \(\mathbf{k} \cdot \mathbf{p}\) model
The final step is to generate the
ElectronKPModel that can be used to
model systems with spatially correlated valley couplings.
3.7.1. Effective-mass prefactors
The valley coupling defines how the two valley components mix. The kinetic
energy still comes from the usual anisotropic effective-mass Hamiltonian. In
QTCAD, the inverse electron effective-mass tensor is stored in the device
parameter Me_inv. The three transform functions below extract the
diagonal components used in
def prefactor_x(me_inv: np.ndarray) -> np.ndarray:
"""Return the x kinetic-energy prefactor.
Args:
me_inv: Inverse electron effective-mass tensor on the mesh nodes.
Returns:
Prefactor multiplying the ``k_x^2`` operator.
"""
return ct.hbar**2 / 2.0 * np.asarray(me_inv)[0, 0]
def prefactor_y(me_inv: np.ndarray) -> np.ndarray:
"""Return the y kinetic-energy prefactor.
Args:
me_inv: Inverse electron effective-mass tensor on the mesh nodes.
Returns:
Prefactor multiplying the ``k_y^2`` operator.
"""
return ct.hbar**2 / 2.0 * np.asarray(me_inv)[1, 1]
def prefactor_z(me_inv: np.ndarray) -> np.ndarray:
"""Return the z kinetic-energy prefactor.
Args:
me_inv: Inverse electron effective-mass tensor on the mesh nodes.
Returns:
Prefactor multiplying the ``k_z^2`` operator.
"""
return ct.hbar**2 / 2.0 * np.asarray(me_inv)[2, 2]
3.7.2. Two-valley \(\mathbf{k} \cdot \mathbf{p}\) model
We now construct the QTCAD
ElectronKPModel. The
constant term contains the local valley matrix
\(H_\mathrm{V}(\mathbf r)\). The word “constant” here means that this term
is \(k\)-independent; it may still be position-dependent.
def build_two_valley_model(off_diagonal: np.ndarray) -> ElectronKPModel:
"""Build an explicit two-valley electron k·p model.
Args:
off_diagonal: Complex intervalley coupling ``Delta`` evaluated on
the mesh nodes.
Returns:
Electron k·p model with anisotropic kinetic terms and local valley
coupling.
"""
delta_nodes = np.asarray(off_diagonal, dtype=np.complex128).reshape(-1)
valley_matrix = np.zeros((delta_nodes.size, 2, 2), dtype=np.complex128)
valley_matrix[:, 0, 1] = delta_nodes
valley_matrix[:, 1, 0] = np.conjugate(delta_nodes)
return ElectronKPModel(
bands=2,
constant=valley_matrix,
quadratic={
("x", "x"): ElectronKPParameter("Me_inv", transform=prefactor_x),
("y", "y"): ElectronKPParameter("Me_inv", transform=prefactor_y),
("z", "z"): ElectronKPParameter("Me_inv", transform=prefactor_z),
},
)
The array valley_matrix has shape (num_nodes, 2, 2). The first axis
indexes the mesh nodes, while the last two axes act on the valley basis.
The bands=2 argument tells the Schrödinger solver that the electron
envelope function has two components. The quadratic dictionary adds the
anisotropic effective-mass kinetic terms for the three Cartesian directions.
The final helper is the only function that later scripts need to call. It builds a regular two-dimensional grid, samples a correlated intervalley-coupling realization, interpolates \(\Delta\) onto the finite-element nodes, and returns the two-valley \(\mathbf{k} \cdot \mathbf{p}\) model.
def generate_correlated_two_valley_model(
xy_nodes: np.ndarray,
*,
mean_valley_splitting: float,
correlation_length_x: float,
correlation_length_y: float,
map_nx: int,
map_ny: int,
random_seed: int,
valley_map_filename: str | Path | None = None,
) -> ElectronKPModel:
"""Build the two-valley k·p model.
Args:
xy_nodes: Mesh-node coordinates where the k·p model is evaluated.
mean_valley_splitting: Mean value used to set the random-field scale.
correlation_length_x: Correlation length along x.
correlation_length_y: Correlation length along y.
map_nx: Number of x-grid points in the sampled map.
map_ny: Number of y-grid points in the sampled map.
random_seed: Random seed used for reproducibility.
valley_map_filename: Optional output path for a diagnostic plot of the
sampled valley-splitting map.
Returns:
Electron k·p model with a correlated intervalley coupling evaluated on
the mesh nodes.
"""
x_coords = np.linspace(
np.min(xy_nodes[:, 0]),
np.max(xy_nodes[:, 0]),
map_nx,
)
y_coords = np.linspace(
np.min(xy_nodes[:, 1]),
np.max(xy_nodes[:, 1]),
map_ny,
)
x_grid, y_grid = np.meshgrid(x_coords, y_coords, indexing="ij")
positions = np.stack((x_grid, y_grid), axis=-1).reshape((-1, 2))
# Sample the valley coupling over the device.
vs, delta = sample_valley_coupling_map(
positions,
mean_valley_splitting=mean_valley_splitting,
correlation_length_x=correlation_length_x,
correlation_length_y=correlation_length_y,
seed=random_seed,
)
if valley_map_filename is not None:
save_valley_map(
file_name=valley_map_filename,
x_coords=x_coords,
y_coords=y_coords,
valley_splitting_grid=vs.reshape(map_nx, map_ny),
)
# Interpolate onto mesh nodes.
delta_mesh = interpolate_map_to_nodes(
x_coords,
y_coords,
delta.reshape(map_nx, map_ny),
xy_nodes,
)
# Create the k·p model.
model = build_two_valley_model(delta_mesh)
return model
3.8. Full code
__copyright__ = "Copyright 2022-2026, Nanoacademic Technologies Inc."
from pathlib import Path
import matplotlib.pyplot as plt
import numpy as np
from mpl_toolkits.axes_grid1 import make_axes_locatable
from scipy.interpolate import RegularGridInterpolator
from qtcad.device import ElectronKPModel
from qtcad.device import ElectronKPParameter
from qtcad.device import constants as ct
def gaussian_covariance(
positions: np.ndarray,
*,
sigma_delta_complex: float,
correlation_length_x: float,
correlation_length_y: float,
) -> np.ndarray:
"""Return the covariance for one component of the complex field Delta.
Args:
positions: Array of ``(x, y)`` positions where the field is sampled.
sigma_delta_complex: Square root of the variance of the complex
valley-coupling field, defined so that
``<|Delta|**2> = sigma_delta_complex**2``.
correlation_length_x: Correlation length along x.
correlation_length_y: Correlation length along y.
Returns:
Gaussian covariance matrix for either the real or imaginary part of
the complex valley-coupling field.
"""
if correlation_length_x <= 0.0:
raise ValueError("correlation_length_x must be positive.")
if correlation_length_y <= 0.0:
raise ValueError("correlation_length_y must be positive.")
dx = positions[:, None, 0] - positions[None, :, 0]
dy = positions[:, None, 1] - positions[None, :, 1]
r2 = (
dx * dx / correlation_length_x**2
+ dy * dy / correlation_length_y**2
)
return 0.5 * sigma_delta_complex**2 * np.exp(-r2)
def sample_valley_coupling_map(
positions: np.ndarray,
*,
mean_valley_splitting: float,
correlation_length_x: float,
correlation_length_y: float,
seed: int,
) -> tuple[np.ndarray, np.ndarray]:
"""Sample a correlated valley-coupling field and its valley splitting.
Args:
positions: Array of ``(x, y)`` positions where the field is sampled.
mean_valley_splitting: Target mean valley splitting ``<E_V>``.
correlation_length_x: Correlation length along x.
correlation_length_y: Correlation length along y.
seed: Random seed used for reproducibility.
Returns:
Tuple containing the derived valley splitting ``E_V = 2*abs(Delta)``
and the sampled complex intervalley coupling ``Delta``.
"""
if mean_valley_splitting <= 0.0:
raise ValueError("mean_valley_splitting must be positive.")
# See Marcks et al., Nat. Commun. 16, 11381 (2025).
sigma_delta_complex = mean_valley_splitting / np.sqrt(np.pi)
# Get covariance matrix assuming Gaussian distribution.
covariance = gaussian_covariance(
positions,
sigma_delta_complex=sigma_delta_complex,
correlation_length_x=correlation_length_x,
correlation_length_y=correlation_length_y,
)
# Sample multivariate Gaussian to account for correlation lengths.
rng = np.random.default_rng(seed)
mean = np.zeros(covariance.shape[0])
delta_real = rng.multivariate_normal(mean=mean, cov=covariance)
delta_imag = rng.multivariate_normal(mean=mean, cov=covariance)
delta = delta_real + 1j * delta_imag
# Compute valley splitting from valley coupling.
valley_splitting = 2.0 * np.abs(delta)
return valley_splitting, delta
def interpolate_map_to_nodes(
x_coords: np.ndarray,
y_coords: np.ndarray,
map_values: np.ndarray,
xy_nodes: np.ndarray,
) -> np.ndarray:
"""Interpolate a regular x-y map onto submesh nodes.
Args:
x_coords: x coordinates of the regular map.
y_coords: y coordinates of the regular map.
map_values: Values sampled on the regular ``(x, y)`` grid.
xy_nodes: Mesh-node coordinates where the map should be evaluated.
Returns:
Interpolated map values on the mesh nodes.
"""
interpolator = RegularGridInterpolator(
(x_coords, y_coords),
map_values,
method="linear",
bounds_error=False,
)
clipped_points = np.column_stack(
(
np.clip(xy_nodes[:, 0], x_coords[0], x_coords[-1]),
np.clip(xy_nodes[:, 1], y_coords[0], y_coords[-1]),
)
)
return interpolator(clipped_points)
def save_valley_map(
file_name: str | Path,
x_coords: np.ndarray,
y_coords: np.ndarray,
valley_splitting_grid: np.ndarray,
) -> None:
"""Save an image of the sampled valley-splitting map.
Args:
file_name: Output image path.
x_coords: x coordinates of the regular map.
y_coords: y coordinates of the regular map.
valley_splitting_grid: Valley splitting sampled on the regular
``(x, y)`` grid.
"""
figure, axis = plt.subplots(figsize=(5.5, 3.5))
color_mesh = axis.pcolormesh(
x_coords / 1e-9,
y_coords / 1e-9,
valley_splitting_grid.T / ct.e * 1e6,
shading="auto",
)
axis.set_aspect("equal")
axis.set_xlabel("x [nm]")
axis.set_ylabel("y [nm]")
divider = make_axes_locatable(axis)
colorbar_axis = divider.append_axes("right", size="5%", pad=0.08)
figure.colorbar(
color_mesh,
cax=colorbar_axis,
label=r"Valley splitting [$\mu$eV]",
)
figure.tight_layout()
figure.savefig(file_name, dpi=200)
plt.close(figure)
def prefactor_x(me_inv: np.ndarray) -> np.ndarray:
"""Return the x kinetic-energy prefactor.
Args:
me_inv: Inverse electron effective-mass tensor on the mesh nodes.
Returns:
Prefactor multiplying the ``k_x^2`` operator.
"""
return ct.hbar**2 / 2.0 * np.asarray(me_inv)[0, 0]
def prefactor_y(me_inv: np.ndarray) -> np.ndarray:
"""Return the y kinetic-energy prefactor.
Args:
me_inv: Inverse electron effective-mass tensor on the mesh nodes.
Returns:
Prefactor multiplying the ``k_y^2`` operator.
"""
return ct.hbar**2 / 2.0 * np.asarray(me_inv)[1, 1]
def prefactor_z(me_inv: np.ndarray) -> np.ndarray:
"""Return the z kinetic-energy prefactor.
Args:
me_inv: Inverse electron effective-mass tensor on the mesh nodes.
Returns:
Prefactor multiplying the ``k_z^2`` operator.
"""
return ct.hbar**2 / 2.0 * np.asarray(me_inv)[2, 2]
def build_two_valley_model(off_diagonal: np.ndarray) -> ElectronKPModel:
"""Build an explicit two-valley electron k·p model.
Args:
off_diagonal: Complex intervalley coupling ``Delta`` evaluated on
the mesh nodes.
Returns:
Electron k·p model with anisotropic kinetic terms and local valley
coupling.
"""
delta_nodes = np.asarray(off_diagonal, dtype=np.complex128).reshape(-1)
valley_matrix = np.zeros((delta_nodes.size, 2, 2), dtype=np.complex128)
valley_matrix[:, 0, 1] = delta_nodes
valley_matrix[:, 1, 0] = np.conjugate(delta_nodes)
return ElectronKPModel(
bands=2,
constant=valley_matrix,
quadratic={
("x", "x"): ElectronKPParameter("Me_inv", transform=prefactor_x),
("y", "y"): ElectronKPParameter("Me_inv", transform=prefactor_y),
("z", "z"): ElectronKPParameter("Me_inv", transform=prefactor_z),
},
)
def generate_correlated_two_valley_model(
xy_nodes: np.ndarray,
*,
mean_valley_splitting: float,
correlation_length_x: float,
correlation_length_y: float,
map_nx: int,
map_ny: int,
random_seed: int,
valley_map_filename: str | Path | None = None,
) -> ElectronKPModel:
"""Build the two-valley k·p model.
Args:
xy_nodes: Mesh-node coordinates where the k·p model is evaluated.
mean_valley_splitting: Mean value used to set the random-field scale.
correlation_length_x: Correlation length along x.
correlation_length_y: Correlation length along y.
map_nx: Number of x-grid points in the sampled map.
map_ny: Number of y-grid points in the sampled map.
random_seed: Random seed used for reproducibility.
valley_map_filename: Optional output path for a diagnostic plot of the
sampled valley-splitting map.
Returns:
Electron k·p model with a correlated intervalley coupling evaluated on
the mesh nodes.
"""
x_coords = np.linspace(
np.min(xy_nodes[:, 0]),
np.max(xy_nodes[:, 0]),
map_nx,
)
y_coords = np.linspace(
np.min(xy_nodes[:, 1]),
np.max(xy_nodes[:, 1]),
map_ny,
)
x_grid, y_grid = np.meshgrid(x_coords, y_coords, indexing="ij")
positions = np.stack((x_grid, y_grid), axis=-1).reshape((-1, 2))
# Sample the valley coupling over the device.
vs, delta = sample_valley_coupling_map(
positions,
mean_valley_splitting=mean_valley_splitting,
correlation_length_x=correlation_length_x,
correlation_length_y=correlation_length_y,
seed=random_seed,
)
if valley_map_filename is not None:
save_valley_map(
file_name=valley_map_filename,
x_coords=x_coords,
y_coords=y_coords,
valley_splitting_grid=vs.reshape(map_nx, map_ny),
)
# Interpolate onto mesh nodes.
delta_mesh = interpolate_map_to_nodes(
x_coords,
y_coords,
delta.reshape(map_nx, map_ny),
xy_nodes,
)
# Create the k·p model.
model = build_two_valley_model(delta_mesh)
return model