We can expand the equations of motion in the adiabatic-Fock basis (Sec. 3.1.1). We define the matrix elements of the polaritonic Hamiltonian including the first-order non-adiabatic couplings as H αn,β m = (H pl ) αn,β m -i ħ d αβ • Rδ nm , where we have used the fact that the Fock states have no dependence on the nuclear coordinates d αn,β m = d αβ δ nm and that R is the classical velocity of the nuclei. Further, the matrix elements of the Lindblad jump operator can be written as L αn,β m = √ mδ n,m -1 δ αβ , which has no affect on the electronic states. Fig. presents results of the computational method by combining (Fig. ) Ehrenfest and (Fig. ) surface hopping dynamics with a stochastic extension to the previously outlined Lindblad decay. For details pertaining to the stochastic analogue to the aforementioned Lindblad loss approach, see Ref. 176. This approach stochastically maps the dynamics in Liouville space to dynamics in the Hilbert space. We refer to this method as the L-MFE approach. The convergence of this approach 176 with respect to the number of trajectories is similar to any mixed quantum-classical method, thus the addition of the stochastic cavity loss does not increase the computational cost. This is a general scheme that can be easily incorporated into any trajectory-based dynamics approaches, such as Ehrenfest, FSSH, or mapping approaches, or additionally to any trajectory-guided wavepacket-based methods, such as ab initio multiple spawning or ab initio multiple cloning.