(26) From Eq. 26, it can be seen that the jump operator L = |0⟩⟨1| causes exponential decay of the excited state to the ground state with a rate of Γ and exponential decay of the excited-ground coherence with a rate of Γ/2. This dynamical picture approximately corresponds to many physical phenomena such as spontaneous emission and photonic cavity loss when a single decay channel dominates. Expressing the MFE approach in the density matrix form (Eq. 17) allows for Lindblad dynamics to be implemented in a straightforward manner. This involves performing MFE dynamics in the Liouville space where each element |ψ a ⟩⟨ψ b | is now treated as a basis vector. Eq. 19 and Eq. 21 describe the nuclear force and equation of motion of the density matrix elements, respectively. However, this density matrix formalism of MFE requires K 2 density matrix elements to be dynamically updated for a dimension K of the original quantum subsystem Hilbert space H Q , with a formal scaling in the range of N • K 3 to N • K 4 in terms of computational cost, where N is the number of trajectories used to perform the MFE part of the Liouvillian in Eq. 21. The coefficient formalism MFE (Eq. 14), on the other hand, only requires K coefficients to be dynamically updated, with a formal cost of N • K 2 . Due to the numerically favorable scaling of a wavefunction based approach, together with the consideration that the trajectory average is required anyway in performing the MQC dynamics governed by Eq. 15, we choose to incorporate Lindblad dynamics into the regular coefficient-based approach of MFE. This requires a formulation of Lindblad dynamics in a coefficient-based formalism, which is derived in the next section.