where Tr ex and Tr ph represent the trace over the excitonic and photonic DOFs, separately. The above procedure can either be terminated after acting μ(t 2 ) to get the 1 st order non-linear spectra (e.g. transient absorption spectra) or can also be continued with further perturbations to get the n th order spectra. The advantage of the PLDM algorithm is its scalability, where (1) each trajectory can be propagated independently, and (2) for each trajectory, any particular branch (t 1 , t 2 , t 3 ) can be propagated independently of any other branch (t 1 ′ , t 2 ′ , t 3 ′ ). Thus, the algorithm is highly parallelizable, and the scaling of the non-linear spectroscopy computation can be efficiently controlled in high-performance computing facilities.