where a short-time evolution by d t can be approximated to first-order by the symmetric, factorized product of individual exponentiated Liouville operators. This leads to a time-reversible and symplectic propagation for the simulated cavity loss, at a similar level to a velocity-Verlet-like scheme. The time-evolution can then be obtained as a step-wise propagation using the following equations of motion,