Effects of non-linear coupling between the system and the bath vibrational modes on the system internal conversion dynamics are investigated using the Dirac–Frenkel variational approach with a newly defined sD₂ ansatz. It explicitly accounts for the entangled system electron-vibrational states, while the bath quantum harmonic oscillator states are expanded in a superposition of quantum coherent states. Using a non-adiabatically coupled three-level model, we show that efficient irreversible internal conversion due to quadratic vibrational-bath coupling occurs when the initially populated system vibrational levels are in resonance with the vibrational levels of a lower energy electronic state, also, a non-Gaussian bath wavepacket representation is required. The quadratic system-bath couplings result in broadened and asymmetrically squeezed bath quantum harmonic oscillator wavepackets in the coordinate-momentum phase space.