with θ(q) being the Heaviside step function. Note that we take the left well i.e., q < 0 to be the reactant region. As in our previous work , the corresponding classical isomerization probabilities are computed by choosing an ensemble of initial conditions N tot sampled from the classical density ρ cl (q, p, x c , p c , 0) corresponding to the initial quantum polaritonic wavepacket. For this study, we take N tot = 50000 and time evolve each initial phase space point by integrating the Hamiltonian equations of motion. Classically, starting from the reactant well, an isomerization event corresponds to the trajectory crossing the saddle (q ‡ , x ‡ c ) = (0, 0) towards the product well. The classical isomerization probability is then calculated by averaging the isomerization events over the initial classical distribution.