where ϵ α is the diabatic polariton energy defined in Eq. 11, λ j = ⟨j| Ĥsb (R)|j⟩ is the pure diagonal contribution of the system-bath coupling (that parametrically depends on R), and c j (t) = ⟨j|Ψ(t)⟩. In principle, one needs to diagonalize Ĥs to obtain ϵ α and {⟨α|j⟩} at a cost of O(K 2 ) (see Ref. 62-64, due to the symmetry of Ĥs , instead of the usual cost of O(K 3 ) for diagonalizing an arbitrary K × K matrix). Nevertheless, because these are diabatic quantities, one can in principle compute them only once and store them.