Competing pathways in polaritonic chemistry and general harmonic environments

Recent experiments suggest that strong light-matter coupling between molecules and cavity modes can significantly modify chemical reaction rates and yields. Despite multiple proposed mechanisms, a clear microscopic understanding of how vibrational strong coupling alters reaction kinetics remains incomplete. Here, we employ a general quantum-dynamical framework for analyzing competing reaction pathways in polaritonic chemistry by considering reaction coordinates with multiple relevant transition frequencies and general harmonic environments. We extend on a tensor-network formulation of the hierarchical equations of motion (HEOM) method to compute probability currents between vibrational states, leveraging a concise formulation for effective spectral densities and a controlled compression of Markovian bath modes into system-level jump operators. This allows us to decompose the resulting state-to-state probability currents, revealing how cavity-induced changes in individual transitions propagate through multistep reaction networks. We find that these local modifications can alter kinetic bottlenecks and selectively redistribute reaction flux between competing pathways, changing branching ratios under irreversible conditions without substantially affecting the equilibrium yield in the reversible regime. Finally, for coupled reaction-coordinate, spectator-mode, and cavity environments, we find structured, asymmetric rate profiles arising from interference between hybridized decay channels. Together, these results provide a mechanistic framework for understanding pathway-specific chemical reactivity under vibrational strong coupling.