  While several numerical techniques are available for predicting the dynamics
of non-Markovian open quantum systems, most struggle with simulations for very
long memory and propagation times, e.g., due to superlinear scaling with the
number of time steps $n$. Here, we introduce a numerically exact algorithm to
calculate process tensors -- compact representations of environmental
influences -- which provides a scaling advantage over previous algorithms by
leveraging self-similarity of the tensor networks that represent Gaussian
environments. Based on a divide-and-conquer strategy, our approach requires
only $\mathcal{O}(n\log n)$ singular value decompositions for environments with
infinite memory. Where the memory can be truncated after $n_c$ time steps, a
scaling $\mathcal{O}(n_c\log n_c)$ is found, which is independent of $n$. This
improved scaling is enabled by identifying process tensors with repeatable
blocks. To demonstrate the power and utility of our approach we provide three
examples. (1) We calculate the fluorescence spectra of a quantum dot under both
strong driving and strong dot-phonon couplings, a task requiring simulations
over millions of time steps, which we are able to perform in minutes. (2) We
efficiently find process tensors describing superradiance of multiple emitters.
(3) We explore the limits of our algorithm by considering coherence decay with
a very strongly coupled environment. The algorithm we present here not only
significantly extends the scope of numerically exact techniques to open quantum
systems with long memory times, but also has fundamental implications for
simulation complexity.
