A. Computational scaling with N For an N -molecule polaritonic system, including the 2 nd excitation manifold leads to a K ∝ O(N 2 ) number of states, and O(N 4 ) numbers of the reduced density matrix elements. The operation of acting ĤQ on ρQ would in principle requires O(N 6 ) operations. The total cost of 2DES simulation for the N -molecule polaritonic system, for M trajectories, will be O(N 6 T 3 M ). This scaling is for both memory and computational cost (run time). However, using the strategies developed in the previous sections, one never needs to generate any O(N 4 ) density matrix, and all of our calculations can be performed with vectors of the size of O(N 2 ). This is the case for any semi-classical or mixed quantum-classical trajectorybased approach, as we extensively discussed in Paper I. For a time step ∆t (where the nuclear configuration is fixed), the propagation of the reduced density matrix can be expressed as applying the unitary propagator Û (∆t) In Fig. , we present the computational scaling cost of different schemes for the operation of Hamiltonian ( ĤQ ) of the reduced density matrix of the system (ρ) or the system state vector (|Ψ⟩). The solid red curve is the computational cost of performing a ĤQ ρ operation. Each of the matrices ĤQ and ρ scales as O(N 4 ) and thus their multiplication here scales as O(N 6 ), as can be seen in the red curve. The solid yellow curve demonstrates the computational cost of performing a ĤQ |Ψ⟩ operation which scales as O(N 4 ). The solid cyan curve is the cost of simulating the vectorized sparse ĤQ |Ψ⟩ operation, and scales as O(N 2 ). The dotted curves represent the computational time taken by these operations, when performed on a GPU . Performing the same calculations on the GPU shows a constant scaling, O(N 0 ), with increasing N . Although the matrix-vector operation, ĤQ ρ, seems very efficient on GPU, these operations occupy a lot of memory, and thus we are restricted to only a few ĤQ ρ operations. ĤQ |Ψ⟩ reduces the memory requirement, but we are still highly restricted in the number of ĤQ |Ψ⟩ operations we can perform. To address this challenge, we move to vectorized ĤQ |Ψ⟩ operation, where the memory requirement is significantly reduced to a more tractable regime, and we can simultaneously perform a lot of Hadamard products with a reduced computational cost.