Gradient-based stochastic estimation of the density matrix

Fast estimation of the single-particle density matrix is key to many applications in quantum chemistry and condensed matter physics. The best numerical methods leverage the fact that the density matrix elements $f(H)_{ij}$ decay rapidly with distance $r_{ij}$ between orbitals. This decay is usually exponential. However, for the special case of metals at zero temperature, algebraic decay of the density matrix appears and poses a significant numerical challenge. We introduce a gradient-based probing method to estimate all local density matrix elements at a computational cost that scales linearly with system size. For zero-temperature metals the stochastic error scales like $S^{-(d+2)/2d}$, where $d$ is the dimension and $S$ is a prefactor to the computational cost. The convergence becomes exponential if the system is at finite temperature or is insulating.