Applying these two theoretical techniques, we extended the L-PLDM method 80 to the collective coupling regime. Using this approach, we simulate the linear spectra of the N -molecule polaritonic system described by the HTC Hamiltonian. The approach we present here will significantly reduce computational costs. To demonstrate the developed methods, we simulate the polaritonic absorption spectra of N = 500 molecules coupled to a cavity under either static energy disorder or dipole orientation disorder, under which the permutation symmetry of the HTC Hamiltonian is broken and a simple mean-field approach is not directly applicable. We observe similar trends of these linear spectra as observed in the previous work. As the disorder in energy σ ε increases, the linewidth increases, and the Rabi splitting also gradually increases due to the disorders in energies. When gradually increasing the orientation disorder, the Rabi splitting decreases, without further change in the line width. In future work, we will outline the theoretical details of stimulating 2DES spectra for the polariton system under the collective coupling regime (Paper II), with new computational techniques for focusing the al-gorithm of the PLDM simulation. and r ∈ [0, 1] is a uniform random number. The phase factor, when averaged over the ensemble of trajectories ensures the exponential population decay of the state |G 1 ⟩⟨G 1 | and a gain in population of the state |G 0 ⟩⟨G 0 | as required by Lindblad dynamics, as analytically proved in our previous work . The rest of the forward and backward coefficients, for a ̸ = G 1 and a ̸ = G 0 are not affected by the Lindblad decay, and the corresponding η a = ηa = 1. When averaging among a reasonable amount of trajectories (typically ∼ 10 2 trajectories, see Appendix B in Ref. 47), Eq. A5-A6 exactly reproduces the Lindblad EOM for reduced density matrix, as demonstrated numerically in our previous work.