A staging Monte Carlo algorithm for sampling off-diagonal density matrix elements via open-chain path integrals

We introduce a computationally simple algorithm for sampling open-chain distributions within the framework of imaginary-time Feynman path integration. The present method is based on the staging algorithm introduced by Pollock and Ceperley [Phys. Rev. B 30, 2555 (1984)] originally developed for computing position-dependent observables. Here, we sample off-diagonal elements of the density matrix, formulated as a distribution describing a linear polymer-like chain of beads, each connected via nearest-neighbor springs to calculate momentum-dependent quantities. This is achieved using a Monte Carlo scheme that ensures efficient and unbiased sampling of all beads along the chain from the free-particle distribution via a staging transformation; we refer to this approach as staging open path integral Monte Carlo (OPIMC). The proposed algorithm is straightforward to implement, as it only involves sampling Gaussian distributions through a transformation defined by a set of recursion relations, followed by a standard Metropolis acceptance/rejection step. The staging OPIMC method accurately reproduces end-to-end and momentum distributions for quantum systems ranging from coupled harmonic oscillators to liquid water.