Linear Scaling Quantum Transport Method

• Linear Scaling Quantum Transport Method is defined by computational approaches that use stochastic trace approximations and polynomial expansions to achieve O(N) scaling.
• It enables the simulation of transport regimes from ballistic to localization with techniques such as Chebyshev expansions and GPU-accelerated sparse matrix operations.
• The method integrates formulations like Kubo-Greenwood and Einstein to accurately compute electronic, spin, and topological transport properties in large-scale materials.

The linear scaling quantum transport method encompasses a family of computational approaches and concrete algorithms for calculating quantum transport properties of large, complex systems with computational complexity scaling linearly or polylogarithmically with system size. These methodologies circumvent the cubic (or worse) scaling of conventional diagonalization and Green's function approaches by leveraging stochastic trace approximations, iterative polynomial expansions, sparse matrix techniques, and more recently quantum algorithmic primitives. Their design enables the paper of materials and devices containing millions of degrees of freedom, capturing regimes from ballistic to diffusive and localization. The approaches are unified by the principle that the dominant computational steps—trace evaluations, time evolutions, spectral projections—require at most $\mathcal{O}(N)$ operations, where $N$ is the number of orbitals or sites.

1. Stochastic Trace and Polynomial Expansion Techniques

The central algorithmic advance is the replacement of explicit basis summations by stochastic trace estimation. For an operator $\hat{A}$, the trace is approximated as

$$\operatorname{Tr}[\hat{A}] \approx \langle \phi | \hat{A} | \phi \rangle,$$

where $|\phi\rangle$ is a normalized random-phase vector. The error decays as $1/\sqrt{N N_r}$ for $N_r$ independent samples, making even a single or few samples sufficient for extensive systems (Fan et al., 2018).

Functions of the Hamiltonian, such as the time evolution operator $U(t) = e^{-i\hat{H}t/\hbar}$ or spectral projectors $\delta(E-\hat{H})$, are efficiently computed using Chebyshev polynomial expansions:

$$U(t) \approx \sum_{m=0}^{M} \bar{u}_m(t)\, T_m(\tilde{\hat{H}})$$

with recurrence relations for the Chebyshev polynomials $T_m$ and coefficients $\bar{u}_m$. This expansion is critical for both the time evolution of quantum states and for implementing the kernel polynomial method for spectral resolution. The sparse and local nature of typical tight-binding Hamiltonians ensures that matrix-vector products, the computational core, scale as $\mathcal{O}(N)$ (Fan et al., 2013, 1705.01387).

2. Representations of Conductivity: Kubo-Greenwood and Einstein Formulations

The formalism for quantum transport in the linear scaling approach most often relies on the Kubo-Greenwood formula for DC conductivity:

$$\sigma(E) = \frac{\pi \hbar e^2}{\Omega} \operatorname{Tr}\!\left[ \delta(E-\hat{H})\, \hat{V}\, \delta(E-\hat{H})\, \hat{V} \right]$$

where $\hat{V}$ is the velocity operator and $\Omega$ is the system volume (Fan et al., 2018).

Two equivalent but complementary temporal representations are used:

• Green-Kubo (velocity autocorrelation):

$$\sigma^{\mathrm{GK}}(E,t) = e^2 \rho(E) \int_0^{t} C_{vv}(E, t')\,dt',$$

with $C_{vv}(E, t)$ the velocity autocorrelation function, and $\rho(E)$ the density of states. This approach is used to obtain running (time-dependent) conductivities in different transport regimes.

• Einstein (mean square displacement):

$$\sigma^{\mathrm{E1}}(E, t) = \frac{e^2 \rho(E)}{2} \frac{d}{dt} \Delta X^2(E, t), \qquad \sigma^{\mathrm{E2}}(E, t) = \frac{e^2 \rho(E)}{2t} \Delta X^2(E, t),$$

where all relevant quantities (MSD, VAC) are obtained by time-evolving random vectors, enabling efficient extraction of conductivity (Fan et al., 2013, 1705.01387).

Both approaches are algorithmically equivalent in the diffusive regime; discrepancies arise in deep localization due to the breakdown of underlying assumptions.

3. Numerical Implementation, GPU Acceleration, and Performance

Production implementations (e.g., GPUQT (1705.01387)) exploit CUDA-based sparse matrix multiplications, recursive Chebyshev and KPM expansions, and fully device-resident data structures. The density of states and related spectral functions are computed using the KPM:

$$\rho(E) \approx \frac{2}{\pi \Omega \Delta E \sqrt{1 - \tilde{E}^2}} \sum_{n=0}^{N_m-1} g_n (2 - \delta_{n0}) T_n(\tilde{E}) C_n^{\mathrm{DOS}}$$

where the KPM moments $C_n$ are computed iteratively.

Comparative benchmarks indicate:

• Linear scaling of computational time and memory with system size.
• Typical speedup factors of up to $16\times$ (double precision) on GPU vs. serial CPU for $>10^6$ orbitals in sparse tight-binding models (Fan et al., 2013).
• Applicability to two-dimensional (e.g. graphene) and one-dimensional (chains or nanoribbons) systems with disorder and defects, supporting systems with up to $N \sim 20$ million orbitals (1705.01387).

The method efficiently covers ballistic, diffusive, and localization regimes; in localization, extraction of the localization length involves exponential decay fits to conductance as a function of system size.

4. Approximation of Delta Functions and Moments: KPM vs. Fourier Transform

The efficiency of the linear-scaling approach depends on accurate and computationally tractable approximations to the Dirac delta operator required for spectral resolution. Two systematic expansions are contrasted:

Method	Expansion Basis	Efficiency (moments per energy point)	Sensitivity to Parameters
KPM	Chebyshev polynomials	Highest (1 expansion per moment)	Stable to time/energy discretization
FTM	Fourier series	Lower (multiple propagations needed)	Time step is bandwidth limited

The KPM, with Jackson damping, is on average 6.4 times more efficient than the FTM for equivalent energy resolution (Fan et al., 2013). This efficiency is vital when spectral functions are needed at many energies or for large-scale parametric studies.

5. Transport Regimes: Ballistic, Diffusive, and Localization

The method enables clear discrimination between transport regimes:

• Ballistic: The MSD grows quadratically ($\Delta X^2 \propto t^2$), and the conductance exhibits quantized plateaus for nanoribbons, matching NEGF predictions. Overshoots near band edges arise due to singularities in $\rho(E)$ and velocity squared at band extrema.
• Diffusive: MSD crosses over to linear growth; Green-Kubo and Einstein conductivity definitions become numerically equivalent. The semi-classical conductivity is extracted from the conductivity plateau.
• Localization: The MSD saturates at long times; the conductance decays exponentially with system size. Caution is required as length definitions become ambiguous and different conductivity formulas can yield diverging results.

In deeply localized regimes, empirical Padé fits to MSD as a function of time are employed to stabilize extraction of relevant quantities.

6. Applications to Graphene and Disordered Systems

The methodology has been validated in:

• Pristine and vacancy-disordered graphene: prediction of minimum conductivity plateaus, localization transitions, and sub-diffusive behavior at Dirac points.
• Armchair and zigzag nanoribbons: quantized conductance, edge state transport, and disorder-induced localization.
• Large-scale 2D models: demonstration of minimum conductivity and effects of grain boundaries, supported by empirical formulas for sheet resistance vs. grain size (Fan et al., 2018).

Extensions include spin and valley Hall conductivities, polaron/phonon-induced transport, and integration with first-principles models for material-specific Hamiltonians.

7. Methodological Extensions and Connections

• Spin and Topological Transport: By initializing stochastic spin- or valley-polarized states, the same linear scaling infrastructure yields spin relaxation times and topological contributions to Hall conductivities (via Kubo-Streda and Kubo-Bastin formulas).
• Wavefunction Matching and NEGF: The order-$N$ principle remains essential in wavefunction matching (WFM) and NEGF approaches, where sparse solvers for the central region and reduction to open channel bases in leads yield $\mathcal{O}(N_SN_P)$ scaling (Santos et al., 2018).
• Quantum-Classical and Quantum Algorithmic Approaches: Linear-scaling principles are now being adopted and extended in variational quantum algorithms, Hamiltonian simulation for PDEs, and hybrid machine learning frameworks for electronic structure and transport (Bengoechea et al., 22 Nov 2024, He et al., 25 Jul 2025, Zylberman et al., 21 Aug 2025).

Summary Table: Core Features

Feature	Approach	Scaling	Key Algorithmic Step
Trace evaluation	Stochastic/random phase	$\mathcal{O}(N)$	Matrix-vector multiplication
Time/energy evolution	Chebyshev/KPM or Fourier expansion	$\mathcal{O}(N)$	Recursive polynomial expansion
Disorder/defects	Directly included in Hamiltonian
Transport regimes	Ballistic/diffusive/localization		MSD/VAC analysis
Parallelization	GPU/CUDA, sparse linear solvers		Shared-memory parallelism

Conclusion

The linear scaling quantum transport methodology employs stochastic trace estimation, Chebyshev polynomial expansions, and sparse matrix multiplication to compute time- and energy-dependent transport properties with computational cost linear in system size. These approaches enable the calculation of electronic, spin, and topological transport in realistic models comprising millions of orbitals, in both ordered and strongly disordered media. The methods are validated against NEGF benchmarks, faithfully capture the physics of ballistic to localized transitions, and serve as a foundation for extensions to correlated, topological, and ab initio-informed quantum transport simulations. Their efficiency, extensibility, and accuracy have established them as essential computational tools for contemporary investigations of mesoscopic electronic transport.