10. Eigenchannel analysis

The transmission $T(E)$ obtained in Chapter 7 is only the total transmission probability per energy; it does not tell us where in the junction the electrons flow or where they are reflected. In this chapter we decompose the total transmission into mutually independent conduction channels (eigenchannels) and render the scattering wavefunction of each channel as a real-space cube, so the current path can be inspected directly. The key point is that this is post-processing, not a recalculation — all that is needed is the trans.TSHS already saved in Chapter 6 and the electrode Electrode.TSHS.

Learning objectives

Understand the theoretical basis for decomposing the transmission into per-channel eigenvalues $\tau_i$ .
Construct Green function → spectral function → eigenchannels by sisl post-processing.
Save channel wavefunctions as three kinds of cube (Re / Im / $|\psi|^2$ ) and visualize them in VESTA.
Learn the criteria for choosing the energy at which to extract eigenchannels, and how to self-verify a poor choice.

Theory: decomposing the total transmission into independent channels

In NEGF the transmission is given by

T(E) = \mathrm{Tr}\!\left[\Gamma_L\, G^r\, \Gamma_R\, G^a\right]

(Chapter 4). This trace can be rewritten as the trace of the Hermitian matrix $t^\dagger t$ (with $t \equiv \Gamma_R^{1/2} G^r \Gamma_L^{1/2}$ ), so it can be diagonalized by a suitable unitary rotation, and that eigenbasis is precisely the set of eigenchannels:

T(E) = \sum_i \tau_i(E), \qquad 0 \le \tau_i \le 1 .

Each $\tau_i$ is "how open the $i$ -th channel is", and a fully open channel ( $\tau_i = 1$ ) carries exactly one conductance quantum $G_0 = 2e^2/h$ . The channel decomposition is not merely a mathematical rewriting; it carries physical information:

Channel counting — the integer plateaus of $T(E)$ in a ballistic 1D conductor equal the number of open channels. A uniform cumulene chain has two degenerate $\pi_x$ , $\pi_y$ bands near $E_F$ , so $T(E_F) \approx 2$ , i.e. $\tau_0 \approx \tau_1 \approx 1$ is expected.
Localization of the scattering — with an impurity present, only the $\tau$ of a particular channel is suppressed, and the $|\psi|^2$ of that channel bending at the impurity site shows "where the reflection happens".
Shot noise — the Fano factor is determined by $\sum_i \tau_i(1-\tau_i) / \sum_i \tau_i$ , so for the same $T$ the noise differs depending on the channel composition.

The construction procedure follows the Paulsson–Brandbyge method. Diagonalizing the left spectral function

A_L(E) = G^r\, \Gamma_L\, G^a

which contains all scattering states injected from the left electrode, yields the injected-state basis; rotating it with $\Gamma_R$ makes the eigenvalues $\tau_i$ and the eigenvectors the channel wavefunctions. Since the PAO basis of SIESTA is non-orthogonal, a Löwdin orthogonalization with the overlap $S$ (the transformation that sandwiches $S^{1/2}$ on both sides) must be performed first.

Method: sisl post-processing, not a recalculation

Three inputs are required.

File	Source	Role
`trans.TSHS`	Chapter 6 device 0 V calculation	Device Hamiltonian + overlap
`Electrode.TSHS`	Chapter 5 electrode calculation	Generating the semi-infinite electrode self-energy
`*.ion.xml`	SIESTA calculation directory	Radial information of the basis orbitals (for cube generation)

TBtrans has a TBT.T.Eig option that prints channel eigenvalues as text, but it is not used. It inflates the output nc file to several GB and causes out-of-memory (OOM) failures, and it does not give the wavefunctions anyway. Obtaining both the channel eigenvalues and the wavefunctions by the post-processing above is the standard.

note

Cube generation requires the radial function of each orbital, which sisl reads from *.ion.xml. Plain *.ion files are not supported by sisl, so keep the *.ion.xml files in the SIESTA calculation directory rather than deleting them. It is also safer to read the Hamiltonian via input.fdf than to read the TSHS directly — the fdf sile finds the *.ion.xml files in the same directory and attaches the basis information to the geometry.

Energy selection — an energy that matches the analysis question

Eigenchannels are defined at any energy. For impurity/molecular junctions one usually selects a resonance or a transmission dip near $E_F$ to analyze the scattering states directly relevant to transport. The largest peak in the PDOS is not always the appropriate choice, so the selected energy must be presented together with the purpose of the analysis.

Always self-verify after extraction:

For a localized resonance, a single channel may dominate — check whether $\tau_0$ accounts for most of $T$ .
Contributions from several channels and $T>1$ can be physically normal for multi-mode electrodes. This is not a criterion for declaring an error; interpret it together with the bands, PDOS, and channel symmetry at the selected energy.
At energies where $T \approx 0$ (effectively disconnected) there is no transmitting channel at all. In that case, render the left-injected scattering state (an eigenstate of $A_L$ ) instead of the $\tau_i$ eigenvectors to visualize "it is blocked here".

Code

The script below is in the repository at code/ch10-eigenchannel/eigenchannel.py. The sisl API differs slightly between versions (this targets 0.14+), so if the function signatures do not match, check the official sisl documentation.

code/ch10-eigenchannel/eigenchannel.py
import numpy as np
import scipy.linalg as sla
import sisl
from sisl.physics.electron import wavefunction

# ---------------------------- user settings ----------------------------
FDF = "input.fdf"                 # device calculation input (for loading geometry + basis)
ELEC_TSHS = "../electrode/Electrode.TSHS"
ENERGIES = {                      # eV, relative to E - E_F
    "EF": 0.0,
    "dip": -0.85,                 # replace with the dip energy read directly from the Chapter 7 T(E)
}
ETA = 1e-3                        # eV. imaginary part of the Green function
NA_ELEC = 4                       # number of atoms corresponding to one electrode cell at each device end (C4)
NCHAN = 2                         # number of channels to save as cubes
GRID_SPACING = 0.2                # Ang
# ---------------------------------------------------------------------

# 1) read the Hamiltonian — reading via the fdf attaches the *.ion.xml basis to the geometry
H = sisl.get_sile(FDF).read_hamiltonian()      # automatically locates the TSHS of SystemLabel (trans)
H_el = sisl.get_sile(ELEC_TSHS).read_hamiltonian()

geom = H.geometry
no = H.no

# 2) Gamma-point dense matrices (1D chain: no transverse k averaging needed)
Hd = H.Hk(format="array")
Sd = H.Sk(format="array")

# 3) semi-infinite electrode self-energies — transport axis z: left -C, right +C
SE_L = sisl.physics.RecursiveSI(H_el, "-C")
SE_R = sisl.physics.RecursiveSI(H_el, "+C")

idx_L = geom.a2o(np.arange(NA_ELEC), all=True)
idx_R = geom.a2o(np.arange(geom.na - NA_ELEC, geom.na), all=True)

def broadening(sig):
    return 1j * (sig - sig.conj().T)

for label, E in ENERGIES.items():
    Z = E + 1j * ETA

    # for the energy-argument convention of self_energy (complex E vs eta keyword) see the sisl docs
    seL = SE_L.self_energy(Z)
    seR = SE_R.self_energy(Z)

    SigL = np.zeros((no, no), dtype=complex)
    SigR = np.zeros((no, no), dtype=complex)
    SigL[np.ix_(idx_L, idx_L)] = seL
    SigR[np.ix_(idx_R, idx_R)] = seR

    GamL = broadening(SigL)
    GamR = broadening(SigR)

    # 4) Green function and left spectral function
    G = np.linalg.inv(Z * Sd - Hd - SigL - SigR)
    A_L = G @ GamL @ G.conj().T
    T_ref = np.trace(GamR @ A_L).real          # total transmission (for comparison)

    # 5) diagonalize A_L after Loewdin orthogonalization (Paulsson-Brandbyge)
    S12 = sla.sqrtm(Sd).real
    S12i = np.linalg.inv(S12)
    A_bar = S12 @ A_L @ S12
    GamR_bar = S12i @ GamR @ S12i

    lam, U = np.linalg.eigh(A_bar)
    lam = lam.clip(min=0.0)

    # scale the injected states, then rotate with Gamma_R -> tau_i
    Ut = U * np.sqrt(lam / (2.0 * np.pi))
    M = 2.0 * np.pi * (Ut.conj().T @ GamR_bar @ Ut)
    tau, W = np.linalg.eigh(M)
    order = np.argsort(tau)[::-1]
    tau, W = tau[order], W[:, order]

    print(f"[{label}]  E-E_F = {E:+.3f} eV   T = {T_ref:.4f}   "
          f"sum(tau) = {tau.sum():.4f}")
    print("        tau_i =", np.round(tau[:6], 4))

    # 6) PAO coefficients of the channel wavefunctions (back to the non-orthogonal basis)
    psi = S12i @ (Ut @ W)

    # 7) three kinds of real-space cube (Re / Im / |psi|^2)
    for i in range(NCHAN):
        c = psi[:, i]
        grid = sisl.Grid(GRID_SPACING, geometry=geom, dtype=np.complex128)
        wavefunction(c, grid, geometry=geom)   # requires *.ion.xml

        for kind, data in (("Re", grid.grid.real),
                           ("Im", grid.grid.imag),
                           ("psi2", np.abs(grid.grid) ** 2)):
            out = sisl.Grid(GRID_SPACING, geometry=geom, dtype=np.float64)
            out.grid[:] = data
            out.write(f"EC_{i}_{kind}_{label}.cube")

To summarize what the script does: it reads the two TSHS files, builds the self-energies from the semi-infinite repetition of the electrode Hamiltonian (RecursiveSI), embeds them in the orbital blocks at the two device ends, and then constructs the Green function and $A_L$ directly as dense matrices. For a 1D chain the matrices are small (a few hundred orbitals), so it runs instantly even on a laptop. Two verification points are built in — sum(tau) must agree with the T from the trace formula, and $|\psi|^2 = \mathrm{Re}^2 + \mathrm{Im}^2$ must hold on the grid.

Running

cd device_0bias            # the directory containing trans.TSHS, input.fdf, and *.ion.xml
python eigenchannel.py

Example output (the numbers depend on the basis and pseudopotentials):

stdout (example)
[EF]   E-E_F = +0.000 eV   T = 1.5321   sum(tau) = 1.5321
        tau_i = [0.9812 0.5509 0.      0.      0.      0.    ]
[dip]  E-E_F = -0.850 eV   T = 0.3130   sum(tau) = 0.3130
        tau_i = [0.2954 0.0176 0.      0.      0.      0.    ]

Two channels × three cube kinds = six cube files are produced per energy: EC_0_Re_EF.cube, EC_0_Im_EF.cube, EC_0_psi2_EF.cube, ... (naming rule: EC_channelindex_kind_energylabel.cube).

Analysis: how to read the cubes

Open EC_0_psi2_EF.cube in VESTA and lower the isosurface level to a few percent of the maximum to see the spatial distribution of the channel. There is a reason all three kinds must be viewed together:

$|\psi|^2$ — the path the current takes. Where the amplitude drops sharply is the reflection point.
Re, Im — a scattering state is a complex propagating wave. A propagating wave appears as a pair with Re and Im shifted by a quarter wavelength, and when reflection is strong enough to form a standing wave the nodes become fixed in space. Looking at $|\psi|^2$ alone, without Re/Im, makes it impossible to distinguish propagating from standing.

The expected picture for the C19N example is as follows. At $E_F$ one of the two $\pi$ channels passes the N site relatively well, and at the dip energy a standing-wave pattern forms on the injection side (left) and the amplitude drops sharply across the N site — reflection at the N impurity is directly visible in real space. The contrast with a uniform C20 chain (Exercise 1) makes the difference clear.

Exercises

Run the same script on the pristine C20 device with N reverted to C. Check whether two channels with $\tau_0 \approx \tau_1 \approx 1$ appear near $E_F$ , and whether the cubes spread uniformly over the whole chain.
From the impurity PDOS of C19N (the PDOS extraction of Chapter 7), find the peak closest to $E_F$ and the largest peak, extract eigenchannels at both energies, and determine which one is a localized resonance from the number of channels and $\sum\tau_i$ .
How do $\tau_0$ and the cube shapes change when ETA is raised from $10^{-3}$ to $10^{-2}$ eV? Read η broadening artifact and interpret the result.

References

M. Paulsson and M. Brandbyge, "Transmission eigenchannels from nonequilibrium Green's functions", Phys. Rev. B 76, 115117 (2007). DOI 10.1103/PhysRevB.76.115117
Official sisl documentation: https://zerothi.github.io/sisl/

Learning objectives​

Theory: decomposing the total transmission into independent channels​

Method: sisl post-processing, not a recalculation​

Code​

Running​

Analysis: how to read the cubes​

Exercises​

References​