Systems in contact with an environment provide a ubiquitous challenge in quantum dynamics. Many fascinating phenomena can arise if the coupling is strong, leading to non-Markovian dynamics of the system, or collective, where the environment can become highly excited. We introduce a significant improvement of the hierarchy of pure states (HOPS) approach, which is able to efficiently deal with such highly excited, non-Markovian environments in a nearly unitary way. As our method relies on quantum trajectories, we can
obtain dynamics efficiently, also for large system sizes, by (1) avoiding the quadratic scaling of a density matrix and (2) exploiting the localization properties of the trajectories with an adaptive basis. We provide the
derivation of the nuHOPS (nearly unitary hierarchy of pure states) method, compare it to the original HOPS, and discuss numerical subtleties based on an illustrative dephasing model. Finally, we show its true power using
the Dicke model as the paradigmatic example of many emitters decaying superradiantly inside a cavity. We reach numerically exact solutions for up to 1000 emitters. We apply our method to study emerging higher-order
correlations in the emitter system or the cavity mode environment and their scaling with the number of emitters.