As discussed above in Fig. , the explicit inclusion of many-molecules yields unfavorable conditions for the direct diagonalization (Sec. 3.1) of the Hamiltonian. In the case of self-consistent approaches (Sec. 3.2), one encounters a similar bottleneck where the inclusion of many sets of one-molecule single-particle orbitals N sp = N occ × N vir leads to a many-molecule ground state density matrix of dim[ ρ ] ∝ N sp × N mol . The off-diagonal blocks of which will include the long-range Coulomb, short-range exchange, and intermolecular DSE interactions. In recent work, 181 it was argued that the off-diagonal blocks of the total density matrix can be decoupled in the cavity Born-Oppenheimer (BO) approximation in which the cavity photon is treated as a parameter (similar to the nuclear BO approximation). Thus, the diagonal blocks of the density matrix and treated as effectively independent, reducing the memory cost of the selfconsistent treatment back to dim[ ρ ] ∝ N sp but with N mol such matrices. However, without approximations, such as cavity BO, treating many molecules coupled to many cavity modes is still beyond the reach of direct simulation. Relying on efficient semi-empirical or tight-binding approaches 182;183;184;185 may provide the most direct route toward the direct simulation of many molecule QED, including both ground and excited states and has already been used to simulate plasmonic systems. 186;187;188;189