  The solution to Maxwell-Bloch systems using an integral-equation-based
framework has proven effective at capturing collective features of laser-driven
and radiation-coupled quantum dots, such as light localization and
modifications of Rabi oscillations. Importantly, it enables observation of the
dynamics of each quantum dot in large ensembles in a rigorous,
error-controlled, and self-consistent way without resorting to spatial
averaging. Indeed, this approach has demonstrated convergence in ensembles
containing up to $10^4$ interacting quantum dots. Scaling beyond $10^4$ quantum
dots tests the limit of computational horsepower, however, due to the
$\mathcal{O}(N_t N_s^2)$ scaling (where $N_t$ and $N_s$ denote the number of
temporal and spatial degrees of freedom). In this work, we present an algorithm
that reduces the cost of analysis to $\mathcal{O}(N_t N_s \log^2 N_s)$. While
the foundations of this approach rely on well-known
particle-particle/particle-mesh and adaptive integral methods, we add
refinements specific to transient systems and systems with multiple spatial and
temporal derivatives. Accordingly, we offer numerical results that validate the
accuracy, effectiveness and utility of this approach in analyzing the dynamics
of large ensembles of quantum dots.
