Markov Inequality as a Tool for Linear-Scaling Estimation of Local Observables

Technical Summary: Markov Inequality as a Tool for Linear-Scaling Estimation of Local Observables

Problem Statement
Revealing spatially resolved quantum observables—such as the Local Density of States (LDoS), local currents, and topological markers—in large disordered lattice systems remains a significant computational challenge. While linear-scaling algorithms based on Chebyshev expansions efficiently estimate volume-averaged spectral functions (operator traces), they struggle when applied to full-sample spatial maps. The generation of these maps typically incurs a computational cost scaling as $O(N^2)$ (where $N$ is the system size) due to the lack of self-averaging in local quantities. In systems with strong spatial fluctuations, such as those exhibiting Anderson localization or multifractal scaling near metal-to-insulator transitions, standard stochastic sampling methods fail to converge uniformly. Specifically, the relative sampling error in naive estimators scales inversely with the local wavefunction amplitude ($|\psi_\alpha(r)|^{-1}$), leading to non-uniform convergence and unphysical results (e.g., negative LDoS values) in regions where the observable is small.

Methodology
The authors introduce a linear-scaling stochastic method designed to compute real-space maps of any positive local spectral operator in tight-binding models with controlled accuracy. The core innovation lies in the use of positive-definite stochastic estimators combined with the Markov inequality.

1. Positive-Definite Estimators: Unlike conventional unbiased estimators that can yield sign-changing values, the proposed method constructs estimators that are strictly non-negative. For a positive-definite operator $X$, a random field $\phi_r = \langle \xi | \circ \sqrt{X} | \xi \rangle$ is introduced. The target diagonal element $X_r$ is then estimated as $|\phi_r|^2$ (for complex vectors) or $\text{Var}(\phi_r)$ (for real vectors satisfying specific moment conditions).
2. Markov Inequality Bound: By ensuring the estimator is non-negative, the Markov inequality ($P(\eta_r \ge C) \le \eta_r C^{-1}$) is employed to rigorously bound the sampling error at each site relative to the local mean. This ensures that the relative fluctuations remain constant across the sample, regardless of the magnitude of the local observable, thereby guaranteeing homogeneous convergence.
3. Extension to Off-Diagonal Observables: The method is extended to non-diagonal elements (e.g., local currents) through local unitary transformations. By applying a local rotation to pairs of sites $\{r, r'\}$, the off-diagonal matrix elements are encoded into the diagonal elements of a rotated basis. This allows the positive-definite estimation technique to be applied to off-diagonal quantities without inducing additional sampling correlations, as the rotations act on disjoint site pairs.

Key Results
The effectiveness of the method is demonstrated through benchmark calculations on the disordered 2D $\pi$-flux model on a square lattice:

• Single Vacancy Case: In a $4096 \times 4096$ supercell with a single central vacancy, the method successfully computed the zero-energy LDoS. It faithfully captured the expected algebraic decay ($\rho_r(0) \propto r^{-2}$) over several orders of magnitude. In contrast, the conventional naive estimator failed in low-amplitude regions, producing non-physical negative values and large relative errors.
• Multiple Vacancies: For a system with 0.5% vacancy concentration on a $2048 \times 2048$ supercell, the method generated converged LDoS maps spanning several orders of magnitude. The 95th percentile of local relative fluctuations was bounded at 3.5%, demonstrating controlled accuracy even in the presence of critical zero-energy states.
• Error Analysis: Comparisons with exact LDoS calculations on a $512 \times 512$ sample showed that the positive-definite estimator reduced the mode of the local relative error by approximately two orders of magnitude compared to the naive estimator. The error distribution scaled as $R^{-1/2}$ with the number of random vectors $R$. Computationally, a stochastic estimation took $\approx 4$ CPU hours, whereas the exact calculation required $\approx 171$ CPU days.
• Local Currents: The method was applied to compute steady-state local current fields in the presence of long-range potential scatterers. The resulting maps revealed intricate vortex–anti-vortex pairs and complex flow patterns, resolving structures associated with the localization transition in graphene (Kosterlitz–Thouless universality class) that were previously difficult to capture with high spatial resolution.

Significance and Claims
The paper claims that this approach enables the simulation of disorder-driven mesoscopic phenomena in realistically large lattices by reducing the computational complexity of generating full-sample maps from quadratic to linear scaling. The rigorous bounding of site-wise stochastic errors via the Markov inequality overcomes the fundamental limitation of lacking self-averaging in local quantities.

The authors assert that this method:

• Enables accurate estimation of local observables even under strong spatial fluctuations.
• Facilitates the study of multifractality and disorder-driven transitions in large systems.
• Accelerates real-space self-consistent mean-field calculations by providing efficient, controlled-accuracy local inputs.
• Provides a general framework applicable to any positive local spectral operator, extendable to off-diagonal observables via unitary transformations.

The work is presented as a tool to bridge the gap between experimental spatial resolution (e.g., STM, SQUID) and theoretical simulation capabilities in condensed-matter physics.