The numerical results obtained from L-MFE were benchmarked against the original Lindblad dynamics in Eq. 21, which are referred to as "exact" results (of performing Lindblad dynamics) in this paper. These numerical results are obtained by using the QuTiP library with the mesolve function, where model Hamiltonian, Lindblad jump operator, and initial wavefunction are entered as arguments to the function. Similar numerical simulations has been recently performed to investigate molecular cavity QED processes as well. For Models 1-8, ĤC + ĤQC = 0 (no nuclear DOF), and the only inputs in the simulation are the matrix elements of ĤQ and the decay rate Γ. For the photo-isomerization coupled to the cavity model, ĤS is described by Eq. 82. The matrix elements of the Hamiltonian in Eq. 82 as well as the jump operator LS in Eq. 83 are evaluated using the basis {|α, n⟩ ⊗ |χ ν ⟩}, where |χ ν ⟩ is the discrete variable representation (DVR) basis for the nuclear DOF R. For the DVR basis |χ ν ⟩, a total of 175 grid points are used in the range of R ∈ [-1.25 a.u., 1.25 a.u.].