only indirectly via virtual processes, and hence, they do not need to be treated explicitly. The corresponding electron-boson couplings are assumed to be integrated out and enter in the mean-field treatment of superconductivity via $W_0$.

For our numerical calculations we have to fix the parameters. Motivated by the numbers for Pb [@pooleBook07] we fix the parameters as follows: gap in the initial state $\Delta (t_{\textrm{i}} ) =1.35$ meV, Fermi energy $E_{\textrm{F}}=9479$ meV, energy cutoff $\hbar\omega_{\textrm{c}}=8.3$ meV, and effective electron mass $m=1.9m_0$, with $m_0$ being the free electron mass. The operator $b^{\dagger}_{\textbf{p}}$ ($b^{\phantom{\dagger}}_{\textbf{p}}$) in Eq.  creates (annihilates) phonons with wave vector $\textbf{p}$ and frequency $\omega_{\text{ph}}$, where $\omega_{\text{ph}}$ is assumed to be constant for simplicity. Similarly, the coupling between electrons and phonons is taken to be of Holstein form \[Eq. \] with a momentum-independent interaction constant $g_{\text{ph}}$.

In the following, we study the nonequilibrium response of Hamiltonian  to a short intense pump pulse which injects a nonthermal distribution of Bogoliubov quasiparticles into the system. The considered pump pulse is of Gaussian shape with photon frequency $\omega_{\textrm{p}}$, photon wave vector ${\bf q}_{\textrm{p}}=|\textbf{q}_{\textrm{p}} | \hat{\textbf{e}}_x$, full width at half maximum (FWHM) $\tau_{\textrm{p}}$, and amplitude ${\bf A}_{\textrm{p}} =|\textbf{A}_{\textrm{p}} | \hat{\textbf{e}}_{y}$. Working within the Coulomb gauge, 