  Time-resolved diffuse scattering experiments have gained increasing attention
due to their potential to reveal non-equilibrium dynamics of crystal lattice
vibrations with full momentum resolution. Although progress has been made in
interpreting experimental data on the basis of one-phonon scattering,
understanding the role of individual phonons can be sometimes hindered by
multi-phonon excitations. In Ref. [{\it arXiv:2103.10108}] we have introduced a
rigorous approach for the calculation of the all-phonon inelastic scattering
intensity of solids from first-principles. In the present work, we describe our
implementation in detail and show that multi-phonon interactions are captured
efficiently by exploiting translational and time-reversal symmetries of the
crystal. We demonstrate its predictive power by calculating the scattering
patterns of monolayer molybdenum disulfide (MoS$_2$), bulk MoS$_2$, and black
phosphorus (bP), and we obtain excellent agreement with our measurements of
thermal electron diffuse scattering. Remarkably, our results show that
multi-phonon excitations dominate in bP across multiple Brillouin zones, while
in MoS$_2$ they play a less pronounced role. We expand our analysis for each
system and examine the effect of individual atomic and interatomic vibrational
motion on the diffuse scattering signals. We further demonstrate the
high-throughput capability of our approach by reporting all-phonon scattering
maps of 2D MoSe2, WSe2, WS2, graphene, and CdI2, rationalizing in each case the
effect of multi-phonon processes. As a side point, we show that the special
displacement method reproduces the thermally distorted configuration that
generates precisely the all-phonon diffuse pattern.
