Small matrix path integral in imaginary time

Thermal equilibrium properties are usually obtained from the imaginary-time path integral representation of the Boltzmann operator in combination with Monte Carlo integration methods. In some situations (identical fermions or frustrated Hamiltonians), the Boltzmann matrix leads to terms of alternating sign, which leads to a sign problem that severely impacts convergence. In this paper, we develop a robust and efficient quadrature-based method suitable for computing the Boltzmann matrix for discrete systems coupled to common or local harmonic baths. By expressing the discretized path integral with the influence functional in terms of a sum of matrix products, we develop a small matrix path integral (SMatPI) decomposition that allows iterative propagation in imaginary time while circumventing the storage of tensors employed in earlier work. The method is illustrated with several examples that involve two- and three-level systems coupled to common or local baths. We show that cyclic tight-binding Hamiltonians with positive coupling parameters give rise to Boltzmann matrix elements with alternating signs, presenting a severe sign problem to Monte Carlo approaches, while the SMatPI algorithm is stable and efficient.