Scattering Systems

Nanokwant supports scattering matrix calculations for 1D systems with attached leads.

Basic Usage

import numpy as np
from nanokwant import scattering_system, compute_smatrix

# Define your system
system = {
    0: {"mu": np.eye(2)},
    1: {"t": np.eye(2)},
}

params = {"mu": 0.5, "t": 1.0}
num_sites = 50
energy = 0.5

# Compute the scattering matrix
S = compute_smatrix(system, num_sites, params, energy, leads="both")

Lead Configuration

You can attach leads at different positions:

• leads="both": Attach leads at both ends (default)
• leads="left": Attach lead only at the left end
• leads="right": Attach lead only at the right end

# Single lead example
S_left = compute_smatrix(system, num_sites, params, energy, leads="left")

Advanced: Direct Access to Linear System

For more control, you can access the linear system directly. The helper returns the banded matrix, bandwidth, right-hand side, and the Kwant propagating modes for each attached lead:

from scipy.linalg import solve_banded

# Build the linear problem
lhs, rhs, modes = scattering_system(system, num_sites, params, energy)
l, u = lhs.bandwidth

# Solve once using the combined RHS (all incoming modes)
solution = solve_banded((l, u), np.asarray(lhs), rhs)

# Assemble the scattering matrix manually from the solved amplitudes
lead_order = [name for name in ("left", "right") if name in modes]
nmodes = {name: len(modes[name].velocities) // 2 for name in lead_order}
iface_dim = {name: modes[name].wave_functions.shape[0] for name in lead_order}

block_offsets = {}
if "left" in modes:
    block_offsets["left"] = 0
if "right" in modes:
    block_offsets["right"] = solution.shape[0] - iface_dim["right"]

rows = [
    solution[block_offsets[name] : block_offsets[name] + nmodes[name]]
    for name in lead_order
]
S_manual = np.vstack(rows) if rows else np.zeros((0, 0))

# Matches the helper result
S_direct = compute_smatrix(system, num_sites, params, energy, leads="both")
assert np.allclose(S_manual, S_direct)

Position-Dependent Parameters

Scattering systems support position-dependent parameters:

def barrier(x):
    """Potential barrier in the middle."""
    return np.where((x >= 20) & (x <= 30), 2.0, 0.5)

params = {"mu": barrier, "t": 1.0}
S = compute_smatrix(system, num_sites, params, energy, leads="both")

Evanescent Modes

The implementation properly handles evanescent modes using Kwant's stabilized mode computation:

# System with different bandwidths for different orbitals
system = {
    0: {"mu": np.eye(2)},
    1: {"t": np.diag([1.0, 0.1])},  # Different hopping strengths
}

# At certain energies, some modes will be evanescent
S = compute_smatrix(system, num_sites, params, energy=0.3, leads="both")

Validation

The scattering matrix should be unitary for Hermitian systems:

# Check unitarity
unitarity_error = np.max(np.abs(S @ S.conj().T - np.eye(S.shape[0])))
print(f"Unitarity error: {unitarity_error:.2e}")

Unitarity error: 1.78e-15