Computational Details. The L-MFE method is implemented by using a symmetrical Trotter decomposition in Eq. 72 to reduce time step error. The unitary dynamics e -i ℏ ( ĤQ+ ĤQC)dt are propagated using the 4th order Runge Kutta (RK4) algorithm, while the Lindblad decay dynamics T are propagated using the coefficient modifications described in Eqs. 69a to 71. When classical nuclear DOFs are present, the velocity Verlet algorithm is used to propagate the nuclear DOFs, with mean-field force given by Eq. 16, alongside the quantum subsystem. For all model calculations, the electronic time step used is dt E = 0.05 a.u. For the photo-isomerization system that contains nuclear DOF, the nuclear time step was dt N = 6 a.u. and the time step used in the Lindblad decay T was dt N /4 (thus the unitary propagation was performed 60 times in between Lindblad decay propagations). A total of 24, 000 trajectories were used to ensure fully converged results, although using only 1000 trajectories already provides a mostly converged result (as shown in Fig. ).