method for this model). Fig. ,c present the population and coherence dynamics, respectively, obtained from Ehrenfest+R (dots), using the same algorithm as described in Ref. 11, where a concise summary of this approach can be found in Appendix A. The +R decay dynamics are implemented as suggested as the original reference 10 (through Eq. A27a and Eq. A27b) using the original γ R (Eq. A22b). The parameters of Model 8 were carefully chosen such that the excited state energy E ′ = 1 a.u. is significantly larger than the decay rate Γ = 0.05 a.u., such that the Ehrenfest+R method should work correctly in the regime E ′ ≫ Γ, while still showing interesting dynamics on the timescale of the population decay (t = 0 ∼ 100 a.u.). In Fig. , the Ehrenfest+R dynamics are qualitatively similar to the exact Lindblad dynamics, but there are some errors in the fluctuations of the state |1⟩ decay as well as in the magnitudes of the oscillations of states |0⟩ and |2⟩. The cause of these deviations can be seen in Fig. , where the magnitudes of the coherences obtained from Ehrenfest+R are larger than those of Lindblad dynamics, causing incorrect population transfer between diagonal populations. This is because the coherence decay rate γ R (Eq. A22b) in the Ehrenfest+R algorithm is smaller than the correct decay rate γ ′ R (Eq. B11), causing artificially large coherences and thus a larger magnitude of population oscillation.