## 2 Density Matrix Simulations  
A non-perturbative, time-dependent density matrix approach is used to simulate the interaction of a femtosecond laser pulse with a four level system consisting of 2 electronic states--ground and excited state, each containing two vibrational levels (Figure 1). All simulations were carried out under the slowly-varying envelope (SVEA) and rotating-wave (RWA) approximations, by using a 4-5 order Runge-Kutta algorithm (see also [65]). The laser temporal pulse shape is a direct input to the code, and the system allows the laser pulse to interact with all electronic transitions. The model can be written as follows:  
\[ \begin{array}{l}\begin{array}{l}\bullet\\ \rho_{22}\end{array}=-\gamma_{22}\rho_{22}+\gamma_{32}\rho_{33}+\gamma_{42} \rho_{44}+i\left\{\mu_{23}(\rho_{32}\Omega^{*}-\rho_{23}\Omega)\,+\,\mu_{24}( \rho_{42}\Omega^{*}-\rho_{24}\Omega)\right\},\\ \begin{array}{l}\bullet\\ \rho_{33}\end{array}=-\gamma_{33}\rho_{33}+\gamma_{43}\rho_{44}+i\left\{\mu_ {13}(\rho_{13}\Omega-\rho_{31}\Omega^{*})\,+\,\mu_{23}(\rho_{23}\Omega-\rho_{ 32}\Omega^{*})\right\},\\ \begin{array}{l}\bullet\\ \rho_{44}\end{array}=-\gamma_{44}\rho_{44}+i\left\{\mu_{14}(\rho_{14}\Omega- \rho_{41}\Omega^{*})\,+\,\mu_{24}(\rho_{24}\Omega-\rho_{42}\Omega^{*})\right\},\\ \begin{array}{l}\bullet\\ \rho_{11}\end{array}=-\begin{array}{l}\bullet\\ \rho_{22}\end{array}-\begin{array}{l}\bullet\\ \rho_{33}\end{array}-\begin{array}{l}\bullet\\ \rho_{44}\end{array}\), \\ \begin{array}{l}\bullet\\ \rho_{12}\end{array}=-\begin{array}{l}\overline{\Gamma}_{12}\rho_{12}+i \omega_{21}\rho_{12}+i\left\{\mu_{13}\rho_{32}\Omega^{*}-\mu_{23}\rho_{13} \Omega+\mu_{14}\rho_{42}\Omega^{*}-\,\mu_{24}\rho_{14}\Omega\right\},\\ \begin{array}{l}\bullet\\ \rho_{13}\end{array}=-\overline{\Gamma}_{13}\rho_{13}+i(\omega_{31}-\omega_{ 1})\rho_{13}+i\Omega^{*}\left\{\mu_{13}(\rho_{33}-\rho_{11})\,+\,\mu_{14}\rho _{43}-\,\mu_{23}\rho_{12}\right\},\\ \begin{array}{l}\bullet\\ \rho_{14}\end{array}=-\overline{\Gamma}_{14}\rho_{14}+i(\omega_{41}-\omega_{ 1})\rho_{14}+i\Omega^{*}\left\{\mu_{14}(\rho_{44}-\rho_{11})\,+\,\mu_{13}\rho _{34}-\,\mu_{24}\rho_{12}\right\},\\ \begin{array}{l}\bullet\\ \rho_{23}\end{array}=-\overline{\Gamma}_{23}\rho_{23}+i(\omega_{32}-\omega_{ 1})\rho_{23}+i\Omega^{*}\left\{\mu_{23}(\rho_{33}-\rho_{22})\,+\,\mu_{24}\rho _{43}-\,\mu_{13}\rho_{21}\right\},\\ \begin{array}{l}\bullet\\ \rho_{24}\end{array}=-\overline{\Gamma}_{24}\rho_{24}+i(\omega_{42}-\omega_{ 1})\rho_{24}+i\Omega^{*}\left\{\mu_{24}(\rho_{44}-\rho_{22})\,+\,\mu_{23}\rho _{34}-\,\mu_{14}\rho_{21}\right\},\\ \begin{array}{l}\bullet\\ \rho_{34}\end{array}=-\overline{\Gamma}_{34}\rho_{34}+i\omega_{43}\rho_{34}+i \left\{\mu_{13}\rho_{14}\Omega-\mu_{14}\rho_{31}\Omega^{*}+\mu_{23}\rho_{24} \Omega-\,\mu_{24}\rho_{32}\Omega^{*}\right\},\end{array}\end{array} \end{array} \tag{1}\]  
where \(\rho_{ii}\) are the populations of the levels, the off-diagonal elements \(\rho_{ij}\) represent the coherences between levels, and \(\mu_{ij}\) are the dipole coupling coefficients of the electronic transitions. \(\omega_{ij}=|E_{j}-E_{i}|/\hbar\) are the angular frequencies of the transitions \(i-j\), with \(E_{i,j}\) being the energies of the quantum states \(i\) and \(j\) (see Figure 1). \(\Omega(t)=E(t)\text{exp}[i\varphi(t)]/2\hbar\) represents the input slowly-varying complex amplitude of the electric field. \(E(t)\) is the slowly varying real field envelope, and \(\varphi(t)\) is the slowly varying phase, which for transform limited pulses is chosen constant. The decay rates of the transitions \(i-j\) are given by \(\gamma_{ij}\), and the relaxation times of the levels satisfy \(\gamma_{22}\)= \(\gamma_{21}\), \(\gamma_{33}\)= \(\gamma_{31}\)+\(\gamma_{32}\) and \(\gamma_{44}\)= \(\gamma_{41}\)+\(\gamma_{42}\)+\(\gamma_{43}\).