SciCode++
Switch benchSciCode++
Request dataset samplesby

Scaled tight binding model for a two dimensional electron gas at the (001) LaAlO₃/SrTiO₃ interfaceMaterial Science

27 verifiers

TaxonomyMaterial ScienceSemiconductor Materials

Task prompt

  • Gemini 3.1 Pro Preview62.5%of 8 runs
  • GPT-5.6 (sol)87.5%of 8 runs

Prompt

Compute the band structure of the (001)(001) LaAlO3_3/SrTiO3_3 two-dimensional electron gas along a specified k-path, then build the corresponding scaled real-space tight-binding Hamiltonian on a finite periodic lattice and verify that both representations yield identical low-energy spectra at scaling factor s=1s=1. The system has 6 internal degrees of freedom per site: 3 t2gt_{2g} orbitals (dxyd_{xy}, dxzd_{xz}, dyzd_{yz}) ×\times 2 spin channels (↑\uparrow, ↓\downarrow), with basis ordering (dxy ⁣↑,dxy ⁣↓,dxz ⁣↑,dxz ⁣↓,dyz ⁣↑,dyz ⁣↓)(d_{xy}\!\uparrow, d_{xy}\!\downarrow, d_{xz}\!\uparrow, d_{xz}\!\downarrow, d_{yz}\!\uparrow, d_{yz}\!\downarrow). The k-space Hamiltonian at each (kx,ky)(k_x, k_y) is H^(k)=H^0(k)+H^SO+H^RSO(k)\hat{H}(\mathbf{k}) = \hat{H}_0(\mathbf{k}) + \hat{H}_{SO} + \hat{H}_{RSO}(\mathbf{k}), where H^0\hat{H}_0 encodes anisotropic nearest-neighbour hopping with light (tlt_l) and heavy (tht_h) bands, crystal-field splitting ΔE\Delta E lowering dxyd_{xy}, and dxzd_{xz}/dyzd_{yz} hybridisation via tdt_d; H^SO\hat{H}_{SO} is the atomic spin-orbit coupling from the t2gt_{2g}-projected L⋅S\mathbf{L}\cdot\mathbf{S} interaction; and H^RSO\hat{H}_{RSO} is the orbital Rashba coupling arising from broken inversion symmetry at the interface. The real-space Hamiltonian on an Nx×NyN_x \times N_y lattice with periodic boundary conditions includes on-site blocks, nearest-neighbour hopping blocks H^x\hat{H}_x and H^y\hat{H}_y, and next-nearest-neighbour diagonal hopping blocks H^mix\hat{H}_{mix}. For scaling factor ss, kinetic hopping parameters scale as 1/s21/s^2, Rashba coupling scales as 1/s1/s, and atomic SOC is completely unscaled. The function must: (1) compute the k-space band structure along the given k-path, (2) build and diagonalise the real-space Hamiltonian on a periodic Nx×NyN_x \times N_y lattice at the specified scaling factor ss, (3) also compute k-space eigenvalues at all commensurate k-points of that lattice, and (4) return a tuple containing: the band structure array (shape (Nk,6)(N_k, 6), energies in meV), the sorted real-space eigenvalue array (length 6NxNy6 N_x N_y, in meV), and the maximum absolute deviation (in meV) between the sorted real-space spectrum and the sorted k-space spectrum evaluated at the commensurate k-points. Implement the following function:

def compute_lao_sto_band_and_realspace(k_path, tl, th, td, dE, dSO, dRSO, a, Nx, Ny, s):
    '''Compute the LAO/STO 2DEG band structure and real-space spectrum.

    Parameters:
      k_path : ndarray, shape (N_k, 2)
          Array of (kx, ky) values in units of 1/nm along the desired path.
      tl : float
          Light-band hopping parameter in meV.
      th : float
          Heavy-band hopping parameter in meV.
      td : float
          Hybridisation (next-nearest-neighbour) hopping in meV.
      dE : float
          Crystal-field splitting in meV (dxy lowered by this amount).
      dSO : float
          Atomic spin-orbit coupling strength in meV.
      dRSO : float
          Rashba spin-orbit coupling strength in meV.
      a : float
          Lattice constant in nm.
      Nx : int
          Number of sites in x-direction for real-space lattice.
      Ny : int
          Number of sites in y-direction for real-space lattice.
      s : float
          Scaling factor for the real-space mesh.

    Returns:
      bands : ndarray, shape (N_k, 6)
          Sorted eigenvalues at each k-point in meV.
      rs_eigs : ndarray, shape (6*Nx*Ny,)
          Sorted eigenvalues from real-space Hamiltonian in meV.
      max_dev : float
          Maximum absolute deviation between sorted k-space and real-space spectra.
    '''
    return bands, rs_eigs, max_dev

Background

The (001)(001) LaAlO3_3/SrTiO3_3 interface hosts a two-dimensional electron gas whose low-energy physics is governed by the t2gt_{2g} orbitals (dxyd_{xy}, dxzd_{xz}, dyzd_{yz}) of the Ti atoms on a square lattice with lattice constant a=0.394a = 0.394 nm. The Hamiltonian combines anisotropic tight-binding hopping (light hopping tl=875t_l = 875 meV, heavy hopping th=40t_h = 40 meV, hybridisation td=40t_d = 40 meV), a crystal-field splitting ΔE=47\Delta E = 47 meV that lowers the dxyd_{xy} band, atomic spin-orbit coupling ΔSO\Delta_{SO} from the t2gt_{2g}-projected L⋅S\mathbf{L}\cdot\mathbf{S} operator, and Rashba spin-orbit coupling ΔRSO\Delta_{RSO} from broken inversion symmetry. The scaled tight-binding approach coarsens the real-space mesh by a factor ss while preserving the low-energy physics through differential rescaling of the Hamiltonian parameters: kinetic terms scale as 1/s21/s^2, Rashba coupling as 1/s1/s, and atomic SOC remains unscaled. Agreement between the k-space and real-space spectra at s=1s=1 serves as a stringent consistency check on the Hamiltonian construction.

Parameters and returns

Compute the LAO/STO 2DEG band structure and real-space spectrum. Parameters: k_path : ndarray, shape (N_k, 2) Array of (kx, ky) values in units of 1/nm along the desired path. tl : float Light-band hopping parameter in meV. th : float Heavy-band hopping parameter in meV. td : float Hybridisation (next-nearest-neighbour) hopping in meV. dE : float Crystal-field splitting in meV (dxy lowered by this amount). dSO : float Atomic spin-orbit coupling strength in meV. dRSO : float Rashba spin-orbit coupling strength in meV. a : float Lattice constant in nm. Nx : int Number of sites in x-direction for real-space lattice. Ny : int Number of sites in y-direction for real-space lattice. s : float Scaling factor for the real-space mesh. Returns: bands : ndarray, shape (N_k, 6) Sorted eigenvalues at each k-point in meV. rs_eigs : ndarray, shape (6NxNy,) Sorted eigenvalues from real-space Hamiltonian in meV. max_dev : float Maximum absolute deviation between sorted k-space and real-space spectra.

Required dependencies

import numpy as np  # numpy==2.4.4
from scipy.linalg import eigvalsh  # scipy==1.17.1
from scipy.sparse import lil_matrix  # scipy==1.17.1
from scipy.sparse import lil_matrix, csr_matrix  # scipy==1.17.1
from scipy.sparse.linalg import eigsh  # scipy==1.17.1
import pytest  # pytest==9.0.2

Published by the benchmark's authors. Opening these shows how the problem is solved and what the graders assert.

Test case 1 of 5
import numpy as np
from scipy.sparse import lil_matrix, csr_matrix
from scipy.sparse.linalg import eigsh
from scipy.linalg import eigvalsh


def test_case_1():
    tl, th, td, dE, dSO, dRSO = 875.0, 40.0, 40.0, 47.0, 10.0, 5.0
    a = 0.394
    Nx, Ny, s = 2, 2, 1.0
    k_path = np.array([[0.0, 0.0], [np.pi/(2*a), 0.0], [np.pi/(2*a), np.pi/(2*a)]])
    result = compute_lao_sto_band_and_realspace(k_path, tl, th, td, dE, dSO, dRSO, a, Nx, Ny, s)
    assert isinstance(result, tuple) and len(result) == 3, \
        f'Must return 3-tuple, got {type(result)} len={len(result) if isinstance(result,tuple) else "N/A"}'
    bands, rs_eigs, max_dev = result
    assert isinstance(bands, np.ndarray)
    assert bands.shape == (3, 6), f'bands shape: expected (3,6), got {bands.shape}'
    assert isinstance(rs_eigs, np.ndarray)
    assert rs_eigs.shape == (6 * Nx * Ny,), \
        f'rs_eigs shape: expected ({6*Nx*Ny},), got {rs_eigs.shape}'
    assert isinstance(max_dev, float), f'max_dev must be float, got {type(max_dev)}'
    assert max_dev >= 0.0
Test case 2 of 5
import numpy as np
from scipy.sparse import lil_matrix, csr_matrix
from scipy.sparse.linalg import eigsh
from scipy.linalg import eigvalsh


def test_case_2():
    tl, th, td, dE, dSO, dRSO = 875.0, 40.0, 40.0, 47.0, 10.0, 5.0
    a = 0.394
    Nx, Ny, s = 2, 2, 1.0
    k_path = np.array([[0.0, 0.0], [np.pi/(2*a), np.pi/(2*a)]])
    _, _, max_dev = compute_lao_sto_band_and_realspace(
        k_path, tl, th, td, dE, dSO, dRSO, a, Nx, Ny, s)
    assert max_dev < 1e-8, \
        f'At s=1, max_deviation must be <1e-8 meV (Bloch theorem), got {max_dev:.4e}'
Test case 3 of 5
import numpy as np
from scipy.sparse import lil_matrix, csr_matrix
from scipy.sparse.linalg import eigsh
from scipy.linalg import eigvalsh


def test_case_3():
    tl, th, td, dE = 875.0, 40.0, 40.0, 47.0
    a = 0.394
    k_path = np.array([[0.0, 0.0]])
    Nx, Ny, s = 1, 1, 1.0
    result = compute_lao_sto_band_and_realspace(k_path, tl, th, td, dE, 0.0, 0.0, a, Nx, Ny, s)
    bands = result[0]
    eigs = np.sort(bands[0])
    assert np.allclose(eigs[:2], -dE, atol=1e-8), \
        f'dxy bands at Gamma should be {-dE} meV, got {eigs[:2]}'
    assert np.allclose(eigs[2:], 0.0, atol=1e-8), \
        f'dxz/dyz bands at Gamma should be 0 meV, got {eigs[2:]}'
Test case 4 of 5
import numpy as np
from scipy.sparse import lil_matrix, csr_matrix
from scipy.sparse.linalg import eigsh
from scipy.linalg import eigvalsh


def test_case_4():
    tl, th, td, dE, dSO, dRSO = 875.0, 40.0, 40.0, 47.0, 10.0, 5.0
    a = 0.394
    k_path = np.linspace(0, np.pi/a, 5)
    k_path_2d = np.column_stack([k_path, np.zeros(5)])
    Nx, Ny, s = 2, 2, 1.0
    bands, rs_eigs, max_dev = compute_lao_sto_band_and_realspace(
        k_path_2d, tl, th, td, dE, dSO, dRSO, a, Nx, Ny, s)
    assert bands.shape == (5, 6), f'bands shape: expected (5,6), got {bands.shape}'
    for i in range(5):
        assert np.all(bands[i, :-1] <= bands[i, 1:]), \
            f'Bands at k-point {i} not sorted: {bands[i,:]}'
    assert np.all(np.isfinite(bands)), 'All bands must be finite'
    assert bands[0, 0] < 0, f'Lowest band at Gamma must be negative, got {bands[0,0]}'
    assert np.all(rs_eigs[:-1] <= rs_eigs[1:]), 'rs_eigs must be sorted ascending'
    assert max_dev < 1e-8, f'At s=1 max_dev={max_dev:.4e} should be <1e-8'
Test case 5 of 5
import numpy as np
from scipy.sparse import lil_matrix, csr_matrix
from scipy.sparse.linalg import eigsh
from scipy.linalg import eigvalsh


def test_case_5():
    tl, th, td, dE, dSO, dRSO = 875.0, 40.0, 40.0, 47.0, 10.0, 5.0
    a = 0.394
    Nx, Ny, s = 3, 3, 1.0
    k_path = np.array([[0.0, 0.0], [np.pi/(2*a), 0.0]])
    bands, rs_eigs, max_dev = compute_lao_sto_band_and_realspace(
        k_path, tl, th, td, dE, dSO, dRSO, a, Nx, Ny, s)
    assert bands.shape == (2, 6), f'Expected bands shape (2,6), got {bands.shape}'
    assert rs_eigs.shape == (6*Nx*Ny,), \
        f'Expected rs_eigs shape ({6*Nx*Ny},), got {rs_eigs.shape}'
    assert np.all(np.isfinite(rs_eigs)), 'Real-space eigenvalues must be finite'
    assert max_dev < 1e-8, \
        f'At s=1, 3x3 lattice: max_dev={max_dev:.4e} must be <1e-8 (Bloch theorem)'
    assert rs_eigs[0] < 0, f'Lowest real-space eigenvalue must be negative, got {rs_eigs[0]}'
Solution
import numpy as np
from scipy.linalg import eigvalsh
from scipy.sparse import lil_matrix

def compute_lao_sto_band_and_realspace(k_path, tl, th, td, dE, dSO, dRSO, a, Nx, Ny, s):
    '''Compute the LAO/STO 2DEG band structure and real-space spectrum.

    Parameters:
      k_path : ndarray, shape (N_k, 2)
          Array of (kx, ky) values in units of 1/nm along the desired path.
      tl : float
          Light-band hopping parameter in meV.
      th : float
          Heavy-band hopping parameter in meV.
      td : float
          Hybridisation (next-nearest-neighbour) hopping in meV.
      dE : float
          Crystal-field splitting in meV (dxy lowered by this amount).
      dSO : float
          Atomic spin-orbit coupling strength in meV.
      dRSO : float
          Rashba spin-orbit coupling strength in meV.
      a : float
          Lattice constant in nm.
      Nx : int
          Number of sites in x-direction for real-space lattice.
      Ny : int
          Number of sites in y-direction for real-space lattice.
      s : float
          Scaling factor for the real-space mesh.

    Returns:
      bands : ndarray, shape (N_k, 6)
          Sorted eigenvalues at each k-point in meV.
      rs_eigs : ndarray, shape (6*Nx*Ny,)
          Sorted eigenvalues from real-space Hamiltonian in meV.
      max_dev : float
          Maximum absolute deviation between sorted k-space and real-space spectra.
    '''
    import numpy as np
    from scipy.linalg import eigvalsh
    from scipy.sparse import lil_matrix

    def compute_kspace_bandstructure(k_path, tl, th, td, dE, dSO, dRSO, a):
        k_path = np.asarray(k_path)
        N_k = k_path.shape[0]
        bands = np.zeros((N_k, 6))
        sx = np.array([[0, 1], [1, 0]], dtype=complex)
        sy = np.array([[0, -1j], [1j, 0]], dtype=complex)
        sz = np.array([[1, 0], [0, -1]], dtype=complex)
        I2 = np.eye(2, dtype=complex)
        H_SO = np.zeros((6, 6), dtype=complex)
        H_SO[0:2, 2:4] = -(dSO / 3.0) * (1j * sx)
        H_SO[0:2, 4:6] = -(dSO / 3.0) * (-1j * sy)
        H_SO[2:4, 0:2] = -(dSO / 3.0) * (-1j * sx)
        H_SO[2:4, 4:6] = -(dSO / 3.0) * (1j * sz)
        H_SO[4:6, 0:2] = -(dSO / 3.0) * (1j * sy)
        H_SO[4:6, 2:4] = -(dSO / 3.0) * (-1j * sz)
        for i, (kx, ky) in enumerate(k_path):
            cxa = np.cos(kx * a)
            cya = np.cos(ky * a)
            sxa = np.sin(kx * a)
            sya = np.sin(ky * a)
            eps_xy = 4*tl - 2*tl*cxa - 2*tl*cya - dE
            eps_xz = 2*tl + 2*th - 2*tl*cxa - 2*th*cya
            eps_yz = 2*tl + 2*th - 2*th*cxa - 2*tl*cya
            eps_h  = 2*td * sxa * sya
            H0_orb = np.array([[eps_xy, 0, 0], [0, eps_xz, eps_h], [0, eps_h, eps_yz]], dtype=complex)
            H0 = np.kron(H0_orb, I2)
            H_RSO_orb = dRSO * np.array([[0, 1j*sya, 1j*sxa], [-1j*sya, 0, 0], [-1j*sxa, 0, 0]], dtype=complex)
            H_RSO = np.kron(H_RSO_orb, I2)
            H = H0 + H_SO + H_RSO
            bands[i] = eigvalsh(H)
        return bands

    def build_scaled_hopping_blocks(tl, th, td, dRSO, s):
        tl_p   = tl   / s**2
        th_p   = th   / s**2
        td_p   = td   / s**2
        dRSO_p = dRSO / s
        I2 = np.eye(2, dtype=complex)
        Hx_orb_kin = np.diag([-tl_p, -tl_p, -th_p]).astype(complex)
        Hx_orb_rso = (dRSO_p / 2.0) * np.array([[0, 0, -1], [0, 0, 0], [1, 0, 0]], dtype=complex)
        Hx = np.kron(Hx_orb_kin + Hx_orb_rso, I2)
        Hy_orb_kin = np.diag([-tl_p, -th_p, -tl_p]).astype(complex)
        Hy_orb_rso = (dRSO_p / 2.0) * np.array([[0, -1, 0], [1, 0, 0], [0, 0, 0]], dtype=complex)
        Hy = np.kron(Hy_orb_kin + Hy_orb_rso, I2)
        Hmix_orb = (td_p / 2.0) * np.array([[0, 0, 0], [0, 0, 1], [0, 1, 0]], dtype=complex)
        Hmix = np.kron(Hmix_orb, I2)
        return {'Hx': Hx, 'Hy': Hy, 'Hmix': Hmix}

    def build_and_diagonalise_realspace(Nx, Ny, tl, th, td, dE, dSO, dRSO, s):
        tl_p = tl / s**2
        th_p = th / s**2
        blocks = build_scaled_hopping_blocks(tl, th, td, dRSO, s)
        Hx = blocks['Hx']
        Hy = blocks['Hy']
        Hmix = blocks['Hmix']
        sx = np.array([[0, 1], [1, 0]], dtype=complex)
        sy = np.array([[0, -1j], [1j, 0]], dtype=complex)
        sz = np.array([[1, 0], [0, -1]], dtype=complex)
        I2 = np.eye(2, dtype=complex)
        H_SO = np.zeros((6, 6), dtype=complex)
        H_SO[0:2, 2:4] = -(dSO / 3.0) * (1j * sx)
        H_SO[0:2, 4:6] = -(dSO / 3.0) * (-1j * sy)
        H_SO[2:4, 0:2] = -(dSO / 3.0) * (-1j * sx)
        H_SO[2:4, 4:6] = -(dSO / 3.0) * (1j * sz)
        H_SO[4:6, 0:2] = -(dSO / 3.0) * (1j * sy)
        H_SO[4:6, 2:4] = -(dSO / 3.0) * (-1j * sz)
        H_onsite_kin_orb = np.diag([4*tl_p - dE, 2*tl_p + 2*th_p, 2*tl_p + 2*th_p]).astype(complex)
        H_on = np.kron(H_onsite_kin_orb, I2) + H_SO
        D = 6 * Nx * Ny
        H = lil_matrix((D, D), dtype=complex)
        for mu in range(Nx):
            for nu in range(Ny):
                n = mu * Ny + nu
                H[6*n : 6*(n+1), 6*n : 6*(n+1)] = H_on
                n_xp1 = ((mu + 1) % Nx) * Ny + nu
                H[6*n : 6*(n+1), 6*n_xp1 : 6*(n_xp1+1)] += Hx
                H[6*n_xp1 : 6*(n_xp1+1), 6*n : 6*(n+1)] += Hx.conj().T
                n_yp1 = mu * Ny + ((nu + 1) % Ny)
                H[6*n : 6*(n+1), 6*n_yp1 : 6*(n_yp1+1)] += Hy
                H[6*n_yp1 : 6*(n_yp1+1), 6*n : 6*(n+1)] += Hy.conj().T
                n_xp1_ym1 = ((mu + 1) % Nx) * Ny + ((nu - 1 + Ny) % Ny)
                H[6*n : 6*(n+1), 6*n_xp1_ym1 : 6*(n_xp1_ym1+1)] += Hmix
                H[6*n_xp1_ym1 : 6*(n_xp1_ym1+1), 6*n : 6*(n+1)] += Hmix.conj().T
                n_xp1_yp1 = ((mu + 1) % Nx) * Ny + ((nu + 1) % Ny)
                H[6*n : 6*(n+1), 6*n_xp1_yp1 : 6*(n_xp1_yp1+1)] -= Hmix
                H[6*n_xp1_yp1 : 6*(n_xp1_yp1+1), 6*n : 6*(n+1)] -= Hmix.conj().T
        eigs = eigvalsh(H.toarray())
        return eigs

    def verify_kspace_realspace_agreement(Nx, Ny, tl, th, td, dE, dSO, dRSO, a, s):
        if s <= 0:
            raise ValueError(f'Scaling factor s must be positive, got s={s}')
        mu_vals = np.arange(Nx)
        nu_vals = np.arange(Ny)
        k_grid = np.array(np.meshgrid(mu_vals, nu_vals)).T.reshape(-1, 2)
        kx_vals = 2 * np.pi * k_grid[:, 0] / (Nx * s * a)
        ky_vals = 2 * np.pi * k_grid[:, 1] / (Ny * s * a)
        commensurate_k_path = np.column_stack([kx_vals, ky_vals])
        k_bands = compute_kspace_bandstructure(commensurate_k_path, tl, th, td, dE, dSO, dRSO, a)
        kspace_eigs = np.sort(k_bands.flatten())
        realspace_eigs = build_and_diagonalise_realspace(Nx, Ny, tl, th, td, dE, dSO, dRSO, s)
        max_deviation = float(np.max(np.abs(kspace_eigs - realspace_eigs)))
        return kspace_eigs, realspace_eigs, max_deviation

    bands = compute_kspace_bandstructure(k_path, tl, th, td, dE, dSO, dRSO, a)
    
    _, rs_eigs, max_dev = verify_kspace_realspace_agreement(
        Nx, Ny, tl, th, td, dE, dSO, dRSO, a, s)

    return bands, rs_eigs, max_dev

Prompt

Compute the band structure of the LAO/STO t2gt_{2g} 2DEG along a k-path by assembling and diagonalising the full 6×66 \times 6 Hamiltonian H^(k)=H^0(k)+H^SO+H^RSO(k)\hat{H}(\mathbf{k}) = \hat{H}_0(\mathbf{k}) + \hat{H}_{SO} + \hat{H}_{RSO}(\mathbf{k}) at each of NkN_k k-points. The basis ordering is (dxy ⁣↑,dxy ⁣↓,dxz ⁣↑,dxz ⁣↓,dyz ⁣↑,dyz ⁣↓)(d_{xy}\!\uparrow, d_{xy}\!\downarrow, d_{xz}\!\uparrow, d_{xz}\!\downarrow, d_{yz}\!\uparrow, d_{yz}\!\downarrow). Kinetic term H^0\hat{H}_0: The 3×33 \times 3 orbital block is H0orb=(ϵxy000ϵxzϵh0ϵhϵyz)H_0^{\text{orb}} = \begin{pmatrix} \epsilon_{xy} & 0 & 0 \\ 0 & \epsilon_{xz} & \epsilon_h \\ 0 & \epsilon_h & \epsilon_{yz} \end{pmatrix} where ϵxy=4tl−2tlcos⁡(kxa)−2tlcos⁡(kya)−ΔE\epsilon_{xy} = 4t_l - 2t_l\cos(k_x a) - 2t_l\cos(k_y a) - \Delta E, ϵxz=2tl+2th−2tlcos⁡(kxa)−2thcos⁡(kya)\epsilon_{xz} = 2t_l + 2t_h - 2t_l\cos(k_x a) - 2t_h\cos(k_y a), ϵyz=2tl+2th−2thcos⁡(kxa)−2tlcos⁡(kya)\epsilon_{yz} = 2t_l + 2t_h - 2t_h\cos(k_x a) - 2t_l\cos(k_y a), and ϵh=2tdsin⁡(kxa)sin⁡(kya)\epsilon_h = 2t_d \sin(k_x a)\sin(k_y a). The full kinetic block is H^0=H0orb⊗σ^0\hat{H}_0 = H_0^{\text{orb}} \otimes \hat{\sigma}_0 (Kronecker product with 2×22 \times 2 identity). Atomic SOC H^SO\hat{H}_{SO}: This is the t2gt_{2g}-projected L⋅S\mathbf{L}\cdot\mathbf{S} interaction. The effective angular momentum operators in the t2gt_{2g} subspace (basis: dxy,dxz,dyzd_{xy}, d_{xz}, d_{yz}) are related to the ℓ=1\ell=1 operators by a sign inversion: L^αt2g=−L^αℓ=1\hat{L}^{t_{2g}}_\alpha = -\hat{L}^{\ell=1}_\alpha. The resulting 6×66 \times 6 SOC matrix in the (orbital⊗spin)(\text{orbital} \otimes \text{spin}) basis is H^SO=−ΔSO3(0iσ^x−iσ^y−iσ^x0iσ^ziσ^y−iσ^z0)\hat{H}_{SO} = -\frac{\Delta_{SO}}{3}\begin{pmatrix} 0 & i\hat{\sigma}_x & -i\hat{\sigma}_y \\ -i\hat{\sigma}_x & 0 & i\hat{\sigma}_z \\ i\hat{\sigma}_y & -i\hat{\sigma}_z & 0 \end{pmatrix} where each entry is a 2×22 \times 2 Pauli matrix block and the rows/columns index orbitals in the order (dxyd_{xy}, dxzd_{xz}, dyzd_{yz}). The Pauli matrices are σ^x=(0110)\hat{\sigma}_x = \begin{pmatrix}0&1\\1&0\end{pmatrix}, σ^y=(0−ii0)\hat{\sigma}_y = \begin{pmatrix}0&-i\\i&0\end{pmatrix}, σ^z=(100−1)\hat{\sigma}_z = \begin{pmatrix}1&0\\0&-1\end{pmatrix}. This matrix is k-independent and Hermitian. Its eigenvalues form a quartet at −ΔSO/3-\Delta_{SO}/3 and a doublet at +2ΔSO/3+2\Delta_{SO}/3. Rashba SOC H^RSO\hat{H}_{RSO}: The Rashba coupling at the LAO/STO interface is purely orbital off-diagonal and does not mix spin. The 3×33 \times 3 orbital block is HRSOorb=ΔRSO(0isin⁡(kya)isin⁡(kxa)−isin⁡(kya)00−isin⁡(kxa)00)H_{RSO}^{\text{orb}} = \Delta_{RSO}\begin{pmatrix} 0 & i\sin(k_y a) & i\sin(k_x a) \\ -i\sin(k_y a) & 0 & 0 \\ -i\sin(k_x a) & 0 & 0 \end{pmatrix} and the full 6×66 \times 6 Rashba matrix is H^RSO=HRSOorb⊗σ^0\hat{H}_{RSO} = H_{RSO}^{\text{orb}} \otimes \hat{\sigma}_0. The sin⁡(kya)\sin(k_y a) entry couples dxy↔dxzd_{xy} \leftrightarrow d_{xz}, and sin⁡(kxa)\sin(k_x a) couples dxy↔dyzd_{xy} \leftrightarrow d_{yz}. There is no dxz↔dyzd_{xz} \leftrightarrow d_{yz} Rashba coupling. At each k-point, sum all three terms to form the 6×66 \times 6 Hermitian Hamiltonian and compute its 6 sorted eigenvalues. Return the full band structure array. Implement the following function:

def compute_kspace_bandstructure(k_path, tl, th, td, dE, dSO, dRSO, a):
    '''Compute band structure along a k-path for the LAO/STO 2DEG.

    Parameters:
      k_path : ndarray, shape (N_k, 2)
          Array of (kx, ky) values in units of 1/nm.
      tl : float
          Light-band hopping in meV.
      th : float
          Heavy-band hopping in meV.
      td : float
          Hybridisation hopping in meV.
      dE : float
          Crystal-field splitting in meV.
      dSO : float
          Atomic spin-orbit coupling in meV.
      dRSO : float
          Rashba spin-orbit coupling in meV.
      a : float
          Lattice constant in nm.

    Returns:
      bands : ndarray, shape (N_k, 6)
          Sorted eigenvalues at each k-point in meV.
    '''
    return bands

Background

The t2gt_{2g} tight-binding model for the LAO/STO 2DEG describes three dd-orbitals (dxyd_{xy}, dxzd_{xz}, dyzd_{yz}) with anisotropic hopping on a square lattice. The kinetic energy is diagonal in orbital space except for the dxzd_{xz}/dyzd_{yz} hybridisation term ϵh=2tdsin⁡(kxa)sin⁡(kya)\epsilon_h = 2t_d \sin(k_x a)\sin(k_y a). The atomic spin-orbit coupling arises from projecting the full L⋅S\mathbf{L}\cdot\mathbf{S} operator onto the t2gt_{2g} subspace, where the effective angular momentum operators pick up a sign inversion relative to the standard ℓ=1\ell=1 matrices, yielding an overall minus sign in the H^SO\hat{H}_{SO} prefactor. The Rashba spin-orbit coupling at this oxide interface has an unusual orbital structure — it couples only dxyd_{xy} to the other two orbitals and is a scalar in spin space, unlike the standard semiconductor Rashba term α(kyσx−kxσy)\alpha(k_y\sigma_x - k_x\sigma_y) which acts purely in spin space.

Function signature

def compute_kspace_bandstructure(k_path, tl, th, td, dE, dSO, dRSO, a):
    return bands

Published by the benchmark's authors. Opening these shows how the problem is solved and what the graders assert.

Test case 1 of 6
import numpy as np
from scipy.sparse import lil_matrix, csr_matrix
from scipy.sparse.linalg import eigsh
from scipy.linalg import eigvalsh


def test_case_sub1_1():
    dSO = 30.0
    a = 0.394
    k_path = np.array([[0.0, 0.0]])
    bands = compute_kspace_bandstructure(k_path, tl=0.0, th=0.0, td=0.0, dE=0.0,
                                         dSO=dSO, dRSO=0.0, a=a)
    expected_low = -dSO / 3.0   # = -10 meV
    expected_high = 2.0 * dSO / 3.0  # = +20 meV
    eigs = np.sort(bands[0])
    # BREAKS if: model uses standard ell=1 angular momentum matrices without t2g sign inversion,
    # producing quartet at +dSO/3 instead of the correct -dSO/3
    assert np.allclose(eigs[:4], expected_low, atol=1e-8), \
        f'SOC quartet should be at {expected_low} meV (4-fold degenerate), got {eigs[:4]}'
    assert np.allclose(eigs[4:], expected_high, atol=1e-8), \
        f'SOC doublet should be at {expected_high} meV (2-fold degenerate), got {eigs[4:]}'
    assert bands.shape == (1, 6), f'Expected shape (1,6), got {bands.shape}'
Test case 2 of 6
import numpy as np
from scipy.sparse import lil_matrix, csr_matrix
from scipy.sparse.linalg import eigsh
from scipy.linalg import eigvalsh


def test_case_sub1_2():
    tl, th, td, dE = 875.0, 40.0, 40.0, 47.0
    a = 0.394
    k_path = np.array([[0.0, 0.0]])
    bands = compute_kspace_bandstructure(k_path, tl=tl, th=th, td=td, dE=dE,
                                         dSO=0.0, dRSO=0.0, a=a)
    eigs = np.sort(bands[0])
    # BREAKS if: kinetic dispersion formula is wrong
    assert np.allclose(eigs[:2], -dE, atol=1e-8), \
        f'Gamma-point dxy bands should be {-dE} meV (spin-degenerate), got {eigs[:2]}'
    assert np.allclose(eigs[2:], 0.0, atol=1e-8), \
        f'Gamma-point dxz/dyz bands should be 0 meV, got {eigs[2:]}'
    assert bands.shape == (1, 6), f'Expected shape (1,6), got {bands.shape}'
Test case 3 of 6
import numpy as np
from scipy.sparse import lil_matrix, csr_matrix
from scipy.sparse.linalg import eigsh
from scipy.linalg import eigvalsh


def test_case_sub1_3():
    tl, th, td, dE, dSO = 100.0, 20.0, 10.0, 30.0, 20.0
    a = 0.394
    k_path = np.array([[0.0, 0.0]])
    bands_norso = compute_kspace_bandstructure(k_path, tl=tl, th=th, td=td, dE=dE,
                                               dSO=dSO, dRSO=0.0, a=a)
    bands_withrso = compute_kspace_bandstructure(k_path, tl=tl, th=th, td=td, dE=dE,
                                                  dSO=dSO, dRSO=500.0, a=a)
    # BREAKS if: Rashba has a k-independent component (wrong implementation)
    assert np.allclose(bands_norso[0], bands_withrso[0], atol=1e-8), \
        f'Rashba must vanish exactly at k=(0,0): no-RSO={bands_norso[0]}, with-RSO={bands_withrso[0]}'
    assert bands_norso.shape == (1, 6)
    assert bands_withrso.shape == (1, 6)
Test case 4 of 6
import numpy as np
from scipy.sparse import lil_matrix, csr_matrix
from scipy.sparse.linalg import eigsh
from scipy.linalg import eigvalsh


def test_case_sub1_4():
    dSO = 10.0
    # Build H_SO with t2g sign inversion (overall minus vs standard ell=1):
    sx = np.array([[0, 1], [1, 0]], dtype=complex)
    sy = np.array([[0, -1j], [1j, 0]], dtype=complex)
    sz = np.array([[1, 0], [0, -1]], dtype=complex)
    H_SO_ref = np.zeros((6, 6), dtype=complex)
    H_SO_ref[0:2, 2:4] = -(1j * sx * dSO / 3.0)
    H_SO_ref[0:2, 4:6] = -(-1j * sy * dSO / 3.0)
    H_SO_ref[2:4, 0:2] = -(-1j * sx * dSO / 3.0)
    H_SO_ref[2:4, 4:6] = -(1j * sz * dSO / 3.0)
    H_SO_ref[4:6, 0:2] = -(1j * sy * dSO / 3.0)
    H_SO_ref[4:6, 2:4] = -(-1j * sz * dSO / 3.0)
    expected_eigs = np.sort(eigvalsh(H_SO_ref).real)
    k_path = np.array([[0.0, 0.0]])
    bands = compute_kspace_bandstructure(k_path, tl=0.0, th=0.0, td=0.0, dE=0.0,
                                         dSO=dSO, dRSO=0.0, a=0.394)
    got = np.sort(bands[0])
    # BREAKS if: model uses standard L.S without t2g sign inversion
    assert np.allclose(got, expected_eigs, atol=1e-8), \
        f'SOC eigenvalues mismatch. Expected {expected_eigs}, got {got}'
    assert np.allclose(got[:4], -dSO/3.0, atol=1e-8), \
        f'Quartet should be at {-dSO/3.0}, got {got[:4]}'
    assert np.allclose(got[4:], 2*dSO/3.0, atol=1e-8), \
        f'Doublet should be at {2*dSO/3.0}, got {got[4:]}'
Test case 5 of 6
import numpy as np
from scipy.sparse import lil_matrix, csr_matrix
from scipy.sparse.linalg import eigsh
from scipy.linalg import eigvalsh


def test_case_sub1_5():
    tl, th, td, dE = 875.0, 40.0, 40.0, 47.0
    a = 0.394
    kx = np.pi / (3 * a)
    k_path = np.array([[kx, 0.0]])
    bands = compute_kspace_bandstructure(k_path, tl=tl, th=th, td=td, dE=dE,
                                         dSO=0.0, dRSO=0.0, a=a)
    eps_xy = 2*tl*(1 - np.cos(kx*a)) - dE
    eps_xz = 2*tl*(1 - np.cos(kx*a))
    eps_yz = 2*th*(1 - np.cos(kx*a))
    eigs = np.sort(bands[0])
    expected = np.sort([eps_xy, eps_xy, eps_xz, eps_xz, eps_yz, eps_yz])
    # BREAKS if: hybridisation is nonzero at ky=0
    assert np.allclose(eigs, expected, atol=1e-8), \
        f'At ky=0, expected {expected} meV, got {eigs}'
    assert bands.shape == (1, 6)
    assert np.all(np.isfinite(eigs))
Test case 6 of 6
import numpy as np
from scipy.sparse import lil_matrix, csr_matrix
from scipy.sparse.linalg import eigsh
from scipy.linalg import eigvalsh


def test_case_sub1_6():
    tl, th, td, dE, dSO, dRSO = 875.0, 40.0, 40.0, 47.0, 10.0, 5.0
    a = 0.394
    kx = np.pi / (2 * a)
    ky = np.pi / (3 * a)
    k_path = np.array([[0.0, 0.0], [kx, 0.0], [kx, ky]])
    bands = compute_kspace_bandstructure(k_path, tl=tl, th=th, td=td, dE=dE,
                                         dSO=dSO, dRSO=dRSO, a=a)
    assert isinstance(bands, np.ndarray), f'Must return ndarray, got {type(bands)}'
    assert bands.shape == (3, 6), f'Expected shape (3,6), got {bands.shape}'
    assert np.all(np.isfinite(bands)), 'All eigenvalues must be finite'
    for i in range(3):
        assert np.all(bands[i, :-1] <= bands[i, 1:]), \
            f'Bands at k-point {i} not sorted: {bands[i,:]}'
    # BREAKS if: crystal field not applied correctly
    assert bands[0, 0] < 0, f'Lowest band at Gamma must be negative, got {bands[0,0]}'
Solution
import numpy as np
from scipy.linalg import eigvalsh

def compute_kspace_bandstructure(k_path, tl, th, td, dE, dSO, dRSO, a):
    '''Compute band structure along a k-path for the LAO/STO 2DEG.

    Parameters:
      k_path : ndarray, shape (N_k, 2)
          Array of (kx, ky) values in units of 1/nm.
      tl : float
          Light-band hopping in meV.
      th : float
          Heavy-band hopping in meV.
      td : float
          Hybridisation hopping in meV.
      dE : float
          Crystal-field splitting in meV.
      dSO : float
          Atomic spin-orbit coupling in meV.
      dRSO : float
          Rashba spin-orbit coupling in meV.
      a : float
          Lattice constant in nm.

    Returns:
      bands : ndarray, shape (N_k, 6)
          Sorted eigenvalues at each k-point in meV.
    '''
    k_path = np.asarray(k_path)
    N_k = k_path.shape[0]
    bands = np.zeros((N_k, 6))

    sx = np.array([[0, 1], [1, 0]], dtype=complex)
    sy = np.array([[0, -1j], [1j, 0]], dtype=complex)
    sz = np.array([[1, 0], [0, -1]], dtype=complex)
    I2 = np.eye(2, dtype=complex)

    # H_SO = -(dSO/3)*block_matrix(Pauli sub-blocks), Eq.(5) of paper
    # Overall minus from t2g projection: L^t2g = -L^(ell=1)
    # Produces quartet at -dSO/3 and doublet at +2*dSO/3
    H_SO = np.zeros((6, 6), dtype=complex)
    H_SO[0:2, 2:4] = -(dSO / 3.0) * (1j * sx)
    H_SO[0:2, 4:6] = -(dSO / 3.0) * (-1j * sy)
    H_SO[2:4, 0:2] = -(dSO / 3.0) * (-1j * sx)
    H_SO[2:4, 4:6] = -(dSO / 3.0) * (1j * sz)
    H_SO[4:6, 0:2] = -(dSO / 3.0) * (1j * sy)
    H_SO[4:6, 2:4] = -(dSO / 3.0) * (-1j * sz)

    for i, (kx, ky) in enumerate(k_path):
        cxa = np.cos(kx * a)
        cya = np.cos(ky * a)
        sxa = np.sin(kx * a)
        sya = np.sin(ky * a)

        eps_xy = 4*tl - 2*tl*cxa - 2*tl*cya - dE  # Eq.(1)
        eps_xz = 2*tl + 2*th - 2*tl*cxa - 2*th*cya  # Eq.(2)
        eps_yz = 2*tl + 2*th - 2*th*cxa - 2*tl*cya  # Eq.(3)
        eps_h  = 2*td * sxa * sya  # Eq.(4)

        H0_orb = np.array([[eps_xy, 0,      0     ],
                            [0,      eps_xz, eps_h ],
                            [0,      eps_h,  eps_yz]], dtype=complex)
        H0 = np.kron(H0_orb, I2)

        # Rashba: orbital off-diagonal, scalar in spin space (⊗σ₀), Eq.(6)
        H_RSO_orb = dRSO * np.array([[0,       1j*sya, 1j*sxa],
                                      [-1j*sya, 0,      0     ],
                                      [-1j*sxa, 0,      0     ]], dtype=complex)
        H_RSO = np.kron(H_RSO_orb, I2)

        H = H0 + H_SO + H_RSO
        bands[i] = eigvalsh(H)

    return bands

Prompt

Construct the 6×66 \times 6 hopping blocks (H^x\hat{H}_x, H^y\hat{H}_y, H^mix\hat{H}_{mix}) for the scaled real-space tight-binding Hamiltonian of the LAO/STO 2DEG. The scaling factor ss coarsens the mesh from lattice constant aa to dx=s⋅adx = s \cdot a. Kinetic hopping parameters are rescaled as tl′=tl/s2t_l' = t_l / s^2, th′=th/s2t_h' = t_h / s^2, td′=td/s2t_d' = t_d / s^2. The Rashba coupling rescales as ΔRSO′=ΔRSO/s\Delta_{RSO}' = \Delta_{RSO} / s. The factor of 1/21/2 in the Rashba hopping entries arises from the discrete representation of sin⁡(ka)\sin(k a) as (eika−e−ika)/(2i)(e^{ika} - e^{-ika})/(2i). x-direction hopping block H^x\hat{H}_x: This connects site (μ,ν)(\mu, \nu) to site (μ+1,ν)(\mu+1, \nu). The orbital-diagonal kinetic part reflects the anisotropy: along xx, the dxyd_{xy} and dxzd_{xz} orbitals have light hopping while dyzd_{yz} has heavy hopping. The Rashba part couples dxy↔dyzd_{xy} \leftrightarrow d_{yz} (the orbital pair linked by sin⁡(kxa)\sin(k_x a) in k-space): H^x=diag(−tl′,−tl′,−th′)⊗σ^0+ΔRSO′2(00−1000100)⊗σ^0\hat{H}_x = \text{diag}(-t_l', -t_l', -t_h') \otimes \hat{\sigma}_0 + \frac{\Delta_{RSO}'}{2}\begin{pmatrix}0 & 0 & -1 \\ 0 & 0 & 0 \\ 1 & 0 & 0\end{pmatrix} \otimes \hat{\sigma}_0 y-direction hopping block H^y\hat{H}_y: This connects site (μ,ν)(\mu, \nu) to site (μ,ν+1)(\mu, \nu+1). Along yy, dxyd_{xy} and dyzd_{yz} have light hopping while dxzd_{xz} has heavy hopping. The Rashba part couples dxy↔dxzd_{xy} \leftrightarrow d_{xz} (linked by sin⁡(kya)\sin(k_y a)): H^y=diag(−tl′,−th′,−tl′)⊗σ^0+ΔRSO′2(0−10100000)⊗σ^0\hat{H}_y = \text{diag}(-t_l', -t_h', -t_l') \otimes \hat{\sigma}_0 + \frac{\Delta_{RSO}'}{2}\begin{pmatrix}0 & -1 & 0 \\ 1 & 0 & 0 \\ 0 & 0 & 0\end{pmatrix} \otimes \hat{\sigma}_0 Diagonal hopping block H^mix\hat{H}_{mix}: This connects next-nearest neighbours along the lattice diagonal. It couples ONLY dxzd_{xz} and dyzd_{yz} orbitals (not dxyd_{xy}): H^mix=td′2(000001010)⊗σ^0\hat{H}_{mix} = \frac{t_d'}{2}\begin{pmatrix}0 & 0 & 0 \\ 0 & 0 & 1 \\ 0 & 1 & 0\end{pmatrix} \otimes \hat{\sigma}_0 The function returns a dictionary with keys ’Hx’\texttt{'Hx'}, ’Hy’\texttt{'Hy'}, ’Hmix’\texttt{'Hmix'}, each mapping to a 6×66 \times 6 complex numpy array. The physical origin of the orbital-direction anisotropy is that dxzd_{xz} extends along xx (light hopping along xx, heavy along yy) while dyzd_{yz} extends along yy (heavy along xx, light along yy). The dxyd_{xy} orbital has equal (light) hopping in both directions. Implement the following function:

def build_scaled_hopping_blocks(tl, th, td, dRSO, s):
    '''Build the scaled real-space hopping blocks.

    Parameters:
      tl : float
          Light-band hopping in meV (unscaled).
      th : float
          Heavy-band hopping in meV (unscaled).
      td : float
          Hybridisation hopping in meV (unscaled).
      dRSO : float
          Rashba coupling in meV (unscaled).
      s : float
          Scaling factor (s=1 means no scaling).

    Returns:
      blocks : dict
          Dictionary with keys 'Hx', 'Hy', 'Hmix', each mapping
          to a complex ndarray of shape (6, 6).
    '''
    return {'Hx': Hx, 'Hy': Hy, 'Hmix': Hmix}

Background

In the scaled real-space representation, hopping between neighbouring sites on the coarsened lattice (spacing dx=s⋅adx = s \cdot a) is captured by 6×66 \times 6 blocks that encode both the kinetic energy and the spatial part of the Rashba coupling. The discrete Fourier representation of sin⁡(kxa)\sin(k_x a) generates a factor of 1/21/2 in the Rashba hopping entries. The three blocks (H^x\hat{H}_x, H^y\hat{H}_y, H^mix\hat{H}_{mix}) have distinct orbital structures reflecting the directional anisotropy of the t2gt_{2g} orbitals and the particular orbital pairs coupled by Rashba terms along each direction.

Function signature

def build_scaled_hopping_blocks(tl, th, td, dRSO, s):
    return {'Hx': Hx, 'Hy': Hy, 'Hmix': Hmix}

Published by the benchmark's authors. Opening these shows how the problem is solved and what the graders assert.

Test case 1 of 6
import numpy as np
from scipy.sparse import lil_matrix, csr_matrix
from scipy.sparse.linalg import eigsh
from scipy.linalg import eigvalsh


def test_case_sub2_1():
    tl, th, td, dRSO, s = 875.0, 40.0, 40.0, 5.0, 1.0
    blocks = build_scaled_hopping_blocks(tl, th, td, dRSO, s)
    Hy = blocks['Hy']
    tl_p = tl / s**2
    th_p = th / s**2
    # BREAKS if: model assigns heavy hopping to dyz in Hy
    assert np.isclose(Hy[0, 0].real, -tl_p, atol=1e-10), \
        f'Hy[0,0] (dxy-up kinetic) should be {-tl_p}, got {Hy[0,0]}'
    assert np.isclose(Hy[2, 2].real, -th_p, atol=1e-10), \
        f'Hy[2,2] (dxz-up kinetic) should be {-th_p} (heavy along y), got {Hy[2,2]}'
    assert np.isclose(Hy[4, 4].real, -tl_p, atol=1e-10), \
        f'Hy[4,4] (dyz-up kinetic) should be {-tl_p} (light along y), got {Hy[4,4]}'
    assert np.isclose(Hy[1, 1].real, -tl_p, atol=1e-10)
    assert np.isclose(Hy[3, 3].real, -th_p, atol=1e-10)
    assert np.isclose(Hy[5, 5].real, -tl_p, atol=1e-10)
Test case 2 of 6
import numpy as np
from scipy.sparse import lil_matrix, csr_matrix
from scipy.sparse.linalg import eigsh
from scipy.linalg import eigvalsh


def test_case_sub2_2():
    tl, th, td, dRSO, s = 0.0, 0.0, 0.0, 10.0, 1.0
    blocks = build_scaled_hopping_blocks(tl, th, td, dRSO, s)
    Hx = blocks['Hx']
    dRSO_p = dRSO / s
    rashba_amp = dRSO_p / 2.0
    # BREAKS if: model uses dxy<->dxz coupling (wrong orbital pair for Hx Rashba)
    assert np.isclose(Hx[0, 4].real, -rashba_amp, atol=1e-10), \
        f'Hx[0,4] (dxy->dyz Rashba) should be {-rashba_amp}, got {Hx[0,4]}'
    assert np.isclose(Hx[4, 0].real, +rashba_amp, atol=1e-10), \
        f'Hx[4,0] (dyz->dxy Rashba) should be {+rashba_amp}, got {Hx[4,0]}'
    assert np.isclose(Hx[0, 2].real, 0.0, atol=1e-10), \
        f'Hx[0,2] (dxy->dxz) must be zero in Hx, got {Hx[0,2]}'
    assert np.isclose(Hx[1, 5].real, -rashba_amp, atol=1e-10)
    assert np.isclose(Hx[5, 1].real, +rashba_amp, atol=1e-10)
Test case 3 of 6
import numpy as np
from scipy.sparse import lil_matrix, csr_matrix
from scipy.sparse.linalg import eigsh
from scipy.linalg import eigvalsh


def test_case_sub2_3():
    tl, th, td, dRSO, s = 100.0, 20.0, 30.0, 40.0, 2.0
    blocks = build_scaled_hopping_blocks(tl, th, td, dRSO, s)
    tl_p = tl / s**2
    th_p = th / s**2
    td_p = td / s**2
    dRSO_p = dRSO / s
    Hx = blocks['Hx']
    Hy = blocks['Hy']
    Hmix = blocks['Hmix']
    # BREAKS if: scaling exponents are wrong
    assert np.isclose(Hx[0, 0].real, -tl_p, atol=1e-10), f'Hx[0,0]=-tl_p expected {-tl_p}, got {Hx[0,0]}'
    assert np.isclose(Hx[2, 2].real, -tl_p, atol=1e-10), f'Hx[2,2]=-tl_p expected {-tl_p}, got {Hx[2,2]}'
    assert np.isclose(Hx[4, 4].real, -th_p, atol=1e-10), f'Hx[4,4]=-th_p expected {-th_p}, got {Hx[4,4]}'
    assert np.isclose(Hy[0, 0].real, -tl_p, atol=1e-10)
    assert np.isclose(Hy[2, 2].real, -th_p, atol=1e-10)
    assert np.isclose(Hy[4, 4].real, -tl_p, atol=1e-10)
    assert np.isclose(abs(Hx[0, 4]), dRSO_p / 2.0, atol=1e-10), \
        f'|Hx[0,4]|=dRSO_p/2={dRSO_p/2}, got {abs(Hx[0,4])}'
    assert np.isclose(Hmix[2, 4].real, td_p / 2.0, atol=1e-10), \
        f'Hmix[2,4]=td_p/2={td_p/2}, got {Hmix[2,4]}'
    assert np.isclose(Hmix[4, 2].real, td_p / 2.0, atol=1e-10)
Test case 4 of 6
import numpy as np
from scipy.sparse import lil_matrix, csr_matrix
from scipy.sparse.linalg import eigsh
from scipy.linalg import eigvalsh


def test_case_sub2_4():
    tl, th, td, dRSO, s = 875.0, 40.0, 40.0, 5.0, 1.0
    blocks = build_scaled_hopping_blocks(tl, th, td, dRSO, s)
    Hmix = blocks['Hmix']
    # BREAKS if: Hmix has nonzero dxy entries
    assert np.allclose(Hmix[0, :], 0.0, atol=1e-15), \
        f'Hmix row 0 (dxy spin-up) must be zero, got {Hmix[0,:]}'
    assert np.allclose(Hmix[1, :], 0.0, atol=1e-15), \
        f'Hmix row 1 (dxy spin-dn) must be zero, got {Hmix[1,:]}'
    assert np.allclose(Hmix[:, 0], 0.0, atol=1e-15), \
        f'Hmix col 0 (dxy spin-up) must be zero, got {Hmix[:,0]}'
    assert np.allclose(Hmix[:, 1], 0.0, atol=1e-15), \
        f'Hmix col 1 (dxy spin-dn) must be zero, got {Hmix[:,1]}'
    assert not np.allclose(Hmix[2, 4], 0.0, atol=1e-10), \
        f'Hmix[2,4] (dxz->dyz) must be nonzero, got {Hmix[2,4]}'
Test case 5 of 6
import numpy as np
from scipy.sparse import lil_matrix, csr_matrix
from scipy.sparse.linalg import eigsh
from scipy.linalg import eigvalsh


def test_case_sub2_5():
    tl, th, td, dRSO, s = 0.0, 0.0, 0.0, 20.0, 1.0
    blocks = build_scaled_hopping_blocks(tl, th, td, dRSO, s)
    Hy = blocks['Hy']
    dRSO_p = dRSO / s
    rashba_amp = dRSO_p / 2.0
    # BREAKS if: Rashba in Hy couples wrong orbital pair
    assert np.isclose(Hy[0, 2].real, -rashba_amp, atol=1e-10), \
        f'Hy[0,2] (dxy->dxz Rashba) should be {-rashba_amp}, got {Hy[0,2]}'
    assert np.isclose(Hy[2, 0].real, +rashba_amp, atol=1e-10), \
        f'Hy[2,0] (dxz->dxy Rashba) should be {+rashba_amp}, got {Hy[2,0]}'
    assert np.isclose(Hy[0, 4].real, 0.0, atol=1e-10), \
        f'Hy[0,4] (dxy->dyz) must be zero in Hy, got {Hy[0,4]}'
    assert np.isclose(Hy[4, 0].real, 0.0, atol=1e-10), \
        f'Hy[4,0] (dyz->dxy) must be zero in Hy, got {Hy[4,0]}'
Test case 6 of 6
import numpy as np
from scipy.sparse import lil_matrix, csr_matrix
from scipy.sparse.linalg import eigsh
from scipy.linalg import eigvalsh


def test_case_sub2_6():
    tl, th, td, dRSO, s = 875.0, 40.0, 40.0, 5.0, 1.0
    blocks = build_scaled_hopping_blocks(tl, th, td, dRSO, s)
    assert isinstance(blocks, dict), f'Must return dict, got {type(blocks)}'
    for key in ['Hx', 'Hy', 'Hmix']:
        assert key in blocks, f'Missing key {key} in blocks dict'
        B = blocks[key]
        assert isinstance(B, np.ndarray), f'blocks["{key}"] must be ndarray, got {type(B)}'
        assert B.shape == (6, 6), f'blocks["{key}"] shape must be (6,6), got {B.shape}'
    Hmix = blocks['Hmix']
    td_p = td / s**2
    assert np.isclose(Hmix[2, 4].real, td_p / 2.0, atol=1e-10), \
        f'Hmix[2,4] should be td_p/2={td_p/2}, got {Hmix[2,4]}'
Solution
import numpy as np

def build_scaled_hopping_blocks(tl, th, td, dRSO, s):
    '''Build the scaled real-space hopping blocks.

    Parameters:
      tl : float
          Light-band hopping in meV (unscaled).
      th : float
          Heavy-band hopping in meV (unscaled).
      td : float
          Hybridisation hopping in meV (unscaled).
      dRSO : float
          Rashba coupling in meV (unscaled).
      s : float
          Scaling factor (s=1 means no scaling).

    Returns:
      blocks : dict
          Dictionary with keys 'Hx', 'Hy', 'Hmix', each mapping
          to a complex ndarray of shape (6, 6).
    '''
    tl_p   = tl   / s**2  # kinetic scales as 1/s^2, Eq.(9)
    th_p   = th   / s**2
    td_p   = td   / s**2
    dRSO_p = dRSO / s     # Rashba scales as 1/s, Eq.(9)

    I2 = np.eye(2, dtype=complex)

    # Hx: along x, dxy and dxz have light hopping, dyz has heavy hopping; Eq.(10)
    Hx_orb_kin = np.diag([-tl_p, -tl_p, -th_p]).astype(complex)
    Hx_orb_rso = (dRSO_p / 2.0) * np.array([[0, 0, -1],
                                              [0, 0,  0],
                                              [1, 0,  0]], dtype=complex)
    Hx = np.kron(Hx_orb_kin + Hx_orb_rso, I2)

    # Hy: along y, dxy and dyz have light hopping, dxz has heavy hopping; Eq.(11)
    Hy_orb_kin = np.diag([-tl_p, -th_p, -tl_p]).astype(complex)
    Hy_orb_rso = (dRSO_p / 2.0) * np.array([[ 0, -1, 0],
                                              [ 1,  0, 0],
                                              [ 0,  0, 0]], dtype=complex)
    Hy = np.kron(Hy_orb_kin + Hy_orb_rso, I2)

    # Hmix: diagonal hopping couples only dxz↔dyz; factor 1/2 from sin discrete FT; Eq.(12)
    Hmix_orb = (td_p / 2.0) * np.array([[0, 0, 0],
                                          [0, 0, 1],
                                          [0, 1, 0]], dtype=complex)
    Hmix = np.kron(Hmix_orb, I2)

    return {'Hx': Hx, 'Hy': Hy, 'Hmix': Hmix}

Prompt

Assemble the full sparse real-space tight-binding Hamiltonian for the LAO/STO 2DEG on an Nx×NyN_x \times N_y periodic lattice with 66 internal degrees of freedom per site, using the hopping blocks from build_scaled_hopping_blocks and the on-site block constructed internally. The total Hamiltonian dimension is D=6NxNyD = 6 N_x N_y. Each site (μ,ν)(\mu, \nu) maps to a global site index n=μ⋅Ny+νn = \mu \cdot N_y + \nu (with μ=0,…,Nx−1\mu = 0, \ldots, N_x - 1 and ν=0,…,Ny−1\nu = 0, \ldots, N_y - 1), and the 6×66 \times 6 block for site nn occupies rows/columns [6n,6n+6)[6n, 6n+6). The on-site block at each site combines the kinetic on-site energy, the crystal-field splitting, and the atomic SOC: H^on=(4tl′−ΔE0002tl′+2th′0002tl′+2th′)⊗σ^0+H^SO\hat{H}_{\text{on}} = \begin{pmatrix} 4t_l' - \Delta E & 0 & 0 \\ 0 & 2t_l' + 2t_h' & 0 \\ 0 & 0 & 2t_l' + 2t_h' \end{pmatrix} \otimes \hat{\sigma}_0 + \hat{H}_{SO} where tl′=tl/s2t_l' = t_l/s^2, th′=th/s2t_h' = t_h/s^2, and H^SO\hat{H}_{SO} is the unscaled atomic SOC matrix (identical to the k-space version defined in compute_kspace_bandstructure, with coupling strength ΔSO\Delta_{SO} and the same overall minus sign in its prefactor: H^SO=−ΔSO3(0iσ^x−iσ^y−iσ^x0iσ^ziσ^y−iσ^z0)\hat{H}_{SO} = -\frac{\Delta_{SO}}{3}\begin{pmatrix} 0 & i\hat{\sigma}_x & -i\hat{\sigma}_y \\ -i\hat{\sigma}_x & 0 & i\hat{\sigma}_z \\ i\hat{\sigma}_y & -i\hat{\sigma}_z & 0 \end{pmatrix}). The atomic SOC is NOT rescaled by 1/s21/s^2 or 1/s1/s — it enters the real-space Hamiltonian with the identical numerical value as in the k-space model. The on-site potential VsiteV_{\text{site}} may be added as Vsite⋅I^6V_{\text{site}} \cdot \hat{I}_6 (here set to zero for a uniform system). The nearest-neighbour hopping: site (μ,ν)(\mu, \nu) connects to (μ+1,ν)(\mu+1, \nu) via H^x\hat{H}_x and to (μ,ν+1)(\mu, \nu+1) via H^y\hat{H}_y, with Hermitian conjugates H^x†\hat{H}_x^\dagger and H^y†\hat{H}_y^\dagger for the reverse directions. Periodic boundary conditions wrap indices modulo NxN_x and NyN_y. The diagonal hopping: site (μ,ν)(\mu, \nu) connects to (μ+1,ν−1)(\mu+1, \nu-1) via +H^mix+\hat{H}_{mix} and to (μ+1,ν+1)(\mu+1, \nu+1) via −H^mix-\hat{H}_{mix}. The two diagonal directions carry opposite signs — this is essential for reproducing the sin⁡(kxa)sin⁡(kya)\sin(k_x a)\sin(k_y a) hybridisation in k-space. The reverse diagonal hops are the Hermitian conjugates: (μ−1,ν+1)(\mu-1, \nu+1) gets +H^mix†+\hat{H}_{mix}^\dagger and (μ−1,ν−1)(\mu-1, \nu-1) gets −H^mix†-\hat{H}_{mix}^\dagger. Assemble the sparse Hamiltonian, convert to a dense or CSR format as appropriate, and compute ALL DD eigenvalues sorted in ascending order. This function must call build_scaled_hopping_blocks to obtain the hopping blocks. Implement the following function:

def build_and_diagonalise_realspace(Nx, Ny, tl, th, td, dE, dSO, dRSO, s):
    """
    Build and diagonalise the real-space Hamiltonian.

    Parameters:
      Nx : int
          Number of sites in x-direction.
      Ny : int
          Number of sites in y-direction.
      tl : float
          Light-band hopping in meV (unscaled).
      th : float
          Heavy-band hopping in meV (unscaled).
      td : float
          Hybridisation hopping in meV (unscaled).
      dE : float
          Crystal-field splitting in meV.
      dSO : float
          Atomic spin-orbit coupling in meV.
      dRSO : float
          Rashba coupling in meV (unscaled).
      s : float
          Scaling factor.

    Returns:
      eigs : ndarray, shape (6*Nx*Ny,)
          All eigenvalues sorted in ascending order, in meV.
    """
    return eigs

Background

The real-space tight-binding Hamiltonian is assembled by placing 6×66 \times 6 blocks at appropriate positions in a D×DD \times D sparse matrix, where D=6NxNyD = 6 N_x N_y. The on-site block includes both scaled kinetic on-site energies and the unscaled atomic SOC. The nearest-neighbour and diagonal hopping blocks connect adjacent sites with periodic boundary conditions. The critical structural feature is that the diagonal hopping block H^mix\hat{H}_{mix} enters with opposite signs for the two diagonal directions (μ+1,ν−1)(\mu+1, \nu-1) and (μ+1,ν+1)(\mu+1, \nu+1), which is the real-space origin of the sin⁡(kxa)sin⁡(kya)\sin(k_x a)\sin(k_y a) hybridisation term in k-space. This step uses build_scaled_hopping_blocks from the previous subproblem and internally constructs the atomic SOC matrix H^SO\hat{H}_{SO} using the same convention as compute_kspace_bandstructure.

Function signature

def build_and_diagonalise_realspace(Nx, Ny, tl, th, td, dE, dSO, dRSO, s):
    return eigs

Published by the benchmark's authors. Opening these shows how the problem is solved and what the graders assert.

Test case 1 of 5
import numpy as np
from scipy.sparse import lil_matrix, csr_matrix
from scipy.sparse.linalg import eigsh
from scipy.linalg import eigvalsh

def build_scaled_hopping_blocks(tl, th, td, dRSO, s):
    '''Build the scaled real-space hopping blocks.

    Parameters:
      tl : float
          Light-band hopping in meV (unscaled).
      th : float
          Heavy-band hopping in meV (unscaled).
      td : float
          Hybridisation hopping in meV (unscaled).
      dRSO : float
          Rashba coupling in meV (unscaled).
      s : float
          Scaling factor (s=1 means no scaling).

    Returns:
      blocks : dict
          Dictionary with keys 'Hx', 'Hy', 'Hmix', each mapping
          to a complex ndarray of shape (6, 6).
    '''
    tl_p   = tl   / s**2  # kinetic scales as 1/s^2, Eq.(9)
    th_p   = th   / s**2
    td_p   = td   / s**2
    dRSO_p = dRSO / s     # Rashba scales as 1/s, Eq.(9)

    I2 = np.eye(2, dtype=complex)

    # Hx: along x, dxy and dxz have light hopping, dyz has heavy hopping; Eq.(10)
    Hx_orb_kin = np.diag([-tl_p, -tl_p, -th_p]).astype(complex)
    Hx_orb_rso = (dRSO_p / 2.0) * np.array([[0, 0, -1],
                                              [0, 0,  0],
                                              [1, 0,  0]], dtype=complex)
    Hx = np.kron(Hx_orb_kin + Hx_orb_rso, I2)

    # Hy: along y, dxy and dyz have light hopping, dxz has heavy hopping; Eq.(11)
    Hy_orb_kin = np.diag([-tl_p, -th_p, -tl_p]).astype(complex)
    Hy_orb_rso = (dRSO_p / 2.0) * np.array([[ 0, -1, 0],
                                              [ 1,  0, 0],
                                              [ 0,  0, 0]], dtype=complex)
    Hy = np.kron(Hy_orb_kin + Hy_orb_rso, I2)

    # Hmix: diagonal hopping couples only dxz↔dyz; factor 1/2 from sin discrete FT; Eq.(12)
    Hmix_orb = (td_p / 2.0) * np.array([[0, 0, 0],
                                          [0, 0, 1],
                                          [0, 1, 0]], dtype=complex)
    Hmix = np.kron(Hmix_orb, I2)

    return {'Hx': Hx, 'Hy': Hy, 'Hmix': Hmix}


def test_case_sub3_1():
    tl, th, td, dE, dSO, dRSO = 100.0, 20.0, 40.0, 30.0, 0.0, 0.0
    a = 0.394
    Nx, Ny, s = 2, 2, 1.0
    eigs_rs = build_and_diagonalise_realspace(Nx, Ny, tl, th, td, dE, dSO, dRSO, s)
    all_kspace_eigs = []
    for mu in range(Nx):
        for nu in range(Ny):
            kx = 2 * np.pi * mu / (Nx * s * a)
            ky = 2 * np.pi * nu / (Ny * s * a)
            eps_xy = 4*tl - 2*tl*np.cos(kx*a) - 2*tl*np.cos(ky*a) - dE
            eps_xz = 2*tl + 2*th - 2*tl*np.cos(kx*a) - 2*th*np.cos(ky*a)
            eps_yz = 2*tl + 2*th - 2*th*np.cos(kx*a) - 2*tl*np.cos(ky*a)
            eps_h = 2*td*np.sin(kx*a)*np.sin(ky*a)
            H0_orb = np.array([[eps_xy, 0, 0], [0, eps_xz, eps_h], [0, eps_h, eps_yz]])
            H0 = np.kron(H0_orb, np.eye(2))
            ek = eigvalsh(H0)
            all_kspace_eigs.extend(ek)
    kspace_sorted = np.sort(all_kspace_eigs)
    rs_sorted = np.sort(eigs_rs)
    # BREAKS if: model uses same sign for both diagonal directions
    assert np.allclose(rs_sorted, kspace_sorted, atol=1e-6), \
        f'RS and k-space spectra must match at s=1. Max diff={np.max(np.abs(rs_sorted - kspace_sorted)):.4e}'
    assert len(rs_sorted) == 6*Nx*Ny
    assert len(kspace_sorted) == 6*Nx*Ny
Test case 2 of 5
import numpy as np
from scipy.sparse import lil_matrix, csr_matrix
from scipy.sparse.linalg import eigsh
from scipy.linalg import eigvalsh

def build_scaled_hopping_blocks(tl, th, td, dRSO, s):
    '''Build the scaled real-space hopping blocks.

    Parameters:
      tl : float
          Light-band hopping in meV (unscaled).
      th : float
          Heavy-band hopping in meV (unscaled).
      td : float
          Hybridisation hopping in meV (unscaled).
      dRSO : float
          Rashba coupling in meV (unscaled).
      s : float
          Scaling factor (s=1 means no scaling).

    Returns:
      blocks : dict
          Dictionary with keys 'Hx', 'Hy', 'Hmix', each mapping
          to a complex ndarray of shape (6, 6).
    '''
    tl_p   = tl   / s**2  # kinetic scales as 1/s^2, Eq.(9)
    th_p   = th   / s**2
    td_p   = td   / s**2
    dRSO_p = dRSO / s     # Rashba scales as 1/s, Eq.(9)

    I2 = np.eye(2, dtype=complex)

    # Hx: along x, dxy and dxz have light hopping, dyz has heavy hopping; Eq.(10)
    Hx_orb_kin = np.diag([-tl_p, -tl_p, -th_p]).astype(complex)
    Hx_orb_rso = (dRSO_p / 2.0) * np.array([[0, 0, -1],
                                              [0, 0,  0],
                                              [1, 0,  0]], dtype=complex)
    Hx = np.kron(Hx_orb_kin + Hx_orb_rso, I2)

    # Hy: along y, dxy and dyz have light hopping, dxz has heavy hopping; Eq.(11)
    Hy_orb_kin = np.diag([-tl_p, -th_p, -tl_p]).astype(complex)
    Hy_orb_rso = (dRSO_p / 2.0) * np.array([[ 0, -1, 0],
                                              [ 1,  0, 0],
                                              [ 0,  0, 0]], dtype=complex)
    Hy = np.kron(Hy_orb_kin + Hy_orb_rso, I2)

    # Hmix: diagonal hopping couples only dxz↔dyz; factor 1/2 from sin discrete FT; Eq.(12)
    Hmix_orb = (td_p / 2.0) * np.array([[0, 0, 0],
                                          [0, 0, 1],
                                          [0, 1, 0]], dtype=complex)
    Hmix = np.kron(Hmix_orb, I2)

    return {'Hx': Hx, 'Hy': Hy, 'Hmix': Hmix}


def test_case_sub3_2():
    tl, th, td, dE, dSO, dRSO = 100.0, 20.0, 30.0, 50.0, 0.0, 0.0
    s = 1.0
    tl_p = tl / s**2
    th_p = th / s**2
    eigs = build_and_diagonalise_realspace(1, 1, tl, th, td, dE, dSO, dRSO, s)
    # For 1x1 periodic lattice, site hops to itself via PBCs:
    # On-site: diag(4tl'-dE, 2tl'+2th', 2tl'+2th') = diag(350, 240, 240)
    # x-hop to self: Hx+Hx^dag = 2*diag(-tl',-tl',-th') = diag(-200,-200,-40)
    # y-hop to self: Hy+Hy^dag = 2*diag(-tl',-th',-tl') = diag(-200,-40,-200)
    # Diagonal hops cancel (+Hmix and -Hmix both to same site)
    # Total: diag(350-200-200, 240-200-40, 240-40-200) = diag(-50, 0, 0) x I2
    expected = np.sort([-50.0, -50.0, 0.0, 0.0, 0.0, 0.0])
    got = np.sort(eigs)
    # BREAKS if: PBC self-hopping not accounted for
    assert np.allclose(got, expected, atol=1e-8), \
        f'Single-site eigenvalues: expected {expected}, got {got}'
    assert eigs.shape == (6,), f'Expected shape (6,), got {eigs.shape}'
    assert np.all(np.isfinite(eigs))
Test case 3 of 5
import numpy as np
from scipy.sparse import lil_matrix, csr_matrix
from scipy.sparse.linalg import eigsh
from scipy.linalg import eigvalsh

def build_scaled_hopping_blocks(tl, th, td, dRSO, s):
    '''Build the scaled real-space hopping blocks.

    Parameters:
      tl : float
          Light-band hopping in meV (unscaled).
      th : float
          Heavy-band hopping in meV (unscaled).
      td : float
          Hybridisation hopping in meV (unscaled).
      dRSO : float
          Rashba coupling in meV (unscaled).
      s : float
          Scaling factor (s=1 means no scaling).

    Returns:
      blocks : dict
          Dictionary with keys 'Hx', 'Hy', 'Hmix', each mapping
          to a complex ndarray of shape (6, 6).
    '''
    tl_p   = tl   / s**2  # kinetic scales as 1/s^2, Eq.(9)
    th_p   = th   / s**2
    td_p   = td   / s**2
    dRSO_p = dRSO / s     # Rashba scales as 1/s, Eq.(9)

    I2 = np.eye(2, dtype=complex)

    # Hx: along x, dxy and dxz have light hopping, dyz has heavy hopping; Eq.(10)
    Hx_orb_kin = np.diag([-tl_p, -tl_p, -th_p]).astype(complex)
    Hx_orb_rso = (dRSO_p / 2.0) * np.array([[0, 0, -1],
                                              [0, 0,  0],
                                              [1, 0,  0]], dtype=complex)
    Hx = np.kron(Hx_orb_kin + Hx_orb_rso, I2)

    # Hy: along y, dxy and dyz have light hopping, dxz has heavy hopping; Eq.(11)
    Hy_orb_kin = np.diag([-tl_p, -th_p, -tl_p]).astype(complex)
    Hy_orb_rso = (dRSO_p / 2.0) * np.array([[ 0, -1, 0],
                                              [ 1,  0, 0],
                                              [ 0,  0, 0]], dtype=complex)
    Hy = np.kron(Hy_orb_kin + Hy_orb_rso, I2)

    # Hmix: diagonal hopping couples only dxz↔dyz; factor 1/2 from sin discrete FT; Eq.(12)
    Hmix_orb = (td_p / 2.0) * np.array([[0, 0, 0],
                                          [0, 0, 1],
                                          [0, 1, 0]], dtype=complex)
    Hmix = np.kron(Hmix_orb, I2)

    return {'Hx': Hx, 'Hy': Hy, 'Hmix': Hmix}


def test_case_sub3_3():
    dSO = 30.0
    s = 2.0
    eigs = build_and_diagonalise_realspace(1, 1, tl=0.0, th=0.0, td=0.0, dE=0.0,
                                            dSO=dSO, dRSO=0.0, s=s)
    got = np.sort(eigs)
    expected_low = -dSO / 3.0    # = -10 meV (unscaled)
    expected_high = 2.0 * dSO / 3.0  # = +20 meV (unscaled)
    # BREAKS if: atomic SOC incorrectly scaled by 1/s^2, giving -2.5 instead of -10
    assert np.allclose(got[:4], expected_low, atol=1e-8), \
        f'SOC quartet must be unscaled at {expected_low} meV, got {got[:4]}. If scaled by 1/s^2 would be {expected_low/s**2}'
    assert np.allclose(got[4:], expected_high, atol=1e-8), \
        f'SOC doublet must be unscaled at {expected_high} meV, got {got[4:]}'
    assert eigs.shape == (6,)
Test case 4 of 5
import numpy as np
from scipy.sparse import lil_matrix, csr_matrix
from scipy.sparse.linalg import eigsh
from scipy.linalg import eigvalsh

def build_scaled_hopping_blocks(tl, th, td, dRSO, s):
    '''Build the scaled real-space hopping blocks.

    Parameters:
      tl : float
          Light-band hopping in meV (unscaled).
      th : float
          Heavy-band hopping in meV (unscaled).
      td : float
          Hybridisation hopping in meV (unscaled).
      dRSO : float
          Rashba coupling in meV (unscaled).
      s : float
          Scaling factor (s=1 means no scaling).

    Returns:
      blocks : dict
          Dictionary with keys 'Hx', 'Hy', 'Hmix', each mapping
          to a complex ndarray of shape (6, 6).
    '''
    tl_p   = tl   / s**2  # kinetic scales as 1/s^2, Eq.(9)
    th_p   = th   / s**2
    td_p   = td   / s**2
    dRSO_p = dRSO / s     # Rashba scales as 1/s, Eq.(9)

    I2 = np.eye(2, dtype=complex)

    # Hx: along x, dxy and dxz have light hopping, dyz has heavy hopping; Eq.(10)
    Hx_orb_kin = np.diag([-tl_p, -tl_p, -th_p]).astype(complex)
    Hx_orb_rso = (dRSO_p / 2.0) * np.array([[0, 0, -1],
                                              [0, 0,  0],
                                              [1, 0,  0]], dtype=complex)
    Hx = np.kron(Hx_orb_kin + Hx_orb_rso, I2)

    # Hy: along y, dxy and dyz have light hopping, dxz has heavy hopping; Eq.(11)
    Hy_orb_kin = np.diag([-tl_p, -th_p, -tl_p]).astype(complex)
    Hy_orb_rso = (dRSO_p / 2.0) * np.array([[ 0, -1, 0],
                                              [ 1,  0, 0],
                                              [ 0,  0, 0]], dtype=complex)
    Hy = np.kron(Hy_orb_kin + Hy_orb_rso, I2)

    # Hmix: diagonal hopping couples only dxz↔dyz; factor 1/2 from sin discrete FT; Eq.(12)
    Hmix_orb = (td_p / 2.0) * np.array([[0, 0, 0],
                                          [0, 0, 1],
                                          [0, 1, 0]], dtype=complex)
    Hmix = np.kron(Hmix_orb, I2)

    return {'Hx': Hx, 'Hy': Hy, 'Hmix': Hmix}


def test_case_sub3_4():
    tl, th, td, dE, dSO, dRSO = 0.0, 0.0, 40.0, 0.0, 0.0, 0.0
    a = 0.394
    Nx, Ny, s = 3, 3, 1.0
    eigs_rs = build_and_diagonalise_realspace(Nx, Ny, tl, th, td, dE, dSO, dRSO, s)
    all_kspace_eigs = []
    for mu in range(Nx):
        for nu in range(Ny):
            kx = 2 * np.pi * mu / (Nx * s * a)
            ky = 2 * np.pi * nu / (Ny * s * a)
            eps_h = 2*td*np.sin(kx*a)*np.sin(ky*a)
            H0_orb = np.array([[0.0, 0, 0], [0, 0, eps_h], [0, eps_h, 0]])
            ek = eigvalsh(np.kron(H0_orb, np.eye(2)))
            all_kspace_eigs.extend(ek)
    kspace_sorted = np.sort(all_kspace_eigs)
    rs_sorted = np.sort(eigs_rs)
    # BREAKS if: diagonal hopping uses same sign for both diagonals
    assert np.allclose(rs_sorted, kspace_sorted, atol=1e-6), \
        f'3x3 RS vs kspace max diff={np.max(np.abs(rs_sorted-kspace_sorted)):.4e}'
    assert eigs_rs.shape == (6*Nx*Ny,), f'Expected shape ({6*Nx*Ny},), got {eigs_rs.shape}'
    assert len(kspace_sorted) == 6*Nx*Ny
Test case 5 of 5
import numpy as np
from scipy.sparse import lil_matrix, csr_matrix
from scipy.sparse.linalg import eigsh
from scipy.linalg import eigvalsh

def build_scaled_hopping_blocks(tl, th, td, dRSO, s):
    '''Build the scaled real-space hopping blocks.

    Parameters:
      tl : float
          Light-band hopping in meV (unscaled).
      th : float
          Heavy-band hopping in meV (unscaled).
      td : float
          Hybridisation hopping in meV (unscaled).
      dRSO : float
          Rashba coupling in meV (unscaled).
      s : float
          Scaling factor (s=1 means no scaling).

    Returns:
      blocks : dict
          Dictionary with keys 'Hx', 'Hy', 'Hmix', each mapping
          to a complex ndarray of shape (6, 6).
    '''
    tl_p   = tl   / s**2  # kinetic scales as 1/s^2, Eq.(9)
    th_p   = th   / s**2
    td_p   = td   / s**2
    dRSO_p = dRSO / s     # Rashba scales as 1/s, Eq.(9)

    I2 = np.eye(2, dtype=complex)

    # Hx: along x, dxy and dxz have light hopping, dyz has heavy hopping; Eq.(10)
    Hx_orb_kin = np.diag([-tl_p, -tl_p, -th_p]).astype(complex)
    Hx_orb_rso = (dRSO_p / 2.0) * np.array([[0, 0, -1],
                                              [0, 0,  0],
                                              [1, 0,  0]], dtype=complex)
    Hx = np.kron(Hx_orb_kin + Hx_orb_rso, I2)

    # Hy: along y, dxy and dyz have light hopping, dxz has heavy hopping; Eq.(11)
    Hy_orb_kin = np.diag([-tl_p, -th_p, -tl_p]).astype(complex)
    Hy_orb_rso = (dRSO_p / 2.0) * np.array([[ 0, -1, 0],
                                              [ 1,  0, 0],
                                              [ 0,  0, 0]], dtype=complex)
    Hy = np.kron(Hy_orb_kin + Hy_orb_rso, I2)

    # Hmix: diagonal hopping couples only dxz↔dyz; factor 1/2 from sin discrete FT; Eq.(12)
    Hmix_orb = (td_p / 2.0) * np.array([[0, 0, 0],
                                          [0, 0, 1],
                                          [0, 1, 0]], dtype=complex)
    Hmix = np.kron(Hmix_orb, I2)

    return {'Hx': Hx, 'Hy': Hy, 'Hmix': Hmix}


def test_case_sub3_5():
    Nx, Ny = 2, 3
    eigs = build_and_diagonalise_realspace(Nx, Ny, tl=875.0, th=40.0, td=40.0,
                                            dE=47.0, dSO=10.0, dRSO=5.0, s=1.0)
    assert isinstance(eigs, np.ndarray), f'Must return ndarray, got {type(eigs)}'
    assert eigs.shape == (6 * Nx * Ny,), f'Expected shape ({6*Nx*Ny},), got {eigs.shape}'
    assert np.all(np.isfinite(eigs)), 'All eigenvalues must be finite'
    assert np.all(eigs[:-1] <= eigs[1:]), 'Eigenvalues must be sorted ascending'
    # BREAKS if: crystal field not applied
    assert eigs[0] < 0, f'Lowest eigenvalue should be negative (crystal field), got {eigs[0]}'
Solution
import numpy as np
from scipy.sparse import lil_matrix
from scipy.linalg import eigvalsh

def build_scaled_hopping_blocks(tl, th, td, dRSO, s):
    '''Build the scaled real-space hopping blocks.

    Parameters:
      tl : float
          Light-band hopping in meV (unscaled).
      th : float
          Heavy-band hopping in meV (unscaled).
      td : float
          Hybridisation hopping in meV (unscaled).
      dRSO : float
          Rashba coupling in meV (unscaled).
      s : float
          Scaling factor (s=1 means no scaling).

    Returns:
      blocks : dict
          Dictionary with keys 'Hx', 'Hy', 'Hmix', each mapping
          to a complex ndarray of shape (6, 6).
    '''
    tl_p   = tl   / s**2  # kinetic scales as 1/s^2, Eq.(9)
    th_p   = th   / s**2
    td_p   = td   / s**2
    dRSO_p = dRSO / s     # Rashba scales as 1/s, Eq.(9)

    I2 = np.eye(2, dtype=complex)

    # Hx: along x, dxy and dxz have light hopping, dyz has heavy hopping; Eq.(10)
    Hx_orb_kin = np.diag([-tl_p, -tl_p, -th_p]).astype(complex)
    Hx_orb_rso = (dRSO_p / 2.0) * np.array([[0, 0, -1],
                                              [0, 0,  0],
                                              [1, 0,  0]], dtype=complex)
    Hx = np.kron(Hx_orb_kin + Hx_orb_rso, I2)

    # Hy: along y, dxy and dyz have light hopping, dxz has heavy hopping; Eq.(11)
    Hy_orb_kin = np.diag([-tl_p, -th_p, -tl_p]).astype(complex)
    Hy_orb_rso = (dRSO_p / 2.0) * np.array([[ 0, -1, 0],
                                              [ 1,  0, 0],
                                              [ 0,  0, 0]], dtype=complex)
    Hy = np.kron(Hy_orb_kin + Hy_orb_rso, I2)

    # Hmix: diagonal hopping couples only dxz↔dyz; factor 1/2 from sin discrete FT; Eq.(12)
    Hmix_orb = (td_p / 2.0) * np.array([[0, 0, 0],
                                          [0, 0, 1],
                                          [0, 1, 0]], dtype=complex)
    Hmix = np.kron(Hmix_orb, I2)

    return {'Hx': Hx, 'Hy': Hy, 'Hmix': Hmix}

def build_and_diagonalise_realspace(Nx, Ny, tl, th, td, dE, dSO, dRSO, s):
    '''Build and diagonalise the real-space Hamiltonian.

    Parameters:
      Nx : int
          Number of sites in x-direction.
      Ny : int
          Number of sites in y-direction.
      tl : float
          Light-band hopping in meV (unscaled).
      th : float
          Heavy-band hopping in meV (unscaled).
      td : float
          Hybridisation hopping in meV (unscaled).
      dE : float
          Crystal-field splitting in meV.
      dSO : float
          Atomic spin-orbit coupling in meV.
      dRSO : float
          Rashba coupling in meV (unscaled).
      s : float
          Scaling factor.

    Returns:
      eigs : ndarray, shape (6*Nx*Ny,)
          All eigenvalues sorted in ascending order, in meV.
    '''
    tl_p = tl / s**2
    th_p = th / s**2

    blocks = build_scaled_hopping_blocks(tl, th, td, dRSO, s)
    Hx = blocks['Hx']
    Hy = blocks['Hy']
    Hmix = blocks['Hmix']

    sx = np.array([[0, 1], [1, 0]], dtype=complex)
    sy = np.array([[0, -1j], [1j, 0]], dtype=complex)
    sz = np.array([[1, 0], [0, -1]], dtype=complex)
    I2 = np.eye(2, dtype=complex)

    H_SO = np.zeros((6, 6), dtype=complex)
    H_SO[0:2, 2:4] = -(dSO / 3.0) * (1j * sx)
    H_SO[0:2, 4:6] = -(dSO / 3.0) * (-1j * sy)
    H_SO[2:4, 0:2] = -(dSO / 3.0) * (-1j * sx)
    H_SO[2:4, 4:6] = -(dSO / 3.0) * (1j * sz)
    H_SO[4:6, 0:2] = -(dSO / 3.0) * (1j * sy)
    H_SO[4:6, 2:4] = -(dSO / 3.0) * (-1j * sz)

    H_onsite_kin_orb = np.diag([4*tl_p - dE, 2*tl_p + 2*th_p, 2*tl_p + 2*th_p]).astype(complex)
    H_on = np.kron(H_onsite_kin_orb, I2) + H_SO

    D = 6 * Nx * Ny
    H = lil_matrix((D, D), dtype=complex)

    for mu in range(Nx):
        for nu in range(Ny):
            n = mu * Ny + nu
            H[6*n : 6*(n+1), 6*n : 6*(n+1)] = H_on

            # Hopping terms (forward directions)
            n_xp1 = ((mu + 1) % Nx) * Ny + nu
            H[6*n : 6*(n+1), 6*n_xp1 : 6*(n_xp1+1)] += Hx
            H[6*n_xp1 : 6*(n_xp1+1), 6*n : 6*(n+1)] += Hx.conj().T

            n_yp1 = mu * Ny + ((nu + 1) % Ny)
            H[6*n : 6*(n+1), 6*n_yp1 : 6*(n_yp1+1)] += Hy
            H[6*n_yp1 : 6*(n_yp1+1), 6*n : 6*(n+1)] += Hy.conj().T

            n_xp1_ym1 = ((mu + 1) % Nx) * Ny + ((nu - 1 + Ny) % Ny)
            H[6*n : 6*(n+1), 6*n_xp1_ym1 : 6*(n_xp1_ym1+1)] += Hmix
            H[6*n_xp1_ym1 : 6*(n_xp1_ym1+1), 6*n : 6*(n+1)] += Hmix.conj().T

            n_xp1_yp1 = ((mu + 1) % Nx) * Ny + ((nu + 1) % Ny)
            H[6*n : 6*(n+1), 6*n_xp1_yp1 : 6*(n_xp1_yp1+1)] -= Hmix
            H[6*n_xp1_yp1 : 6*(n_xp1_yp1+1), 6*n : 6*(n+1)] -= Hmix.conj().T

    eigs = eigvalsh(H.toarray())
    return eigs

Prompt

Verify the consistency of the k-space and real-space tight-binding Hamiltonians by computing the maximum absolute deviation between their sorted eigenvalue spectra on a finite periodic Nx×NyN_x \times N_y lattice. For a periodic lattice with Nx×NyN_x \times N_y sites and scaling factor ss, the commensurate k-points are kx(μ)=2πμNx⋅s⋅a,ky(ν)=2πνNy⋅s⋅ak_x^{(\mu)} = \frac{2\pi \mu}{N_x \cdot s \cdot a}, \quad k_y^{(\nu)} = \frac{2\pi \nu}{N_y \cdot s \cdot a} for μ=0,1,…,Nx−1\mu = 0, 1, \ldots, N_x - 1 and ν=0,1,…,Ny−1\nu = 0, 1, \ldots, N_y - 1. The k-space reference spectrum is obtained by calling compute_kspace_bandstructure at these Nx×NyN_x \times N_y commensurate k-points, collecting all 6NxNy6 N_x N_y eigenvalues, and sorting them. The real-space spectrum is obtained by calling build_and_diagonalise_realspace. Both spectra have length 6NxNy6 N_x N_y and, at scaling factor s=1s=1, must agree to machine precision (the maximum deviation should be <10−8< 10^{-8} meV). At scaling factors s>1s > 1, exact agreement is not expected over the full spectrum, but the low-energy eigenvalues (near the band bottom) should agree to within an error that scales with the higher-order terms dropped by the coarse-graining approximation. The function returns both sorted spectra and the maximum absolute deviation. This function must call both compute_kspace_bandstructure (subproblem 1) and build_and_diagonalise_realspace (subproblem 3). Implement the following function:

def verify_kspace_realspace_agreement(Nx, Ny, tl, th, td, dE, dSO, dRSO, a, s):
    '''Compare k-space and real-space spectra.

    Parameters:
      Nx : int
          Number of sites in x-direction.
      Ny : int
          Number of sites in y-direction.
      tl : float
          Light-band hopping in meV.
      th : float
          Heavy-band hopping in meV.
      td : float
          Hybridisation hopping in meV.
      dE : float
          Crystal-field splitting in meV.
      dSO : float
          Atomic spin-orbit coupling in meV.
      dRSO : float
          Rashba coupling in meV.
      a : float
          Lattice constant in nm.
      s : float
          Scaling factor.

    Returns:
      kspace_eigs : ndarray, shape (6*Nx*Ny,)
          Sorted eigenvalues from k-space at commensurate k-points, in meV.
      realspace_eigs : ndarray, shape (6*Nx*Ny,)
          Sorted eigenvalues from real-space Hamiltonian, in meV.
      max_deviation : float
          Maximum absolute difference in meV.
    '''
    return kspace_eigs, realspace_eigs, max_deviation

Background

The Bloch theorem guarantees that the eigenvalues of a periodic real-space tight-binding Hamiltonian on a finite Nx×NyN_x \times N_y lattice are exactly the union of the 66 eigenvalues at each of the Nx×NyN_x \times N_y commensurate k-points in the first Brillouin zone. The commensurate k-points are determined by the lattice size and the mesh spacing dx=s⋅adx = s \cdot a. Agreement at s=1s=1 serves as a stringent check on the correctness of both the k-space Hamiltonian assembly (including the atomic SOC sign pattern and Rashba orbital structure) and the real-space hopping blocks (including the opposite-sign diagonal hopping). This subproblem integrates compute_kspace_bandstructure and build_and_diagonalise_realspace from earlier subproblems.

Function signature

def verify_kspace_realspace_agreement(Nx, Ny, tl, th, td, dE, dSO, dRSO, a, s):
    return kspace_eigs, realspace_eigs, max_deviation

Published by the benchmark's authors. Opening these shows how the problem is solved and what the graders assert.

Test case 1 of 5
import numpy as np
from scipy.sparse import lil_matrix, csr_matrix
from scipy.sparse.linalg import eigsh
from scipy.linalg import eigvalsh
import pytest

def compute_kspace_bandstructure(k_path, tl, th, td, dE, dSO, dRSO, a):
    '''Compute band structure along a k-path for the LAO/STO 2DEG.

    Parameters:
      k_path : ndarray, shape (N_k, 2)
          Array of (kx, ky) values in units of 1/nm.
      tl : float
          Light-band hopping in meV.
      th : float
          Heavy-band hopping in meV.
      td : float
          Hybridisation hopping in meV.
      dE : float
          Crystal-field splitting in meV.
      dSO : float
          Atomic spin-orbit coupling in meV.
      dRSO : float
          Rashba spin-orbit coupling in meV.
      a : float
          Lattice constant in nm.

    Returns:
      bands : ndarray, shape (N_k, 6)
          Sorted eigenvalues at each k-point in meV.
    '''
    k_path = np.asarray(k_path)
    N_k = k_path.shape[0]
    bands = np.zeros((N_k, 6))

    sx = np.array([[0, 1], [1, 0]], dtype=complex)
    sy = np.array([[0, -1j], [1j, 0]], dtype=complex)
    sz = np.array([[1, 0], [0, -1]], dtype=complex)
    I2 = np.eye(2, dtype=complex)

    # H_SO = -(dSO/3)*block_matrix(Pauli sub-blocks), Eq.(5) of paper
    # Overall minus from t2g projection: L^t2g = -L^(ell=1)
    # Produces quartet at -dSO/3 and doublet at +2*dSO/3
    H_SO = np.zeros((6, 6), dtype=complex)
    H_SO[0:2, 2:4] = -(dSO / 3.0) * (1j * sx)
    H_SO[0:2, 4:6] = -(dSO / 3.0) * (-1j * sy)
    H_SO[2:4, 0:2] = -(dSO / 3.0) * (-1j * sx)
    H_SO[2:4, 4:6] = -(dSO / 3.0) * (1j * sz)
    H_SO[4:6, 0:2] = -(dSO / 3.0) * (1j * sy)
    H_SO[4:6, 2:4] = -(dSO / 3.0) * (-1j * sz)

    for i, (kx, ky) in enumerate(k_path):
        cxa = np.cos(kx * a)
        cya = np.cos(ky * a)
        sxa = np.sin(kx * a)
        sya = np.sin(ky * a)

        eps_xy = 4*tl - 2*tl*cxa - 2*tl*cya - dE  # Eq.(1)
        eps_xz = 2*tl + 2*th - 2*tl*cxa - 2*th*cya  # Eq.(2)
        eps_yz = 2*tl + 2*th - 2*th*cxa - 2*tl*cya  # Eq.(3)
        eps_h  = 2*td * sxa * sya  # Eq.(4)

        H0_orb = np.array([[eps_xy, 0,      0     ],
                            [0,      eps_xz, eps_h ],
                            [0,      eps_h,  eps_yz]], dtype=complex)
        H0 = np.kron(H0_orb, I2)

        # Rashba: orbital off-diagonal, scalar in spin space (⊗σ₀), Eq.(6)
        H_RSO_orb = dRSO * np.array([[0,       1j*sya, 1j*sxa],
                                      [-1j*sya, 0,      0     ],
                                      [-1j*sxa, 0,      0     ]], dtype=complex)
        H_RSO = np.kron(H_RSO_orb, I2)

        H = H0 + H_SO + H_RSO
        bands[i] = eigvalsh(H)

    return bands


def test_case_sub4_1():
    import pytest
    tl, th, td, dE, dSO, dRSO = 875.0, 40.0, 40.0, 47.0, 10.0, 5.0
    a = 0.394
    Nx, Ny, s = 2, 2, 1.0
    kspace_eigs, rs_eigs, max_dev = verify_kspace_realspace_agreement(
        Nx, Ny, tl, th, td, dE, dSO, dRSO, a, s)
    # BREAKS if: model uses wrong k-point formula or inconsistent SOC signs
    assert max_dev < 1e-8, \
        f'At s=1, max_deviation must be <1e-8 meV (machine precision), got {max_dev:.4e}'
    assert np.all(kspace_eigs[:-1] <= kspace_eigs[1:]), 'kspace_eigs must be sorted'
    assert np.all(rs_eigs[:-1] <= rs_eigs[1:]), 'rs_eigs must be sorted'
Test case 2 of 5
import numpy as np
from scipy.sparse import lil_matrix, csr_matrix
from scipy.sparse.linalg import eigsh
from scipy.linalg import eigvalsh
import pytest

def compute_kspace_bandstructure(k_path, tl, th, td, dE, dSO, dRSO, a):
    '''Compute band structure along a k-path for the LAO/STO 2DEG.

    Parameters:
      k_path : ndarray, shape (N_k, 2)
          Array of (kx, ky) values in units of 1/nm.
      tl : float
          Light-band hopping in meV.
      th : float
          Heavy-band hopping in meV.
      td : float
          Hybridisation hopping in meV.
      dE : float
          Crystal-field splitting in meV.
      dSO : float
          Atomic spin-orbit coupling in meV.
      dRSO : float
          Rashba spin-orbit coupling in meV.
      a : float
          Lattice constant in nm.

    Returns:
      bands : ndarray, shape (N_k, 6)
          Sorted eigenvalues at each k-point in meV.
    '''
    k_path = np.asarray(k_path)
    N_k = k_path.shape[0]
    bands = np.zeros((N_k, 6))

    sx = np.array([[0, 1], [1, 0]], dtype=complex)
    sy = np.array([[0, -1j], [1j, 0]], dtype=complex)
    sz = np.array([[1, 0], [0, -1]], dtype=complex)
    I2 = np.eye(2, dtype=complex)

    # H_SO = -(dSO/3)*block_matrix(Pauli sub-blocks), Eq.(5) of paper
    # Overall minus from t2g projection: L^t2g = -L^(ell=1)
    # Produces quartet at -dSO/3 and doublet at +2*dSO/3
    H_SO = np.zeros((6, 6), dtype=complex)
    H_SO[0:2, 2:4] = -(dSO / 3.0) * (1j * sx)
    H_SO[0:2, 4:6] = -(dSO / 3.0) * (-1j * sy)
    H_SO[2:4, 0:2] = -(dSO / 3.0) * (-1j * sx)
    H_SO[2:4, 4:6] = -(dSO / 3.0) * (1j * sz)
    H_SO[4:6, 0:2] = -(dSO / 3.0) * (1j * sy)
    H_SO[4:6, 2:4] = -(dSO / 3.0) * (-1j * sz)

    for i, (kx, ky) in enumerate(k_path):
        cxa = np.cos(kx * a)
        cya = np.cos(ky * a)
        sxa = np.sin(kx * a)
        sya = np.sin(ky * a)

        eps_xy = 4*tl - 2*tl*cxa - 2*tl*cya - dE  # Eq.(1)
        eps_xz = 2*tl + 2*th - 2*tl*cxa - 2*th*cya  # Eq.(2)
        eps_yz = 2*tl + 2*th - 2*th*cxa - 2*tl*cya  # Eq.(3)
        eps_h  = 2*td * sxa * sya  # Eq.(4)

        H0_orb = np.array([[eps_xy, 0,      0     ],
                            [0,      eps_xz, eps_h ],
                            [0,      eps_h,  eps_yz]], dtype=complex)
        H0 = np.kron(H0_orb, I2)

        # Rashba: orbital off-diagonal, scalar in spin space (⊗σ₀), Eq.(6)
        H_RSO_orb = dRSO * np.array([[0,       1j*sya, 1j*sxa],
                                      [-1j*sya, 0,      0     ],
                                      [-1j*sxa, 0,      0     ]], dtype=complex)
        H_RSO = np.kron(H_RSO_orb, I2)

        H = H0 + H_SO + H_RSO
        bands[i] = eigvalsh(H)

    return bands


def test_case_sub4_2():
    import pytest
    tl, th, td, dE, dSO, dRSO = 875.0, 40.0, 40.0, 47.0, 10.0, 5.0
    a = 0.394
    Nx, Ny, s = 2, 2, 1.0
    result = verify_kspace_realspace_agreement(Nx, Ny, tl, th, td, dE, dSO, dRSO, a, s)
    assert isinstance(result, tuple), f'Must return tuple, got {type(result)}'
    assert len(result) == 3, f'Tuple must have 3 elements, got {len(result)}'
    kspace_eigs, rs_eigs, max_dev = result
    assert isinstance(kspace_eigs, np.ndarray)
    assert isinstance(rs_eigs, np.ndarray)
    assert kspace_eigs.shape == (6 * Nx * Ny,), \
        f'kspace_eigs shape: expected ({6*Nx*Ny},), got {kspace_eigs.shape}'
    assert rs_eigs.shape == (6 * Nx * Ny,), \
        f'rs_eigs shape: expected ({6*Nx*Ny},), got {rs_eigs.shape}'
    assert isinstance(max_dev, float), f'max_deviation must be float, got {type(max_dev)}'
    assert max_dev >= 0.0, f'max_deviation must be non-negative, got {max_dev}'
Test case 3 of 5
import numpy as np
from scipy.sparse import lil_matrix, csr_matrix
from scipy.sparse.linalg import eigsh
from scipy.linalg import eigvalsh
import pytest

def compute_kspace_bandstructure(k_path, tl, th, td, dE, dSO, dRSO, a):
    '''Compute band structure along a k-path for the LAO/STO 2DEG.

    Parameters:
      k_path : ndarray, shape (N_k, 2)
          Array of (kx, ky) values in units of 1/nm.
      tl : float
          Light-band hopping in meV.
      th : float
          Heavy-band hopping in meV.
      td : float
          Hybridisation hopping in meV.
      dE : float
          Crystal-field splitting in meV.
      dSO : float
          Atomic spin-orbit coupling in meV.
      dRSO : float
          Rashba spin-orbit coupling in meV.
      a : float
          Lattice constant in nm.

    Returns:
      bands : ndarray, shape (N_k, 6)
          Sorted eigenvalues at each k-point in meV.
    '''
    k_path = np.asarray(k_path)
    N_k = k_path.shape[0]
    bands = np.zeros((N_k, 6))

    sx = np.array([[0, 1], [1, 0]], dtype=complex)
    sy = np.array([[0, -1j], [1j, 0]], dtype=complex)
    sz = np.array([[1, 0], [0, -1]], dtype=complex)
    I2 = np.eye(2, dtype=complex)

    # H_SO = -(dSO/3)*block_matrix(Pauli sub-blocks), Eq.(5) of paper
    # Overall minus from t2g projection: L^t2g = -L^(ell=1)
    # Produces quartet at -dSO/3 and doublet at +2*dSO/3
    H_SO = np.zeros((6, 6), dtype=complex)
    H_SO[0:2, 2:4] = -(dSO / 3.0) * (1j * sx)
    H_SO[0:2, 4:6] = -(dSO / 3.0) * (-1j * sy)
    H_SO[2:4, 0:2] = -(dSO / 3.0) * (-1j * sx)
    H_SO[2:4, 4:6] = -(dSO / 3.0) * (1j * sz)
    H_SO[4:6, 0:2] = -(dSO / 3.0) * (1j * sy)
    H_SO[4:6, 2:4] = -(dSO / 3.0) * (-1j * sz)

    for i, (kx, ky) in enumerate(k_path):
        cxa = np.cos(kx * a)
        cya = np.cos(ky * a)
        sxa = np.sin(kx * a)
        sya = np.sin(ky * a)

        eps_xy = 4*tl - 2*tl*cxa - 2*tl*cya - dE  # Eq.(1)
        eps_xz = 2*tl + 2*th - 2*tl*cxa - 2*th*cya  # Eq.(2)
        eps_yz = 2*tl + 2*th - 2*th*cxa - 2*tl*cya  # Eq.(3)
        eps_h  = 2*td * sxa * sya  # Eq.(4)

        H0_orb = np.array([[eps_xy, 0,      0     ],
                            [0,      eps_xz, eps_h ],
                            [0,      eps_h,  eps_yz]], dtype=complex)
        H0 = np.kron(H0_orb, I2)

        # Rashba: orbital off-diagonal, scalar in spin space (⊗σ₀), Eq.(6)
        H_RSO_orb = dRSO * np.array([[0,       1j*sya, 1j*sxa],
                                      [-1j*sya, 0,      0     ],
                                      [-1j*sxa, 0,      0     ]], dtype=complex)
        H_RSO = np.kron(H_RSO_orb, I2)

        H = H0 + H_SO + H_RSO
        bands[i] = eigvalsh(H)

    return bands


def test_case_sub4_3():
    import pytest
    tl, th, td, dE, dSO, dRSO = 875.0, 40.0, 40.0, 47.0, 10.0, 5.0
    a = 0.394
    Nx, Ny, s = 3, 3, 1.0
    kspace_eigs, rs_eigs, max_dev = verify_kspace_realspace_agreement(
        Nx, Ny, tl, th, td, dE, dSO, dRSO, a, s)
    # BREAKS if: any sign convention is inconsistent between k-space and real-space
    assert kspace_eigs.shape == (54,), f'Expected 54 eigenvalues, got {kspace_eigs.shape}'
    assert rs_eigs.shape == (54,), f'Expected 54 eigenvalues, got {rs_eigs.shape}'
    assert max_dev < 1e-8, \
        f'Bloch theorem: max_dev must be <1e-8 meV at s=1, got {max_dev:.4e}'
    assert kspace_eigs[0] < 0, f'Lowest eigenvalue must be negative, got {kspace_eigs[0]}'
Test case 4 of 5
import numpy as np
from scipy.sparse import lil_matrix, csr_matrix
from scipy.sparse.linalg import eigsh
from scipy.linalg import eigvalsh
import pytest

def compute_kspace_bandstructure(k_path, tl, th, td, dE, dSO, dRSO, a):
    '''Compute band structure along a k-path for the LAO/STO 2DEG.

    Parameters:
      k_path : ndarray, shape (N_k, 2)
          Array of (kx, ky) values in units of 1/nm.
      tl : float
          Light-band hopping in meV.
      th : float
          Heavy-band hopping in meV.
      td : float
          Hybridisation hopping in meV.
      dE : float
          Crystal-field splitting in meV.
      dSO : float
          Atomic spin-orbit coupling in meV.
      dRSO : float
          Rashba spin-orbit coupling in meV.
      a : float
          Lattice constant in nm.

    Returns:
      bands : ndarray, shape (N_k, 6)
          Sorted eigenvalues at each k-point in meV.
    '''
    k_path = np.asarray(k_path)
    N_k = k_path.shape[0]
    bands = np.zeros((N_k, 6))

    sx = np.array([[0, 1], [1, 0]], dtype=complex)
    sy = np.array([[0, -1j], [1j, 0]], dtype=complex)
    sz = np.array([[1, 0], [0, -1]], dtype=complex)
    I2 = np.eye(2, dtype=complex)

    # H_SO = -(dSO/3)*block_matrix(Pauli sub-blocks), Eq.(5) of paper
    # Overall minus from t2g projection: L^t2g = -L^(ell=1)
    # Produces quartet at -dSO/3 and doublet at +2*dSO/3
    H_SO = np.zeros((6, 6), dtype=complex)
    H_SO[0:2, 2:4] = -(dSO / 3.0) * (1j * sx)
    H_SO[0:2, 4:6] = -(dSO / 3.0) * (-1j * sy)
    H_SO[2:4, 0:2] = -(dSO / 3.0) * (-1j * sx)
    H_SO[2:4, 4:6] = -(dSO / 3.0) * (1j * sz)
    H_SO[4:6, 0:2] = -(dSO / 3.0) * (1j * sy)
    H_SO[4:6, 2:4] = -(dSO / 3.0) * (-1j * sz)

    for i, (kx, ky) in enumerate(k_path):
        cxa = np.cos(kx * a)
        cya = np.cos(ky * a)
        sxa = np.sin(kx * a)
        sya = np.sin(ky * a)

        eps_xy = 4*tl - 2*tl*cxa - 2*tl*cya - dE  # Eq.(1)
        eps_xz = 2*tl + 2*th - 2*tl*cxa - 2*th*cya  # Eq.(2)
        eps_yz = 2*tl + 2*th - 2*th*cxa - 2*tl*cya  # Eq.(3)
        eps_h  = 2*td * sxa * sya  # Eq.(4)

        H0_orb = np.array([[eps_xy, 0,      0     ],
                            [0,      eps_xz, eps_h ],
                            [0,      eps_h,  eps_yz]], dtype=complex)
        H0 = np.kron(H0_orb, I2)

        # Rashba: orbital off-diagonal, scalar in spin space (⊗σ₀), Eq.(6)
        H_RSO_orb = dRSO * np.array([[0,       1j*sya, 1j*sxa],
                                      [-1j*sya, 0,      0     ],
                                      [-1j*sxa, 0,      0     ]], dtype=complex)
        H_RSO = np.kron(H_RSO_orb, I2)

        H = H0 + H_SO + H_RSO
        bands[i] = eigvalsh(H)

    return bands


def test_case_sub4_4():
    import pytest
    tl, th, td, dE, dSO, dRSO = 875.0, 40.0, 40.0, 47.0, 10.0, 5.0
    a = 0.394
    Nx, Ny, s = 2, 2, 1.0
    kspace_eigs, rs_eigs, max_dev = verify_kspace_realspace_agreement(
        Nx, Ny, tl, th, td, dE, dSO, dRSO, a, s)
    expected_max_dev = float(np.max(np.abs(kspace_eigs - rs_eigs)))
    # BREAKS if: model computes max_dev from unsorted arrays or uses wrong metric
    assert np.isclose(max_dev, expected_max_dev, atol=1e-10), \
        f'max_deviation must equal max|sorted_kspace - sorted_rs|. Expected {expected_max_dev}, got {max_dev}'
    assert max_dev >= 0.0
    assert isinstance(max_dev, float)
Test case 5 of 5
import numpy as np
from scipy.sparse import lil_matrix, csr_matrix
from scipy.sparse.linalg import eigsh
from scipy.linalg import eigvalsh
import pytest

def compute_kspace_bandstructure(k_path, tl, th, td, dE, dSO, dRSO, a):
    '''Compute band structure along a k-path for the LAO/STO 2DEG.

    Parameters:
      k_path : ndarray, shape (N_k, 2)
          Array of (kx, ky) values in units of 1/nm.
      tl : float
          Light-band hopping in meV.
      th : float
          Heavy-band hopping in meV.
      td : float
          Hybridisation hopping in meV.
      dE : float
          Crystal-field splitting in meV.
      dSO : float
          Atomic spin-orbit coupling in meV.
      dRSO : float
          Rashba spin-orbit coupling in meV.
      a : float
          Lattice constant in nm.

    Returns:
      bands : ndarray, shape (N_k, 6)
          Sorted eigenvalues at each k-point in meV.
    '''
    k_path = np.asarray(k_path)
    N_k = k_path.shape[0]
    bands = np.zeros((N_k, 6))

    sx = np.array([[0, 1], [1, 0]], dtype=complex)
    sy = np.array([[0, -1j], [1j, 0]], dtype=complex)
    sz = np.array([[1, 0], [0, -1]], dtype=complex)
    I2 = np.eye(2, dtype=complex)

    # H_SO = -(dSO/3)*block_matrix(Pauli sub-blocks), Eq.(5) of paper
    # Overall minus from t2g projection: L^t2g = -L^(ell=1)
    # Produces quartet at -dSO/3 and doublet at +2*dSO/3
    H_SO = np.zeros((6, 6), dtype=complex)
    H_SO[0:2, 2:4] = -(dSO / 3.0) * (1j * sx)
    H_SO[0:2, 4:6] = -(dSO / 3.0) * (-1j * sy)
    H_SO[2:4, 0:2] = -(dSO / 3.0) * (-1j * sx)
    H_SO[2:4, 4:6] = -(dSO / 3.0) * (1j * sz)
    H_SO[4:6, 0:2] = -(dSO / 3.0) * (1j * sy)
    H_SO[4:6, 2:4] = -(dSO / 3.0) * (-1j * sz)

    for i, (kx, ky) in enumerate(k_path):
        cxa = np.cos(kx * a)
        cya = np.cos(ky * a)
        sxa = np.sin(kx * a)
        sya = np.sin(ky * a)

        eps_xy = 4*tl - 2*tl*cxa - 2*tl*cya - dE  # Eq.(1)
        eps_xz = 2*tl + 2*th - 2*tl*cxa - 2*th*cya  # Eq.(2)
        eps_yz = 2*tl + 2*th - 2*th*cxa - 2*tl*cya  # Eq.(3)
        eps_h  = 2*td * sxa * sya  # Eq.(4)

        H0_orb = np.array([[eps_xy, 0,      0     ],
                            [0,      eps_xz, eps_h ],
                            [0,      eps_h,  eps_yz]], dtype=complex)
        H0 = np.kron(H0_orb, I2)

        # Rashba: orbital off-diagonal, scalar in spin space (⊗σ₀), Eq.(6)
        H_RSO_orb = dRSO * np.array([[0,       1j*sya, 1j*sxa],
                                      [-1j*sya, 0,      0     ],
                                      [-1j*sxa, 0,      0     ]], dtype=complex)
        H_RSO = np.kron(H_RSO_orb, I2)

        H = H0 + H_SO + H_RSO
        bands[i] = eigvalsh(H)

    return bands


def test_case_sub4_constraint_error_handling():
    import pytest
    # Constraint coverage: error_handling
    # verify_kspace_realspace_agreement raises ValueError when s <= 0
    with pytest.raises(ValueError):
        verify_kspace_realspace_agreement(2, 2, 875.0, 40.0, 40.0, 47.0, 10.0, 5.0, 0.394, -1.0)
    # Also verify that a valid call succeeds and returns correct types
    kspace_eigs, rs_eigs, max_dev = verify_kspace_realspace_agreement(
        2, 2, 875.0, 40.0, 40.0, 47.0, 10.0, 5.0, 0.394, 1.0)
    assert isinstance(max_dev, float)
    assert max_dev >= 0.0
Solution
import numpy as np
from scipy.sparse import lil_matrix
from scipy.linalg import eigvalsh

def compute_kspace_bandstructure(k_path, tl, th, td, dE, dSO, dRSO, a):
    '''Compute band structure along a k-path for the LAO/STO 2DEG.

    Parameters:
      k_path : ndarray, shape (N_k, 2)
          Array of (kx, ky) values in units of 1/nm.
      tl : float
          Light-band hopping in meV.
      th : float
          Heavy-band hopping in meV.
      td : float
          Hybridisation hopping in meV.
      dE : float
          Crystal-field splitting in meV.
      dSO : float
          Atomic spin-orbit coupling in meV.
      dRSO : float
          Rashba spin-orbit coupling in meV.
      a : float
          Lattice constant in nm.

    Returns:
      bands : ndarray, shape (N_k, 6)
          Sorted eigenvalues at each k-point in meV.
    '''
    k_path = np.asarray(k_path)
    N_k = k_path.shape[0]
    bands = np.zeros((N_k, 6))

    sx = np.array([[0, 1], [1, 0]], dtype=complex)
    sy = np.array([[0, -1j], [1j, 0]], dtype=complex)
    sz = np.array([[1, 0], [0, -1]], dtype=complex)
    I2 = np.eye(2, dtype=complex)

    # H_SO = -(dSO/3)*block_matrix(Pauli sub-blocks), Eq.(5) of paper
    # Overall minus from t2g projection: L^t2g = -L^(ell=1)
    # Produces quartet at -dSO/3 and doublet at +2*dSO/3
    H_SO = np.zeros((6, 6), dtype=complex)
    H_SO[0:2, 2:4] = -(dSO / 3.0) * (1j * sx)
    H_SO[0:2, 4:6] = -(dSO / 3.0) * (-1j * sy)
    H_SO[2:4, 0:2] = -(dSO / 3.0) * (-1j * sx)
    H_SO[2:4, 4:6] = -(dSO / 3.0) * (1j * sz)
    H_SO[4:6, 0:2] = -(dSO / 3.0) * (1j * sy)
    H_SO[4:6, 2:4] = -(dSO / 3.0) * (-1j * sz)

    for i, (kx, ky) in enumerate(k_path):
        cxa = np.cos(kx * a)
        cya = np.cos(ky * a)
        sxa = np.sin(kx * a)
        sya = np.sin(ky * a)

        eps_xy = 4*tl - 2*tl*cxa - 2*tl*cya - dE  # Eq.(1)
        eps_xz = 2*tl + 2*th - 2*tl*cxa - 2*th*cya  # Eq.(2)
        eps_yz = 2*tl + 2*th - 2*th*cxa - 2*tl*cya  # Eq.(3)
        eps_h  = 2*td * sxa * sya  # Eq.(4)

        H0_orb = np.array([[eps_xy, 0,      0     ],
                            [0,      eps_xz, eps_h ],
                            [0,      eps_h,  eps_yz]], dtype=complex)
        H0 = np.kron(H0_orb, I2)

        # Rashba: orbital off-diagonal, scalar in spin space (⊗σ₀), Eq.(6)
        H_RSO_orb = dRSO * np.array([[0,       1j*sya, 1j*sxa],
                                      [-1j*sya, 0,      0     ],
                                      [-1j*sxa, 0,      0     ]], dtype=complex)
        H_RSO = np.kron(H_RSO_orb, I2)

        H = H0 + H_SO + H_RSO
        bands[i] = eigvalsh(H)

    return bands

def build_scaled_hopping_blocks(tl, th, td, dRSO, s):
    '''Build the scaled real-space hopping blocks.

    Parameters:
      tl : float
          Light-band hopping in meV (unscaled).
      th : float
          Heavy-band hopping in meV (unscaled).
      td : float
          Hybridisation hopping in meV (unscaled).
      dRSO : float
          Rashba coupling in meV (unscaled).
      s : float
          Scaling factor (s=1 means no scaling).

    Returns:
      blocks : dict
          Dictionary with keys 'Hx', 'Hy', 'Hmix', each mapping
          to a complex ndarray of shape (6, 6).
    '''
    tl_p   = tl   / s**2  # kinetic scales as 1/s^2, Eq.(9)
    th_p   = th   / s**2
    td_p   = td   / s**2
    dRSO_p = dRSO / s     # Rashba scales as 1/s, Eq.(9)

    I2 = np.eye(2, dtype=complex)

    # Hx: along x, dxy and dxz have light hopping, dyz has heavy hopping; Eq.(10)
    Hx_orb_kin = np.diag([-tl_p, -tl_p, -th_p]).astype(complex)
    Hx_orb_rso = (dRSO_p / 2.0) * np.array([[0, 0, -1],
                                              [0, 0,  0],
                                              [1, 0,  0]], dtype=complex)
    Hx = np.kron(Hx_orb_kin + Hx_orb_rso, I2)

    # Hy: along y, dxy and dyz have light hopping, dxz has heavy hopping; Eq.(11)
    Hy_orb_kin = np.diag([-tl_p, -th_p, -tl_p]).astype(complex)
    Hy_orb_rso = (dRSO_p / 2.0) * np.array([[ 0, -1, 0],
                                              [ 1,  0, 0],
                                              [ 0,  0, 0]], dtype=complex)
    Hy = np.kron(Hy_orb_kin + Hy_orb_rso, I2)

    # Hmix: diagonal hopping couples only dxz↔dyz; factor 1/2 from sin discrete FT; Eq.(12)
    Hmix_orb = (td_p / 2.0) * np.array([[0, 0, 0],
                                          [0, 0, 1],
                                          [0, 1, 0]], dtype=complex)
    Hmix = np.kron(Hmix_orb, I2)

    return {'Hx': Hx, 'Hy': Hy, 'Hmix': Hmix}

def build_and_diagonalise_realspace(Nx, Ny, tl, th, td, dE, dSO, dRSO, s):
    '''Build and diagonalise the real-space Hamiltonian.

    Parameters:
      Nx : int
          Number of sites in x-direction.
      Ny : int
          Number of sites in y-direction.
      tl : float
          Light-band hopping in meV (unscaled).
      th : float
          Heavy-band hopping in meV (unscaled).
      td : float
          Hybridisation hopping in meV (unscaled).
      dE : float
          Crystal-field splitting in meV.
      dSO : float
          Atomic spin-orbit coupling in meV.
      dRSO : float
          Rashba coupling in meV (unscaled).
      s : float
          Scaling factor.

    Returns:
      eigs : ndarray, shape (6*Nx*Ny,)
          All eigenvalues sorted in ascending order, in meV.
    '''
    tl_p = tl / s**2
    th_p = th / s**2

    blocks = build_scaled_hopping_blocks(tl, th, td, dRSO, s)
    Hx = blocks['Hx']
    Hy = blocks['Hy']
    Hmix = blocks['Hmix']

    sx = np.array([[0, 1], [1, 0]], dtype=complex)
    sy = np.array([[0, -1j], [1j, 0]], dtype=complex)
    sz = np.array([[1, 0], [0, -1]], dtype=complex)
    I2 = np.eye(2, dtype=complex)

    H_SO = np.zeros((6, 6), dtype=complex)
    H_SO[0:2, 2:4] = -(dSO / 3.0) * (1j * sx)
    H_SO[0:2, 4:6] = -(dSO / 3.0) * (-1j * sy)
    H_SO[2:4, 0:2] = -(dSO / 3.0) * (-1j * sx)
    H_SO[2:4, 4:6] = -(dSO / 3.0) * (1j * sz)
    H_SO[4:6, 0:2] = -(dSO / 3.0) * (1j * sy)
    H_SO[4:6, 2:4] = -(dSO / 3.0) * (-1j * sz)

    H_onsite_kin_orb = np.diag([4*tl_p - dE, 2*tl_p + 2*th_p, 2*tl_p + 2*th_p]).astype(complex)
    H_on = np.kron(H_onsite_kin_orb, I2) + H_SO

    D = 6 * Nx * Ny
    H = lil_matrix((D, D), dtype=complex)

    for mu in range(Nx):
        for nu in range(Ny):
            n = mu * Ny + nu
            H[6*n : 6*(n+1), 6*n : 6*(n+1)] = H_on

            # Hopping terms (forward directions)
            n_xp1 = ((mu + 1) % Nx) * Ny + nu
            H[6*n : 6*(n+1), 6*n_xp1 : 6*(n_xp1+1)] += Hx
            H[6*n_xp1 : 6*(n_xp1+1), 6*n : 6*(n+1)] += Hx.conj().T

            n_yp1 = mu * Ny + ((nu + 1) % Ny)
            H[6*n : 6*(n+1), 6*n_yp1 : 6*(n_yp1+1)] += Hy
            H[6*n_yp1 : 6*(n_yp1+1), 6*n : 6*(n+1)] += Hy.conj().T

            n_xp1_ym1 = ((mu + 1) % Nx) * Ny + ((nu - 1 + Ny) % Ny)
            H[6*n : 6*(n+1), 6*n_xp1_ym1 : 6*(n_xp1_ym1+1)] += Hmix
            H[6*n_xp1_ym1 : 6*(n_xp1_ym1+1), 6*n : 6*(n+1)] += Hmix.conj().T

            n_xp1_yp1 = ((mu + 1) % Nx) * Ny + ((nu + 1) % Ny)
            H[6*n : 6*(n+1), 6*n_xp1_yp1 : 6*(n_xp1_yp1+1)] -= Hmix
            H[6*n_xp1_yp1 : 6*(n_xp1_yp1+1), 6*n : 6*(n+1)] -= Hmix.conj().T

    eigs = eigvalsh(H.toarray())
    return eigs

def verify_kspace_realspace_agreement(Nx, Ny, tl, th, td, dE, dSO, dRSO, a, s):
    '''Compare k-space and real-space spectra.

    Parameters:
      Nx : int
          Number of sites in x-direction.
      Ny : int
          Number of sites in y-direction.
      tl : float
          Light-band hopping in meV.
      th : float
          Heavy-band hopping in meV.
      td : float
          Hybridisation hopping in meV.
      dE : float
          Crystal-field splitting in meV.
      dSO : float
          Atomic spin-orbit coupling in meV.
      dRSO : float
          Rashba coupling in meV.
      a : float
          Lattice constant in nm.
      s : float
          Scaling factor.

    Returns:
      kspace_eigs : ndarray, shape (6*Nx*Ny,)
          Sorted eigenvalues from k-space at commensurate k-points, in meV.
      realspace_eigs : ndarray, shape (6*Nx*Ny,)
          Sorted eigenvalues from real-space Hamiltonian, in meV.
      max_deviation : float
          Maximum absolute difference in meV.
    '''
    if s <= 0:
        raise ValueError(f'Scaling factor s must be positive, got s={s}')

    # K-space calculation
    mu_vals = np.arange(Nx)
    nu_vals = np.arange(Ny)
    k_grid = np.array(np.meshgrid(mu_vals, nu_vals)).T.reshape(-1, 2)
    
    kx_vals = 2 * np.pi * k_grid[:, 0] / (Nx * s * a)
    ky_vals = 2 * np.pi * k_grid[:, 1] / (Ny * s * a)
    commensurate_k_path = np.column_stack([kx_vals, ky_vals])
    
    k_bands = compute_kspace_bandstructure(commensurate_k_path, tl, th, td, dE, dSO, dRSO, a)
    kspace_eigs = np.sort(k_bands.flatten())
    
    # Real-space calculation
    realspace_eigs = build_and_diagonalise_realspace(Nx, Ny, tl, th, td, dE, dSO, dRSO, s)
    
    # Comparison
    max_deviation = float(np.max(np.abs(kspace_eigs - realspace_eigs)))
    
    return kspace_eigs, realspace_eigs, max_deviation