Quadratic scaling bosonic path integral molecular dynamics

Bosonic exchange symmetry leads to fascinating quantum phenomena, from exciton condensation in quantum materials to the superfluidity of liquid ⁴He. Unfortunately, path integral molecular dynamics (PIMD) simulations of bosons are computationally prohibitive beyond ~100 particles, due to a cubic scaling with the system size. We present an algorithm that reduces the complexity from cubic to quadratic, allowing the first simulations of thousands of bosons using PIMD. Our method is orders of magnitude faster, with a speedup that scales linearly with the number of particles and the number of imaginary time slices (beads). Simulations that would have otherwise taken decades can now be done in days. In practice, the new algorithm eliminates most of the added computational cost of including bosonic exchange effects, making them almost as accessible as PIMD simulations of distinguishable particles. This work is an important step towards computing transport properties of bosonic condensed phases using approximations based on path integral molecular dynamics, such as ring polymer molecular dynamics and centroid molecular dynamics.