To reduce this unfavorable scaling related to density matrix propagation in the Liouville space, one can equivalently describes dynamics as a stochastic wavefunction in the open system's Hilbert space that exactly reproduces the reduced density matrix of the open system, thus reducing the cost of solving the same dynamics by utilizing an ensemble of trajectories in the Hilbert space. This is referred to as the unravelling of the master equation. It is well known that the Lindblad mater equation can be equivalently expressed as the stochastic Schrödinger equation (SSE) with jump processes that captures the projective action of the environment on the system, which randomly collapses the system wavefunction into a pure system state. Through the trajectory average of a large number of realizations of the jump trajectories, the SSE generates identical reduced density matrix dynamics as the Lindblad master equation. However, due to the large fluctuations in population of the system, the SSE encounters numerical instability when combined with mixed quantum-classical approaches. The Ehrenfest+R method was recently developed to simulate the electronic quantum subsystem coupled to the classical electromagnetic field in order to accurately describe spontaneous emission processes. It effectively captures Lindblad dynamics with a deterministic change of the magnitude of the quantum coefficients, and stochastic changes of the phases. The Ehrenfest+R method, however, cannot exactly reproduce the Lindblad master equation due to a specific choice of the off-diagonal decay rate.