Modular path integral for discrete systems with non-diagonal couplings

The modular decomposition of the path integral, which leads to linear scaling with the system length, is extended to Hamiltonians with intermonomer couplings that are not diagonalizable in any single-particle basis. An optimal factorization of the time evolution operator is identified, which minimizes the number of path integral variables while ensuring high accuracy and preservation of detailed balance. The modular path integral decomposition is described, along with a highly efficient tensor factorization of the path linking process. The algorithm is illustrated with applications to a model of coupled spins and a Frenkel exciton chain.