where c(t) is the vector of the quantum coefficients, T describes the coefficients update due to the jump operator using Eqs. 69a-69c, and e -i ℏ ( ĤQ+ ĤQC)dt describes the unitary evolution of the coefficients due the MFE dynamics described in Eq. 15. Note that ĤQC (R) depends on the nuclear DOFs R, where R(t) evolves based on the force in Eq. 16. Here, we use a symmetrical Trotter decomposition in Eq. 72 to reduce error due to a finite time step dt, and we put the Ehrenfest propagation in the middle because it is computational expensive (which scales as ∼ K 2 ) compared to the Lindblad decay part (denoted by T ) which only requires the update of two coefficients c 1,ξ (t) and c 0,ξ (t).