\chapter{Non-equilibrium interacting theory} \label{sec:thz-nca} The driven-dissipative lattice problem is straightforward: a pulse of electromagnetic radiation arrives, interacts with a heavy-fermion system and propagates away, carrying some information acquired in the interaction. The model is given by the auxiliary-particle light-matter Anderson lattice Hamiltonian (\cref{sec:alm-light}) \begin{equation} \label{eq:full-model} \hat H = - \sum_{\ev{i,j}} t^c_{ij}\hat c^\dagger_{i\sigma}\hat c_{j\sigma} % - \sum_{\ev{i,j}} t^f_{ij}\hat f^\dagger_{i\sigma} \hat b_{i} \hat b^\dagger_{j} \hat f_{j\sigma} + \sum_{i \sigma} \varepsilon_0^f \hat f^\dagger_{i\sigma}\hat f_{i\sigma} + \sbr{V_0 - i g_0 \del{\hat a - \hat a^\dagger}} \sum_{i \sigma} \del{\hat c^\dagger_{i \sigma} \hat b^\dagger_i \hat f_{i \sigma} + \hc} \end{equation} where the hopping of the $f$ electrons was set to zero and $\hat a$ is the annihilation operator of the photons of the pulse of radiation that drives the system (\cref{sec:travelling-pulse}). Furthermore, the conduction electrons $c$ are coupled to a fermionic heat bath (\cref{sec:fermionic-bath}) at the temperature of the cryostat. Upon tracing out the reservoir, an additional non-interacting effective action term appears \begin{equation} i S_{\mathrm{eff}}^{\mathrm{bath}} = \int_\gamma \dif z \dif z'\, \sum_{i \sigma} c^*_{i \sigma}(z) \sbr{\alpha^2 \Delta_{c}(z, z')} c_{i \sigma}(z') \tag{\ref{eq:fermionic-bath-action}}~, \end{equation} where the density of states of the bath is taken to be the same as the $c$ electrons'. The term $\alpha^2 \Delta_c(z, z')$ can be regarded as an additional hybridisation function of conduction electrons, containing thermal correlations from the interaction with the heat bath and serving as a dissipative/thermalising channel for the lattice system. The remaining (non-pulse) modes of the electromagnetic field are also treated as a reservoir (\cref{sec:photonic-bath}), for which the interacting action of the system gains the additional term upon tracing out the reservoir's degrees of freedom \begin{equation} i S_{\mathrm{eff}}^{\mathrm{bath}} = \resizebox{0.85\linewidth}{!}{$\int_\gamma \dif z\,\dif z' \sum_{i \sigma} \sbr{c_{i \sigma}^*(z) b_i^*(z) f_{i\sigma}(z) + \hc} i g_0^2 \eta \sbr{\Pi_a(z, z') + \Pi_a(z', z)} \sbr{c_{i \sigma}^*(z') b_i^*(z') f_{i\sigma}(z') + \hc}$} \tag{\ref{eq:photonic-bath-action}}~. \end{equation} In free space, the coupling with the environment is typically much stronger than with the pulse, as there are many more modes with which the system can couple. This is controlled by the dimensionless parameter $\eta \approx \frac{4\pi}{\Omega_\mathrm{p}} \del{1 - \frac{\Omega_\mathrm{p}}{4\pi}}$, which can be thought as a geometric factor that covers the solid angle of all modes except the pulse's~\cite{Silberfarb_2004}. After the excitation, the system will hence predominantly spontaneously emit into the modes of the environment~\cite{Albarelli_2022}. \begin{figure}[!htb] \centering \includegraphics[width=\figwidth]{diagrams/sketches/lattice} \caption{Driven-dissipative lattice problem. Figure adapted from~\cite{Scarlatella_2021}.} \label{fig:lattice} \end{figure} In its current form, it is hopeless to solve the problem as it is a lattice of interacting particles, for which not only is there no closed-form solution, brute-forcing perturbative expansions would unearth a diagrammatic hydra, requiring careful summations of non-local interacting terms. For this reason, two \textit{herculean} approximations are yet to be employed: dynamical mean-field theory, which will map the problem to an effective single-site problem and the non-crossing approximation, which will truncate the perturbative expansion of the single-site, or local, problem. \section{Dynamical mean-field theory} \label{sec:dmft} Dynamical mean-field theory (DMFT) heavily contributed to the understanding of many-body strongly-interacting lattice systems in thermal equilibrium~\cite{Georges_1996}, especially Mott physics. In its heart lies the observation~\cite{Metzner_1989} that the self-energy of itinerant systems becomes local in the limit of infinite dimensions. This, in turn, allows the mapping~\cite{Georges_1992} of a lattice problem into an Anderson-impurity problem of a single-site embedded in a renormalised conduction electron sea. Even though the renormalised sea neglects \textit{all} quantum spatial fluctuations\footnote{Refer to~\cite{Rohringer_2018} for a review on the inclusion of non-local effects.}, it retains the full information about quantum temporal fluctuations. Fortuitously, due to the DMFT construction being solely dependent on the real-space properties of $n$-point functions, the extension of DMFT to non-equilibrium settings~\cite{Schmidt_2002, Aoki2014} is straightforward. \subsection{The limit of infinite dimensions} In the limit of infinite dimensions $d\to\infty$, the hopping amplitude $t$ in Hamiltonians such as the Hubbard or Anderson lattice model must be re-scaled such that its competition with local effects such as the Coulomb interaction remains non-trivial, as otherwise, the kinetic energy density would be infinite. The re-scaling \begin{equation} \label{eq:dmft-rescale} t \to \frac{t}{\sqrt{Z}}~, \end{equation} where $Z$ is the lattice connectivity/coordination-number ($Z=2d$ in a $d$-dimensional hyper-cubic lattice) has its most important consequence on the scaling of the non-interacting lattice Green function \begin{equation} \label{eq:dmft-gf-scaling} G_{0,ij} \stackrel{d\to\infty}{=}\mathcal{O}(\sqrt{1/Z}^{d(i,j)})~, \end{equation} where $d(i,j)$ denotes the distance between sites $i$ and $j$. The scaling, however, does not imply that the system becomes localised -- from a probabilistic perspective, it simply means that for a lattice with infinite connectivity, the probability density of going from the site $i$ to $j$ is infinitesimal -- and summing over all lattice sites can counteract the decay. However, it does suppress all kinds of self-energy processes beyond the local term \begin{equation} \label{eq:dmft-self-energy-scaling} \Sigma_{ij} \stackrel{d\to\infty}{=} \delta_{ij} \Sigma_{ii} + \mathcal{O}(\sqrt{1/Z}^{d(i,j)})~. \end{equation} This can be shown via simple power-counting: any two internal vertices with a lattice index summation scale as $\mathcal{O}(Z^{d(i,j)})$ are connected by at least three Green functions, each scaling as~\eqref{eq:dmft-gf-scaling}. Thus, only $i=j$ contributions are non-vanishing. Despite the then \textit{local} nature of the self-energy, without further considerations, it would still be required to sum up all diagrams that could contribute to the local $\Sigma_{ii}$. However, via the same power-counting argument, and noting that, due to~\eqref{eq:dmft-gf-scaling} and~\eqref{eq:dmft-self-energy-scaling}, the lattice Green function must also scale as \begin{equation} G_{ij} \stackrel{d\to\infty}{=}\mathcal{O}(\sqrt{1/Z}^{d(i,j)})~, \end{equation} for self-energies derived through~\eqref{eq:2pi-self-energy}, only $\Gamma_2$ expansions with local Green functions are non-vanishing \begin{equation} \label{eq:dmft-gamma2} \Gamma_2\sbr{G_{ij}} \stackrel{d\to\infty}{=} \Gamma_2\sbr{G_{ii}}~. \end{equation} Hence, all contributions from~\eqref{eq:gamma2} have the same internal and external lattice indices, and the complex lattice interactions reduce to single-site interactions. This is the critical idea~\cite{Georges_1992} behind DMFT, where a lattice problem is mapped\footnote{Regarding the treatment of Anderson lattice models, the works~\cite{Grewe_1987, Kuramoto_1987} preceded the DMFT, ending up with descriptions of single-impurities interacting with an effective conduction sea. However, these overestimated the role of local interactions~\cite{Pruschke_1993} in comparison with DMFT, which is exact for any lattice problem with on-site interactions in infinite dimensions.} to an effective single-site problem. \subsection{The cavity construction} The cavity construction of DMFT is based on a formulation where the Dyson series is written in terms of Green functions where a single site is removed from the lattice. Consider a general lattice Hamiltonian $H(z) = \sum_{ij} \sbr{h_{ii}(z) \delta_{ij} + h_{ij}(z) \del{1 - \delta_{ij}}}$ composed of non-interacting and interacting local terms $h_{ii}(z) = h_{0,ii}(z) + h_{\mathrm{int},ii}(z)$, diagonal in real-space indices, as well as non-local non-interacting hopping terms $h_{ij}(z)$, with real-space matrix elements $t_{ij}$. By defining the \textit{single-site} Green function \begin{equation} \label{eq:single-site-gf} g^{-1}_{ij}(z,z') = \delta_{ij} \delta_\gamma(z,z') \sbr{i \partial_z - h_{ii}(z)} \equiv \delta_{ij} \cbr{\delta_\gamma(z,z') \sbr{i \partial_z - h_{0,ii}(z)} - \bar{\Sigma}_{ii}(z,z')}~, \end{equation} which is diagonal in real-space $g_{ij}(z,z') = \delta_{ij} g_{ii}(z,z')$ by definition, the Dyson series for the \textit{lattice} Green function of the system reads \begin{equation} G_{ij} = \delta_{ij} g_{ii} + g_{ii} \sum_{k_1 \ldots k_n} t_{i k_1} g_{k_2 k_2}\,t_{k_2 k_3}\,g_{k_3, k_3} \ldots \,t_{k_n j}\,g_{j j}~, \end{equation} where $\sum_{k_1 \ldots k_n} = \sum_{k_1} + \sum_{k_1 k_2} + \ldots$. The Dyson series is expressed in terms of the \textit{cavity} Green function $G^{[i]}_{jk}$, which is the lattice Green function with site $i$ removed \begin{equation} G_{ij} = \delta_{ij} g_{ii} + \underbrace{\del{\sum_{k_1 \ldots k_n} g_{ii}\,t_{i k_1 }\,g_{k_1 k_1}\,\ldots\,t_{k_n i}\,g_{i i}}}_{G_{ii}} \sum_{q_1 \neq i} t_{i q_1} \underbrace{\del{\sum_{q_2 \ldots q_n\neq i} g_{q_1 q_1}\,t_{q_1 q_2} \ldots \,g_{q_n q_n}}}_{G^{[i]}_{q_1 q_n}} \,t_{q_n j}\,g_{j j}~. \end{equation} Re-iterating the same procedure on the cavity Green function yields \begin{equation} \label{eq:dmft-hybridisation} G_{ij} = \delta_{ij} g_{ii} + G_{ii} \underbrace{\sum_{q_1} t_{i q_1} \del{\sum_{q_2 \ldots q_n} G^{[i]}_{q_1 q_1}\,t_{q_1 q_2}\,G^{[i, q_1]}_{q_2 q_2}\,t_{q_2 q_3}\,G^{[i, q_1, q_2]}_{q_3 q_3} \ldots}\,t_{q_n j}}_{\Delta_{ij}}\,g_{j j}~, \end{equation} for which the \textit{hybridisation} function $\Delta_{ii}$ enters the Dyson series for the \textit{local} Green function as \begin{equation} \label{eq:dmft-local-dyson} G_{ii} = g_{ii} + G_{ii}\,\Delta_{ii}\,g_{ii}~, \end{equation} and hence \begin{equation} \label{eq:dmft-local-gf} G_{ii}^{-1}(z,z') \stackrel{\mathrm{(DMFT)}}{=} \delta_\gamma(z,z') \sbr{i \partial_z - h_{0,ii}(z)} - \bar{\Sigma}_{ii}(z,z') - \Delta_{ii}(z,z')~, \end{equation} the local Green function is formulated as a local problem interacting with a renormalised sea. The hybridisation function $\Delta_{ii}$ contains the temporal \textit{correlations} of hopping from the site $i$ to the rest of the lattice and then back to site $i$ at a different time. Within DMFT, the lattice problem is reduced to a set of effective sites which act as independent local scattering centres~\cite{Grewe_1987}. The equivalence of the \textit{lattice} Green function \begin{equation} \label{eq:dmft-lattice-gf} \begin{split} G_{ij}^{-1}(z,z') &= \delta_\gamma(z,z') \sbr{\delta_{ij} \del{i \partial_z - h_{0,ii}(z)} - t_{ij}} - \Sigma_{ij}(z,z') \\ &\stackrel{\mathrm{(DMFT)}}{=} \delta_\gamma(z,z') \sbr{\delta_{ij} \del{i \partial_z - h_{0,ii}(z)} - t_{ij}} - \bar{\Sigma}_{ii}(z,z')~, \end{split} \end{equation} commonly expressed in the momentum basis, with~\eqref{eq:dmft-local-dyson} or~\eqref{eq:dmft-local-gf} composes the set of DMFT equations. Note that the self-energies $\Sigma_{ii}$ and $\bar{\Sigma}_{ii}$ of~\eqref{eq:dmft-lattice-gf} are two very different objects, which are only equivalent if~\eqref{eq:dmft-gamma2} holds: whereas the former is the local self-energy of a lattice problem at site $i$, the latter is the self-energy of a local problem defined at site $i$~\eqref{eq:single-site-gf}. \subsection{Formulation on the Bethe lattice} A Bethe lattice is an infinitely connected (cycle-free) graph (\cref{fig:bethe-lattice}) with connectivity~$Z$. The density of states (distribution of the momentum energy in crystal lattices) for a Bethe lattice with hopping matrix element $t$ and infinite connectivity \begin{equation} \rho_{Z\to\infty}(\varepsilon) = \frac{1}{\pi t} \sqrt{1 - \del{\frac{\varepsilon}{2t}}^2} \Theta(4t^2 - \varepsilon^2)~, \end{equation} or any $Z>2$, shares similar features with the ones found in three-dimensional lattices, namely the square-root behaviour near the band edges. Even though specific features of the density of states are of great importance in quantitative calculations, capturing its broader features is generally enough for qualitative analysis, which makes the Bethe lattice a suitable approximation for 3-dimensional systems~\cite{Economou2006}. \begin{figure}[!htb] \centering \includegraphics[width=0.5\figwidth]{diagrams/sketches/bethe-lattice} \caption{A Bethe lattice with connectivity $Z=3$. For $Z>2$, the ratio of boundary sites and the total number of sites becomes one when the latter approaches the thermodynamic limit. For this reason, the Bethe lattice is considered to be formed by the sites deep within the $Z$-Cayley tree, i.e., infinitely far away from the boundary. In this region, all sites become equivalent, with connectivity~$Z$.} \label{fig:bethe-lattice} \end{figure} A remarkable property of Bethe lattices is their self-similarity, resulting in equivalence between the \textit{cavity} and the \textit{local} Green function. Consider a Bethe lattice with a re-scaled~\eqref{eq:dmft-rescale} hopping matrix element $t$ and the contour-ordered local Green function~\eqref{eq:dmft-local-gf} at the site~$i$ \begin{equation} \label{eq:dmft-bethe-local} G_{ii}^{-1}(z,z') = \delta_\gamma (z,z')i \partial_z - \bar{\Sigma}_{ii}(z,z') - Z \frac{t}{\sqrt{Z}} G^{[i]}_{j_\alpha j_\alpha}(z,z') \frac{t}{\sqrt{Z}}~, \end{equation} Due to being a cycle-free graph, the hybridisation function $\Delta_{ii}(z,z')$~\eqref{eq:dmft-hybridisation} simplifies since there is just one possible path of connecting the site $i$ to each of its $Z$ neighbours. The cavity Green function reads \begin{equation} \label{eq:dmft-bethe-weiss} G^{{[i]}^{-1}}_{j_\alpha j_\alpha}(z,z') = \delta_\gamma (z,z')i \partial_z - \bar{\Sigma}_{j_\alpha j_\alpha}(z,z') - \del{Z-1}\frac{t}{\sqrt{Z}} G^{[i,j_\alpha]}_{k_\alpha k_\alpha}(z,z') \frac{t}{\sqrt{Z}}~. \end{equation} In the $Z\to\infty$ limit, due to the self-similarity and cycle-free properties of the lattice, $G^{[i]}_{j_\alpha j_\alpha}(z,z') = G^{[i,j_\alpha]}_{k_\alpha k_\alpha}(z,z')$ and~\eqref{eq:dmft-bethe-weiss} can be inverted \begin{equation} G^{{[i]}}_{j_\alpha j_\alpha}(z,z') = \frac{y \pm \sqrt{y^2 - 4t^2}}{2 t^2}~,\quad\text{with }y=i \partial_z - \bar{\Sigma}(z,z')~, \end{equation} where $\bar{\Sigma}_{ii}(z,z') = \bar{\Sigma}_{j_\alpha j_\alpha}(z,z') = \bar{\Sigma}(z,z')$. Inserting this result back in~\eqref{eq:dmft-bethe-local} yields \begin{equation} G_{ii}(z,z') = G^{{[i]}}_{j_\alpha j_\alpha}(z,z')~, \end{equation} and, ultimately, a closed expression for the hybridisation function \begin{equation} \label{eq:bethe-hybridisation} \Delta_{ii}(z,z') = t\,G_{ii}(z,z')\,t~, \end{equation} thus reducing the DMFT equations to a single-equation for the local problem~\eqref{eq:dmft-local-gf}. \section{The non-crossing approximation} \label{sec:nca} Given that in infinite dimensions the lattice problem is reduced to a local problem~\eqref{eq:dmft-gamma2}, the lowest-order terms of the 2PI loop expansion~\eqref{eq:gamma2} of~\eqref{eq:full-model} and~\eqref{eq:photonic-bath-action} read \begin{equation} \label{eq:nca-gamma2} \begin{split} \Gamma_2 \stackrel{\mathrm{(DMFT)}}{=} &-i\int \dif z \dif z' \sum_{\sigma} i G_{c_\sigma}(z,z') i G_b(z,z') i G_{f_\sigma}(z',z) \sbr{V_0^2 + i g_0^2 \Xi_a(z,z')} \\ &= -i \sum_{\sigma} \sbr{ \vcenter{\hbox{\input{diagrams/feynman-diagrams/gamma2_nophoton}}} + \vcenter{\hbox{\input{diagrams/feynman-diagrams/gamma2_photon}}} }~. \end{split} \end{equation} The driven-dissipative lattice problem is hence reduced to an effective single-impurity Anderson problem, where the bare coupling vertex $V_0^2$ becomes dressed by interactions with the electromagnetic field \begin{equation} \label{eq:vertex-dressing} V_0^2 \to V_0^2 + i g_0^2 \Xi_a(z, z') = V_0^2 + i g_0^2 \sbr{G_a(z, z') + G_a(z', z) + \eta \Pi_a(z, z') + \eta \Pi_a(z', z)}~, \end{equation} containing the incident pulse and photonic reservoir dynamics. %Note that $G_c$ also contains a bath-induced hybridisation~\eqref{eq:fermionic-bath-action} as well as the DMFT hybridisation~\eqref{eq:dmft-local-gf} functions. Due to the absence of \textit{crossing} lines, the approximation was coined~\cite{Kuramoto_1984} the non-crossing approximation (NCA\footnote{This technique is also cynically known as the \enquote{Never Correct Approximation}, due to its shortcomings versus behemoths such as Quantum Monte Carlo. However, as foretold by Philip Anderson~\cite{Jones_2015} \begin{displayquote} The better the machinery, the more likely it is to conceal the workings of nature, in the sense that it simply gives you the experimental answer without telling you why the experimental answer is true \end{displayquote} the weaknesses of the NCA in \textit{exact} numerics strengthen the fine control and understanding of the underlying physics. And despite requiring a numerical solution due to its non-linear nature, its solutions are well understood, and extensions beyond the NCA have shown good quantitative agreement with exact methods.}). This approximation is the simplest term arising from the loop expansion. It involves the maximum number of intermediate states for each order of the perturbation~\cite{Kuramoto_1984} due to the non-crossing lines. Similarly to the dimensionality scaling of DMFT, the NCA becomes exact in the limit of infinite $f$-level degeneracy $N$. By rescaling the hybridisation strength, $V \to V/\sqrt{N}$, as the number of crossing lines increases, less intermediate spin states are summed over, and the $\frac{1}{\sqrt{N}^n}$ becomes dominant for an $n$-order expansion. The NCA qualitatively predicts the emergence of the Kondo effect and some properties of the Kondo resonance~\cite{Grewe_1987} and is accurate for high-energy features. However, it does not fully capture all low energy properties and can pathologically diverge at the Fermi energy in the limit of vanishing temperature. \subsection{\texorpdfstring{$\zeta$}{Zeta}-scaling of contour-ordered self-energies} Due to the projection requirements (\cref{sec:ap-zeta-scaling}), the Green functions are to be calculated in the $\zeta \to 0$ limit of the enlarged Hilbert space. The self-energies are obtained by taking the functional derivative~\eqref{eq:2pi-self-energy} of~\eqref{eq:nca-gamma2} \begin{equation} \label{eq:nca-self-energies} \begin{split} \lim_{\zeta\to0}\Sigma_{b_\zeta}(z, z') &= \lim_{\zeta\to0} -i\, \sbr{V_0^2 + i g_0^2 \Xi(z,z')} \sum_\sigma G_{f_{\zeta,\sigma}}(z, z') G_{c_{\zeta,\sigma}}(z', z)\\ \lim_{\zeta\to0}\Sigma_{f_{\zeta,\sigma}}(z, z') &= \lim_{\zeta\to0} +i\, \sbr{V_0^2 + i g_0^2 \Xi(z,z')} G_{c_{\zeta,\sigma}}(z, z') G_{b_\zeta}(z, z')\\ \lim_{\zeta\to0}\Sigma_{c_{\zeta,\sigma}}(z, z') &= \lim_{\zeta\to0} +i\, \sbr{V_0^2 + i g_0^2 \Xi(z,z')} G_{f_{\zeta,\sigma}}(z, z') G_{b_\zeta}(z', z)\\ \lim_{\zeta\to0}\Sigma_{a_\zeta}(z, z') &= \lim_{\zeta\to0}-i\, g_0^2 \sum_{\sigma} \cbr{G_{c_{\zeta,\sigma}}(z,z') \sbr{i G_{f_{\zeta,\sigma}}(z',z)G_{b_\zeta}(z,z')} + \del{(z,z') \to (z',z)}}~. \end{split} \end{equation} Due to the self-energies of the conduction electrons $c$ and the pulse-photons $a$ containing auxiliary-particle loops, their components have an \textit{additional} $\zeta$ pre-factor in comparison with the auxiliary-particle self-energies \begin{equation} \label{eq:self-energy-zeta-scaling} \begin{split} \lim_{\zeta\to0} \Sigma_{b_\zeta}(z, z') &\sim \Theta_\gamma(z,z') \mathcal{O}(1) + \Theta_\gamma(z',z) \mathcal{O}(\zeta)\\ \lim_{\zeta\to0} \Sigma_{f_{\zeta,\sigma}}(z, z') &\sim \Theta_\gamma(z,z') \mathcal{O}(1) + \Theta_\gamma(z',z) \mathcal{O}(\zeta)\\ \lim_{\zeta\to0} \Sigma_{c_{\zeta,\sigma}}(z, z') &\sim \mathcal{O}(\zeta)\\ \lim_{\zeta\to0} \Sigma_{a_\zeta}(z, z') &\sim \mathcal{O}(\zeta)~. \end{split} \end{equation} Since the conduction and photon-pulse Green functions do not share the same scaling as auxiliary-particle Green functions and scale as $\mathcal{O}(1)$, their self-energies arising from~\eqref{eq:nca-gamma2} are vanishing in the $\zeta\to0$ limit. Consequently, the $c$ and $a$ Green functions are not \textit{renormalised} by local interactions with the auxiliary particles. Note that for conduction electrons $c$ within DMFT, despite the \textit{local} self-energy vanishing, there are still temporal fluctuations encoded in the hybridisation function $\Delta(z,z')$ that may renormalise $\lim_{\zeta\to0} G_{c_{\zeta,\sigma}}(z,z')$. For the photon pulse, unlike in a closed cavity where coherent excitations between the electronic system and the electromagnetic field can build up, in free space, the photon pulse propagates past the lattice, interacts with it and flies away. Physically, it is then expected that the pulse photons are not renormalised due to the lack of a DMFT-like procedure, as there is no mechanism to build up temporal correlations. Therefore, the photon-pulse Green function is reduced to the non-interacting Green function (\cref{sec:travelling-pulse}) $\lim_{\zeta\to0} G_{a_\zeta}(z,z') \to G_{a_0}(z,z')$. Note, however, that this only holds in the $\zeta\to0$ limit, which is used for calculations. The conduction and the photon-pulse Green functions will be renormalised locally by interactions with \textit{physical} particles. However, these renormalisation effects, encoded in the local self-energies, do not play a role in the $\zeta\to0$ limit. \subsection{Projection of physical observables} \label{sec:nca-projection} Apart from the self-energies generated by the approximation, one must also formulate how to project to the original Hilbert space within the NCA~\eqref{eq:nca-gamma2}. The Green function for physical $f$ electrons can be obtained by the projection~\eqref{eq:projection} of the auxiliary particles \begin{equation} \label{eq:projected-f-gf} \begin{split} G_{f_\sigma}(z,z') \coloneqq -i \ev{\mathds{T}_\gamma \hat f^\dagger_\sigma(z) \hat f_\sigma(z')}_{\mathrm{physical}} &\equiv \lim_{\zeta\to0} \frac{-i \ev{\mathds{T}_\gamma \hat f^\dagger_\sigma(z) \hat b(z) \hat b^\dagger(z') \hat f_\sigma(z)}_\zeta}{\ev{\hat{Q}}_\zeta} \\ &\stackrel{\mathrm{(NCA)}}{=} % -i \lim_{\zeta\to0} \frac{\ev{\hat f^\dagger(z) \hat f(z')}_\zeta \ev{\hat b(z) \hat b^\dagger(z')}_\zeta} {\ev{\hat{Q}}_\zeta} = \lim_{\zeta\to0} \frac{i G_{f_{\zeta,\sigma}}(z,z')G_{b_\zeta}(z',z)}{\ev{Q}_\zeta}~. \end{split} \end{equation} This is because the auxiliary-particle 4-point function factorises at the NCA level \begin{equation} \ev{\mathds{T}_\gamma \hat f_\sigma^\dagger(z) \hat b(z) \hat b^\dagger(z') \hat f_\sigma(z)}_\zeta = \input{diagrams/feynman-diagrams/gd-projection}~, \end{equation} since the non-trivial \textit{connected} vertex contributions, described by a Bethe-Salpeter equation\footnote{The Bethe-Salpeter equation is the analogue of a Dyson series for a 4-point function $X$ \begin{equation*} X = X_0 + X_0\, \Gamma_\mathrm{irreducible}\, X~, \end{equation*} where $X_0$ denotes a non-interacting ensemble average of the 4-point function and $\Gamma_\mathrm{irreducible}$ is the 4-point \textit{irreducible} vertex of the theory. Similar to a t-matrix series~\eqref{eq:t-matrix}, the Bethe-Salpeter equation can be related to a \textit{reducible} vertex $\Gamma_\mathrm{reducible}$ via \begin{equation*} X = X_0 + X_0\, \Gamma_\mathrm{reducible}\, X_0~, \end{equation*} where \begin{equation*} \Gamma_\mathrm{reducible} = \Gamma_\mathrm{irreducible} + \Gamma_\mathrm{irreducible}\, X_0\, \Gamma_\mathrm{reducible}~. \end{equation*} For more details, refer to~\cite{Abrikosov1975, Haussmann_1999, Bickers2004}.} containing \textit{crossing} diagrams beyond the NCA~\eqref{eq:nca-gamma2} \begin{equation} \input{diagrams/feynman-diagrams/gd-projection-vertex}~. \end{equation} The literature~\cite{Kroha2004} had taken a different route for calculating other physical observables. For an Anderson lattice model without photon coupling, the only other physical observables are related to the conduction electron. In this model, there is an exact known expression for the t-matrix $T_{c_\sigma}$ of the local conduction electrons \begin{equation} \label{eq:t-matrix} G_{c_\sigma}(z, z') = \lim_{\zeta\to0} G_{c_{\zeta,\sigma}}(z,z') + \sbr{\del{\lim_{\zeta\to0} G_{c_{\zeta,\sigma}}} * T_{c_\sigma} * \del{\lim_{\zeta\to0} G_{c_{\zeta,\sigma}}}}(z, z')~, \end{equation} with $T_{c_\sigma} = V_0^2 G_{f_\sigma}$. However, the situation is far more complicated when photons are also coupled to the system, as it is unclear how to obtain a similar t-matrix expression for conduction electron or photon observables compatible with the NCA~\eqref{eq:nca-gamma2}. Noting that $\hat Q$ does \textit{not} factorise in interacting theories, these observables can be calculated via the \textit{t-matrix} formulation of the Bethe-Salpeter equation and the projection~\eqref{eq:projection-Q} \begin{equation} \label{eq:nca-c-projected} \begin{split} G_{c_\sigma}(z,z') &\equiv -i \lim_{\zeta\to0} \frac{\ev{\mathds{T}_\gamma \hat c^\dagger_\sigma(z) \hat c_\sigma(z') \hat Q}_\zeta}{\ev{\hat Q}_\zeta} \\ &\stackrel{\mathrm{(NCA)}}{=} \lim_{\zeta\to0} \cbr{G_{c_{\zeta,\sigma}}(z,z') + \sbr{G_{c_{\zeta,\sigma}} * \frac{\Sigma_{c_{\zeta,\sigma}}}{\ev{\hat{Q}}_\zeta} * G_{c_{\zeta,\sigma}}}(z,z')}~. % &\stackrel{\mathrm{(NCA)}}{=} \lim_{\zeta\to0} \sbr{\ev{\hat c(z) \hat c^\dagger(z')}_\zeta + \int_\gamma \dif \bar{z}\,\dif\bar{z}'\,\ev{\hat c(z) \hat c^\dagger(z')}_\zeta \frac{\Sigma_{c_\zeta}(\bar{z},\bar{z}')}{\ev{\hat{Q}}_\zeta}\ev{\hat c(\bar{z}') \hat c^\dagger(z')}_\zeta}~, \end{split} \end{equation} The t-matrix formulation Bethe-Salpeter equation for the 4-point function \begin{equation} \ev{\mathds{T}_\gamma \hat c^\dagger(z) \hat c(z') \hat Q}_\zeta = \input{diagrams/feynman-diagrams/gc-projection}~, \end{equation} is a function of the Bethe-Salpeter equation for the 4-point \textit{reducible} vertex \begin{equation} \input{diagrams/feynman-diagrams/gc-vertex-reducible}~. \end{equation} Due to $\ev{\hat Q}_\zeta \sim \mathcal{O}(\zeta)$, the denominator in~\eqref{eq:nca-c-projected} cancels at most one $Q$ loop. Auxiliary-particle loops also scale as $\mathcal{O}(\zeta)$ and hence, to leading order of the $\zeta$-expansion, the reducible vertex reduces to \begin{equation} \input{diagrams/feynman-diagrams/gc-vertex-irreducible}~, \end{equation} where the $Q$ insertion is dropped by noting that the resulting interaction, composed of contractions of auxiliary-particle operators with $\hat Q$, is equivalent to the same contraction without $\hat Q$~\eqref{eq:projection}. The projected photon-pulse Green function is similar to~\eqref{eq:nca-c-projected} and reads \begin{equation} \label{eq:nca-a-projected} G_a(z,z') \stackrel{\mathrm{(NCA)}}{=} \lim_{\zeta\to0} \cbr{G_{a_\zeta}(z,z') + \sbr{G_{a_\zeta} * \frac{\Sigma_{a_\zeta}}{\ev{\hat{Q}}_\zeta} * G_{a_\zeta}}(z,z')}~. \end{equation} \subsection{Dynamical mean-field theory and auxiliary particles} \label{sec:dmft-aux-particles} The usual non-equilibrium DMFT formulation results in a set of integral equations~\cite{Aoki2014}. However, semi-analytic methods such as the NCA make it possible to express these equations in a Kadanoff-Baym form, despite generally not being possible to obtain an expression for the local self-energies -- it would be required to find an expression for their projection (\cref{sec:nca-projection}). For the sake of brevity, it is assumed that the $f$ electrons do not hop, and all Green functions concerning the DMFT are related to the $c$ conduction electrons. Since auxiliary-particle observables are calculated as grand-canonical ensemble averages in the $\zeta\to0$ limit (\cref{sec:ap-zeta-scaling}), the DMFT-derived observables entering such calculations must be considered in the same limit. \paragraph{Bethe-lattice DMFT Kadanoff-Baym form} In a Bethe lattice, the explicit dependency of the hybridisation can be removed through~\eqref{eq:bethe-hybridisation}. For an effective single-impurity Anderson model defined at some site $i$, the auxiliary particles are defined \textit{only} at said site. By definition, the hybridisation function $\Delta_{c_{ii\sigma}}$ is traced over all states \textit{not} pertaining to site $i$, and hence, the canonical ensemble average of $\Delta_{c_{ii\sigma}}$ is the same as the grand-canonical one \begin{equation} \lim_{\zeta\to0} \Delta_{c_{\zeta,ii\sigma}}(z,z') = \Delta_{c_{ii\sigma}}(z,z')~. \end{equation} Assuming translational symmetry, the site index subscripts can now be dropped. The equation of motion for the local conduction electron grand-canonical Green function that enters the NCA self-energies~\eqref{eq:nca-self-energies} is directly obtained from~\eqref{eq:dmft-local-gf} \begin{equation} \label{eq:dmft-cavity-gf} \sbr{i \partial_z - h_{c_{0}}(z)} \del{\lim_{\zeta\to0} G_{c_{\zeta,\sigma}}(z,z')} = \delta_\gamma(z,z') + \sbr{t\,G_{c_\sigma}\,t * \del{\lim_{\zeta\to0} G_{c_{\zeta,\sigma}}}}(z,z')~, \end{equation} where the self-energy of the local problem is vanishing in the $\zeta\to0$ limit~\eqref{eq:self-energy-zeta-scaling}. The equation that generates the canonical local conduction electron Green function is obtained by convolving the local Green function~\eqref{eq:nca-c-projected} with the inverse Green function of~\eqref{eq:dmft-cavity-gf}, resulting in \begin{equation} \label{eq:dmft-local-gf-2} \sbr{i \partial_z - h_{c_{0}}(z)} G_{c_\sigma}(z,z') = \delta_\gamma(z,z') + \sbr{\lim_{\zeta\to0}\del{\frac{\bar{\Sigma}_{c_{\zeta,\sigma}}}{\ev{Q}_\zeta} * G_{c_{\zeta,\sigma}}} + t\,G_{c_{\sigma}}\,t * G_{c_{\sigma}}}(z,z')~. \end{equation} \paragraph{General DMFT Kadanoff-Baym form} For a general problem in a crystal lattice, consider the equation of motion of the lattice Green function~\eqref{eq:dmft-lattice-gf} in the momentum basis \begin{equation} \label{eq:dmft-lattice-gf-k} \sbr{i \partial_z - \varepsilon^c_{\bm{k}}} G_{c_{\bm{k}\sigma}}(z,z') = \delta_\gamma(z,z') + \sbr{\bar{\Sigma}_{c_\sigma} * G_{c_{\bm{k}\sigma}}}(z,z')~, \end{equation} where tight-binding was assumed $\varepsilon^c_{\bm{k}} = \sum_{ij} t^c_{ij} e^{i \bm{k} \cdot (\bm{R}_i - \bm{R}_j)}$, $h_{c_{0}}(z) = 0$. The local Green function is obtained by summing\footnote{ Discretising the momenta creates an infrared momentum cut-off, equivalent to considering a finite but periodic lattice. Despite such discretisations being far away from the discretisation found in the real world, where typically one finds $10^{23}$ cm$^{-3}$ atoms, it is generally good enough to capture most physics, including phase transitions -- which are defined only in the thermodynamic limit. However, when summing up plane-wave-like states, such as when calculating local Green functions, one can find an unexpected coherence revival at long times. Consider summing $n$ 1-dimensional plane waves distributed over the first Brillouin zone equidistantly \begin{equation*} \frac{1}{n} \sum_{k \in \underbrace{\sbr{-\pi, \ldots, \pi}}_{n\textrm{ elements}}} e^{i k \tau} = \cos\del{\pi \tau} + \cot\del{\frac{\pi \tau}{n-1}}\sin\del{\pi \tau}~. \end{equation*} Due to the finite number of $k$ states, the plane waves become resonant and return to their initial state after a \textit{finite} time, a phenomenon known as the \textit{Poincaré recurrence time}.} all the $\bm{k}$-modes \begin{equation} G_{c_{\sigma}}(z,z') = \sum_{\bm{k}} G_{c_{\bm{k}\sigma}}(z,z')~, %e^{-i \bm{k} \cdot \bm{R}_i}~, \end{equation} and the equation of motion for the local conduction electron grand-canonical Green function, that enters the NCA self-energies~\eqref{eq:nca-self-energies}, is obtained by the Dyson series of the local conduction electron Green function \begin{equation} i \partial_z \del{\lim_{\zeta\to0}G_{c_{\zeta,\sigma}}(z,z')} = i \partial_z \sbr{G_{c_{\sigma}} * \del{\delta_\gamma - \bar{\Sigma}_{c_\sigma} * \del{\lim_{\zeta\to0}{G_{c_{\zeta,\sigma}}}}}}(z,z')~. \end{equation} Finally, the self-energy is either obtained by the relation between the t-matrix equation (cf.~\eqref{eq:t-matrix}) and the Dyson series \begin{equation} \label{eq:dmft-self-energy} \sbr{T_{c_\sigma} * \del{\lim_{\zeta\to0} G_{c_{\zeta,\sigma}}}}(z,z') = \sbr{\bar{\Sigma}_{c_\sigma} * G_{c_{\sigma}}}(z, z')~, \end{equation} or, e.g., for $\bar{\Sigma}_{c_\sigma}(z,z') = V_0 \, g_{f_\sigma}(z,z') \, V_0$, by differentiating the Dyson series for $f$ electrons \begin{equation} i \partial_z \bar{\Sigma}_{c_\sigma}(z,z') = i \partial_z \sbr{\del{V_0\, G_{f_\sigma}\, V_0} * \del{\delta_\gamma - G_{c_{\sigma}} * \bar{\Sigma}_{c_\sigma}}}(z,z')~. \end{equation} \section{Numerical procedure} \label{sec:numerical-procedure} The non-equilibrium driven-dissipative lattice is solved at the level of 2-point functions (\cref{sec:neqft-gf}), and all presented results were obtained using the same numerical procedure -- the only difference being the model parameters. The auxiliary-particle Hamiltonian~\eqref{eq:full-model} is solved in a Bethe lattice at the level of DMFT (\cref{sec:dmft}) and NCA (\cref{sec:nca}), and is always coupled to fermionic and photonic heat baths (\cref{sec:heat-baths}) as well as the photon-pulse mode (\cref{sec:travelling-pulse}). Because of the auxiliary particles, the 2-point functions are calculated in the $\zeta\to0$ limit (\cref{sec:ap-zeta-scaling}) in both thermal equilibrium and non-equilibrium regimes and only after time-integration are projected to the physical Hilbert space. Solutions of the 2-point functions in thermal equilibrium are used as initial conditions (\cref{sec:gf-initial-conditions}) for the non-equilibrium time-evolution. Hence, both regimes are described by the same effective action, with the only difference being in the time contour (\cref{sec:contours}). This is then reflected in the resulting equations of motion of the 2-point functions and the occupation of the photon pulse, which is taken as constant for all times before the beginning of the time evolution, at time $t=t_0$. \subsubsection{Solutions in thermal equilibrium} Due to the inherent time-translational invariance of systems in thermal equilibrium, the equations of motion of the 2-point functions are transformed into a problem in real-frequency, resulting in a series of self-consistent Dyson equations for auxiliary-particle (\cref{sec:project-dyson}) and DMFT observables (\cref{sec:dmft-aux-particles}). Due to the fluctuation-dissipation relation (\cref{sec:neqft-interpretability}), thermal distributions can be enforced for \textit{all} observables. However, to ensure that the heat baths can thermalise the system, the distributions are enforced \textit{only} for bath observables. This safeguard ensures the system can thermalise \textit{numerically}, since too small bath coupling strengths can be washed out due to finite precision and tolerances of the equilibrium and non-equilibrium solvers. Note that the occupation of the fields can be indeterminate when solving self-consistent equations without enforcing distribution functions. For example, the occupation of the electrons is only (indirectly) fixed if the system is coupled to a fermionic heat bat. Otherwise, there is no mechanism to fix the chemical potential (this problem is not present in the time evolution since the occupation number is determined by the system's initial condition and conserved by the dynamics generated by the Hamiltonian). The self-consistent equations are solved via a simple fixed-point iteration scheme, namely Anderson mixing~\cite{Walker_2011}, which has shown to be quite robust, fast, and more than adequate for this class of problems. Convergence of the self-consistent equations is achieved when the infinity-norm between the residuals is below the square of the tolerances of the non-equilibrium solver. \paragraph{Leveraging time and frequency representations in thermal equilibrium} Problems of the NCA family in thermal equilibrium are commonly solved on a fixed basis, typically in real-frequency or imaginary time. However, the locality of the Dyson equations in frequency or of the self-energies in time can be maximally exploited with a Fourier transformation between bases. Solutions for the NCA in thermal equilibrium are hence be found by finding the fixed-point of \begin{enumerate} \item Compute the Dyson equations~\eqref{eq:dyson-real-freq-aux} and~\eqref{eq:lesser-real-freq-aux} in real frequency, \item Inverse Fourier transform~\eqref{eq:wigner-ville-inv} the Green functions to relative time, \item Compute the self-energies~\eqref{eq:nca-self-energies} in relative time, \item Fourier transform the self-energies to real frequency. \end{enumerate} Finding solutions of NCA equations in real frequency is notoriously difficult~\cite{Kroha2004, Costi_1996} due to the sharp features of the functions at low temperatures. Finding solutions in their Fourier dual space is, however, a more stable procedure as the highly-peaked functions in frequency are transformed into slow decaying functions in relative time. Furthermore, relations such as Kramers-Kronig \begin{subequations} \begin{align} \begin{split} G^R(\tau) &= % \int \frac{\dif \omega}{2\pi} e^{-i \omega \tau} G^R(\omega) % = \int \frac{\dif \omega}{2\pi} e^{-i \omega \tau} \int \frac{\dif \varepsilon}{\pi} \frac{- \im G^R(\varepsilon)}{\omega - \varepsilon + i \eta} % = % - \int \frac{\dif \omega}{2\pi i} \frac{e^{-i \omega \tau}}{\omega + i \eta} 2i \int \frac{\dif \varepsilon}{2\pi} e^{-i \varepsilon \tau} \im G^R(\varepsilon) = \Theta(\tau) \int \frac{\dif \varepsilon}{2\pi} e^{-i \varepsilon \tau} \sbr{2i \im G^R(\varepsilon)} \end{split} \\ \begin{split} G^R(\omega) &= \int \dif \tau\, e^{i \omega \tau} \Theta(\tau) \int \frac{\dif \varepsilon}{2\pi} e^{-i \varepsilon \tau} \sbr{2i \im G^R(\varepsilon)}~. \end{split} \end{align} \end{subequations} are elegantly satisfied when changing between bases. \subsubsection{Solutions in non-equilibrium} Instead of using mixed Green functions (\cref{sec:neqft-components}) to encode vestigial thermal correlations, the thermal equilibrium Green functions are inverse Wigner-Ville FFT'ed\footnote{For a function $f(x)$ discretised on an equidistant grid $x_i = x_0 + i \Delta x$, for $i \in \sbr{0, n-1}$, its discretised Fourier transform (cf.~\eqref{eq:wigner-ville-inv}) \begin{equation*} \hat{f}(\xi_k) = \frac{\Delta x}{2\pi} e^{-i x_0 \xi_k} \operatorname{FFT}\del{f(x)}_k~, \end{equation*} where $\xi_k = \frac{2\pi}{\Delta x} k / n$, can efficiently be calculated with recourse to a fast Fourier transformation (FFT) \begin{equation*} \operatorname{FFT}(f(x))_k \coloneqq \sum_{j=0}^{n-1} f(x_j) e^{-2\pi i j k / n}~. \end{equation*}} (\cref{sec:gf-initial-conditions}), and encode the system in thermal equilibrium for all times \textit{before} the beginning of the time-integration. The resulting initial conditions for the Green functions have an infinite size and must be truncated at some time $t \ll t_0$. The truncation is only validated if the right-hand side of the diagonal equations of motion (\cref{sec:stepping-scheme}) at $t = t_0$ are approximately zero within the numerical tolerances of the non-equilibrium solver. This ensures that despite possibly not retaining the whole information of the system in thermal equilibrium (at the level of 2-point functions), enough information is retained such that the system is integrated as if it had started with the full information. Another possible validation is confirming that the system stays in thermal equilibrium for all times $t > t_0$ when the coupling strength to the driving field is zero. The non-equilibrium equations of motion for auxiliary-particle (\cref{sec:project-kb}) and DMFT observables (\cref{sec:dmft-aux-particles}) are Kadanoff-Baym equations (\cref{sec:dyson-kbe}) which are solved via an adaptive Kadanoff-Baym solver~\cite{Meirinhos_2022} with tolerances $\texttt{rtol}=10^{-5}$ and $\texttt{atol}=10^{-8}$ (\cref{sec:adams}). \paragraph{Why truncation is ruled out} Truncation of the integral kernels (\cref{sec:memory-truncation}) found in the Kadanoff-Baym equations of motion~\eqref{eq:kb} can greatly reduce the computational effort, however, are not viable for the total integration-time of interest $t_\mathrm{max}$. In quantum many-body systems, temperature provides a natural cut-off timescale \begin{equation} \tau_\mathrm{cutoff} \sim \frac{\hbar}{k_\mathrm{B}T}~, \end{equation} beyond which coherent quantum processes are washed out by thermal fluctuations~\cite{Coleman_2001}. Given the Kondo lattice coherence developing only when the system is at a temperature below the Kondo coherence temperature (\cref{sec:kondo-pheno}) and $t_\mathrm{max} \sim \frac{\hbar}{k_\mathrm{B} T^*_K}$, any advantageous truncation would spuriously induce high-temperature effects. Note that for very low temperatures -- lower than the smallest intrinsic energy scale of the problem, the Green function components displaying the most extended tails in the relative-time direction (\cref{sec:wigner-basis}) are not spectral components but components related with the occupation of particles. From the viewpoint of thermal equilibrium, this can be understood by the sharp features of the distribution functions at low temperatures in frequency being translated to very long tails in the dual (relative-)time basis. \subsubsection{Projection and Wigner transformation} After time integration, the physical Green functions are computed by projecting the auxiliary particles out (\cref{sec:nca-projection}). A naive application of the Langreth rules (\cref{sec:langreths-rules}) would be very time-consuming, and these are accelerated by expressing the time-integrals as matrix products. The Green functions are then linearly interpolated into an equidistant time grid and Wigner-Ville FFT'ed. The interpolation is the only step in the whole procedure that could introduce significant errors in the data, typically as spurious small oscillatory terms, which are easily identified visually. \section{Kondo collapse and revival by pulsed light} Consider~\eqref{eq:full-model}, \eqref{eq:fermionic-bath-action} and \eqref{eq:photonic-bath-action} with the parameters of \cref{tab:parameters}. These parameters constitute the base of the following analyses and are kept constant throughout unless explicitly denoted. A large value of the inverse temperature $\beta$ and the $f$-electron ground-state energy $\varepsilon^f_0$ and hybridisation $V_0$ are associated with the Kondo regime ($J_K \sim 0.5$). Due to the DMFT freezing out all spatial fluctuations, the system is in a heavy Fermi liquid phase, very far away from magnetic instabilities brought by criticality (\cref{fig:qcp}). The system is dipole-coupled to the electromagnetic field with coupling strength $g_0$. The coupling to the electromagnetic-field vacuum modes is $\eta$ times stronger than the coupling to the external pulse, where $\Lambda$ is a soft cut-off of the coupling to high-energy vacuum modes. The coupling strengths to the fermionic and photonic reservoirs, $\alpha$ and $\sqrt{\eta} g_0$, respectively, are small enough that the system's qualitative properties remain unchanged. For example, too large $\alpha$ can change the quasiparticles' lifetime and wash out the Kondo lattice coherence of the conduction electrons, and too large $\sqrt{\eta} g_0$ can induce a significant Lamb shift of the single-particle peak at $\varepsilon^f_0$. However, $\alpha$ is chosen sufficiently large for thermalisation times -- associated with the timescale $\sim \frac{\pi \upsilon}{\alpha^2}$ -- to be accessible through numerical integration, where $\upsilon$ is the Bethe lattice hopping. The external pulse has a small central frequency $\omega_0$ and frequency-bandwidth $\Omega_0$, both associated with the THz regime, and a maximum photon intensity of $\bar{n}_a$. The effective coupling strength of the external pulse to matter, $g_0^2 \bar{n}_a$, is small since the non-equilibrium dynamics are intended not to be very violent. \begin{table}[!htb] \centering \resizebox{\linewidth}{!}{% \begin{tabular}{|clllllllll|} \hline \multicolumn{2}{|c|}{Anderson model} &\multicolumn{2}{c|}{Fermionic reservoir} &\multicolumn{3}{c|}{Electromagnetic field} &\multicolumn{3}{c|}{External pulse} \\ \hline \multicolumn{1}{|l}{$\varepsilon^f_0 / \upsilon =-0.35$} &\multicolumn{1}{l|}{$V_0 / \upsilon = 0.3$} &$\beta \upsilon = 150$ &\multicolumn{1}{l|}{$\alpha / \upsilon = 0.04$} &$g_0 / \upsilon = 0.04$ &$\eta = 10$ &\multicolumn{1}{l|}{$\Lambda / \upsilon = 0.25$} &$\omega_0 / \upsilon= 0.001$ &$\Omega_0 / \upsilon = 0.06$ &$\bar{n}_a = 10$ \\ \hline \end{tabular} } \caption{Parameters in units of the Bethe lattice hopping $\upsilon$.} \label{tab:parameters} \end{table} To understand the collapse\footnote{There had been some previous attempts~\cite{Fauseweh_2020, Zhu_2021} at observing some sort of Kondo collapse on a Kondo lattice driven by a strong pulse of radiation. However, the studies suffered several limitations. First, the Kondo lattice model contains no charge fluctuations, which, according to the investigations of~\cref{sec:thz-mf}, appear critical in non-equilibrium dynamics. Second, long-time dynamics could not be investigated since only short-time evolutions could be resolved. Additionally, only large pulse frequencies could be considered, which are too energetic and strongly excite the system, which ends up at a very high temperature due to being a closed formulation. Lastly, the light was treated classically and hence could not induce intra-band transitions and spontaneous emission into experimentally-accessibly electromagnetic field modes.} and revival of the Kondo coherence, one must first know how to identify it. In a system with many $f$ sites, below the Kondo coherence temperature $T^*_K$, the $f$ electrons which scatter into the conduction band and later scatter back into the same site have partaken in the build-up of Kondo resonances (\cref{sec:kondo-pheno}) in other lattice sites. Such excitations can be detected by an enhancement of the local density-of-states of the \textit{cavity} Green function~\eqref{eq:dmft-cavity-gf} near the Fermi energy, which influences the Kondo temperature~\eqref{eq:kondo-temp} exponentially\footnote{Within DMFT,~\eqref{eq:kondo-temp} could be used to estimate the Kondo coherence temperature, as the problem is reduced to a single-site problem -- with the non-interacting density of states entering the equation being replaced the density of states of the \textit{cavity} Green function~\eqref{eq:dmft-cavity-gf}. However, this is a crude approximation since the cavity density of states is energy- and temperature-dependent and all \textit{spatial} fluctuations were neglected. Within this approximation, the Kondo lattice-coherence temperature $T^*_K$ is expected to be \textit{higher} than the single-impurity $T_K$ due to the spectral enhancement of the cavity Green function at the Fermi energy, in agreement with several compounds~\cite{Kirchner_2020}}. Similarly, the \textit{local} Green function of the $f$ electrons~\eqref{eq:projected-f-gf} exhibits a Kondo resonance at the Fermi energy and the \textit{local} Green function of the $c$ electrons~\eqref{eq:dmft-local-gf-2} a dip, as conduction electrons are removed from the Fermi surface as they hybridise with the $f$ electrons. \subsection{Time-resolved collapse and revival of the heavy quasiparticles} The time evolution of the spectral function of the $f$ electrons in \cref{fig:gd-time-trace} shows in great detail how a THz light pulse induces strong non-equilibrium dynamics and leads to the collapse of Kondo coherence. At first, the THz pulse has not interacted with the system, which remains in thermal equilibrium at a low temperature, as evidenced by the sharp Kondo peak and thermal distribution function. However, as the THz pulse starts interacting with the system, the Kondo peak is immediately affected as the pulse excites the low-energy states near zero frequency. As the intensity of the external THz pulse grows, two main physical mechanisms of action of the pulse on the electronic system can be identified: enhanced hybridisation and correlation decoherence. \paragraph{Enhanced hybridisation} Without coupling to the electromagnetic fields, the width $\Delta$ of the single-particle peak around $\varepsilon^f_0$ is given by $\Delta \sim \pi V_0^2 \rho^c_0(\varepsilon_F)$. The dressing of the bare hybridisation vertex~\eqref{eq:vertex-dressing} by the coupling to the electromagnetic field can considerably change $\Delta$ near the THz pulse maximum $\Delta \sim \pi \del{V_0^2 + g_0^2 \bar{n}_a} \rho^c_0(\varepsilon_F)$, where the last term is obtained by convolving a sharp THz spectral function with $\rho^c_0$. The enhanced hybridisation shifts the problem from the Kondo regime to a mixed-valence regime, where the average occupation of the $f$ level can drop significantly, highlighting the importance of including charge fluctuations (encoded in the dynamics of the $b$ auxiliary particle) in the dynamics. Apart from the direct effect on the spectral content of the fields, a larger hybridisation strength results in a larger Kondo temperature and hence a broader and smaller Kondo peak. The coupling to the electromagnetic-field reservoir has a negligible effect on the single-particle peak's frequency (Lamb) shift and width. \paragraph{Correlation decoherence} Despite DMFT signifying the presence of an \textit{infinite} reservoir of particles, its temperature is only fixed at all times before $t_0$. At any time $T > t_0$, the temperature is an indefinite quantity and is only possible to infer through a generalisation of the fluctuation-dissipation relation (\cref{sec:neqft-interpretability}). While driven out of equilibrium by the incident THz pulse, there is a light-mediated energy transfer from the valence band to the conduction band -- interband heating, creating hot charge carriers. As possible to infer from the distributions functions in \cref{fig:gd-time-trace}, the THz pulse induces a non-thermal hot distribution function which washes out the Kondo coherence that was present in the system before the interaction with the THz pulse. This process induces decoherence of the Kondo correlations and fundamentally differs from the heating encountered experimentally when using long light pulses. A pulse interacting for a long time with a material excites all the sub-bands and phononic degrees of freedom, resulting in a heat-up of the reservoirs, which also provides correlation decoherence through thermal fluctuations. However, the correlation decoherence described here is simply an extreme deviation of the microscopic system from a thermal state, as high-energy states from the THz pulse mix with the low-energy thermal electronic states. Without a coupling to a reservoir, the system would reach a stationary-state effective temperature far above the initial temperature. By coupling to a reservoir at the same temperature as the system at initial times, the system first reaches a quasi-equilibrium state at some high temperature\footnote{The behaviour of the distribution function around zero frequency is in agreement with a hot temperature, but is in itself a weak heuristic, as the definition of temperature or thermal distribution can only re-emerge at long times, after thermalisation.} due to the coherence time of conduction electrons being $\sim \upsilon^{-1}$, before a proper transfer of the excess energy to the environment, with a rate entirely determined by the coupling to the reservoir(s). This can be evidenced by the fast recovery of the single-particle peak at $\varepsilon^f_0$ and a Fermi-Dirac-like distribution well before the Kondo peak is reformed, which takes longer because of two reasons. One, the effective temperature must drop below the Kondo temperature -- through interband cooling, with electronic transitions from the conduction to the valence band -- for the Kondo effect to become dominant. And second, because Kondo coherence -- associated with a sharp resonance in frequency around the Fermi energy -- takes a long time to build up because of the intrinsic low-energy effects associated with Kondo physics. \begin{figure}[!htb] \centering \includesvg[width=\linewidth]{figs/gd} \caption{Time evolution of the spectral $-2 \im G^R_{f_\sigma}(T, \omega)_{\mathrm{W}}$ and distribution~\eqref{eq:fdr} (insets) function components of the local $f$ Green function upon drive by a Gaussian THz pulse. The bottom right plot denotes the intensity profile of the THz pulse, and the vertical lines show the time at which the different plots take place. The second to last plot exhibits oscillations due to not existing enough points in the $\tau$-direction~\eqref{eq:wigner-ville} at late centre-of-mass times $T$ to resolve the sharp spectral features, a region shaded in red. A similar effect is not seen in the first plot because the Green functions are defined for times $T < 0$, and there are enough points to resolve the Green functions properly.} \label{fig:gd-time-trace} \end{figure} \subsubsection{Time-resolved view of the heavy quasiparticles} A solution to the lattice Green function~\eqref{eq:dmft-lattice-gf-k} for all momenta is necessary to obtain a time-resolved view of the heavy quasiparticles. However, the DMFT problem was formulated in a Bethe lattice to avoid the elevated computational cost of separately resolving each $\bm{k}$-point. Since the data is interpreted in a Wigner-Ville transformed basis, which already has some caveats in strong non-equilibrium regimes (\cref{sec:wigner-basis}), to add insult to injury, the lattice Green functions are calculated as \begin{equation} G^{-1}_{c_{\bm{k}\sigma}}(T, \omega)_{\tilde{\mathrm{W}}} \stackrel{!}{=} \sbr{\omega - \varepsilon^c_{\bm{k}} - \Sigma_{c_\sigma}(T, \omega)_{\tilde{\mathrm{W}}}}, \end{equation} where the local self-energy of conduction electrons is calculated through \begin{equation} T_{c_\sigma}(T, \omega)_{\tilde{\mathrm{W}}} \del{\lim_{\zeta\to0} G_{c_{\zeta,\sigma}}(T, \omega)_{\tilde{\mathrm{W}}}} \stackrel{!}{=} \Sigma_{c_\sigma}(T, \omega)_{\tilde{\mathrm{W}}} G_{c_\sigma}(T, \omega)_{\tilde{\mathrm{W}}}~. \end{equation} The self-energy $\Sigma_{c_\sigma}(T, \omega)_{\tilde{\mathrm{W}}}$ is, of course, not the \textit{true} local self-energy of the system~\eqref{eq:dmft-self-energy}. However, for qualitative assessment, it should contain the most important features, especially since when Kondo coherence collapses, the system loses its low-energy features and the timescales of the pulse drive are much slower than the intrinsic timescales of the system (\cref{sec:separation-timescales}). Away from transient effects, it is expected\footnote{ Consider the rotation $G_{c_{\bm{k}\sigma}}(t, t') = e^{-i \varepsilon^c_{\bm{k}} (t - t')} \bar{G}_{c_{\bm{k}\sigma}}(t, t')$, where~\eqref{eq:dmft-lattice-gf-k} reads \begin{equation*} i \partial_t \bar{G}_{c_{\bm{k}\sigma}}(t, t') = \int_\gamma \dif\bar{t}\, \Sigma_{c_\sigma}(t, \bar{t}) e^{+i \varepsilon^c_{\bm{k}} (t - \bar{t})} \bar{G}_{c_{\bm{k}\sigma}}(\bar{t}, t')~. \end{equation*} The Bogolyubov principle states that the temporal correlations in a system decay within a period with characteristic time $\Lambda$. This is translated directly into 2-point functions such as $G_{c_{\bm{k}\sigma}}(t, t')$ or self-energies $\Sigma_{c_\sigma}(t, t')$, which should have negligible values for outside the strip $\envert{t - t'} \lesssim \Lambda$. Hence, \begin{equation*} i \partial_t \bar{G}_{c_{\bm{k}}}(t, t') \approx \int_\gamma \dif \bar{t}\, \bar{\Sigma}(t, \bar{t}) e^{-\frac{(t-\bar{t})^2}{2 \Lambda^2}} e^{+i \varepsilon_{\bm{k}} (t - \bar{t}) } \bar{G}_{\bm{k}}(\bar{t}, t') \approx \bar{\Sigma}(t, t) \bar{G}_{c_{\bm{k}}}(t', t') e^{-\frac{1}{2}\varepsilon_{\bm{k}}^2 \Lambda^2} \sqrt{2\pi \Lambda^2} \end{equation*} Given that $\Lambda \sim \frac{\hbar}{k_\mathrm{B} T^*_K} \gg 1$, the conduction electrons that get renormalised are the ones with small kinetic energy, namely the ones with energy close to $k_\mathrm{B} T^*_K$.} that the self-energy renormalises the momenta close to the Fermi energy strongly. For vanishing light-matter coupling, the lattice Green function of $f$ electrons can be calculated as \begin{equation} G_{f_{\bm{k}\sigma}}(T, \omega)_{\tilde{\mathrm{W}}} \stackrel{!}{=} \frac{1}{V_0^2}\sbr{\Sigma_{c_\sigma}(T, \omega)_{\tilde{\mathrm{W}}} + \Sigma_{c_\sigma}(T, \omega)_{\tilde{\mathrm{W}}} G_{c_{\bm{k}\sigma}}(T, \omega)_{\tilde{\mathrm{W}}}\Sigma_{c_\sigma}(T, \omega)_{\tilde{\mathrm{W}}}}~. \end{equation} In \cref{fig:kondo-collapse-k} the collapse and revival of the heavy quasiparticles can be seen in detail. Before the interaction of the THz pulse, the system exhibits heavy quasiparticles, identified by the flat and thin spectral intensity near zero frequency $\omega$ for most values of $\kappa$. Despite NCA + DMFT being unable to capture a Kondo insulating phase ($J_K \gtrsim 1.0$) without introducing doubly-occupied $f$-level states, the heavy quasiparticles form a many-body indirect gap near $\omega = 0$, bordered by two peaks of width $\sim T^*_K$. Furthermore, a single-particle hybridisation gap near $\omega = \varepsilon^f_0$, also arising from the hybridisation of $f$ and $c$ particles, is also visible. The THz pulse induces a momentary mixed-valence regime and subsequent loss of the large effective electron mass as the hybridisation and indirect gap smear out entirely, and the latter shifts towards the Fermi energy. This signals the collapse of the heavy quasiparticles as Kondo coherence melts, the large Fermi surface shrinks and the system is characterised by a higher metallic character. Upon the disappearance of the THz pulse, the system starts its thermalisation process. The single-particle hybridisation gap recovers first, as it is a single-particle effect with an associated fast timescale $\sim V_0^2 \rho^0_c(\varepsilon_F)$. Only long after the interaction of the THz pulse is the Kondo effect recovered, and the heavy quasiparticles reformed again around the edges of the indirect gap. \begin{figure}[!htb] \centering \includesvg[width=\linewidth]{figs/collapse-revival-heatmap} \caption{Time-resolved collapse and revival of the heavy quasiparticles. In the top panel, the spectral function of the charge carriers $-2 \im \sbr{G_{c_{\bm{k}}}(T, \omega)_{\tilde{\mathrm{W}}} + G_{f_{\bm{k}}}(T, \omega)_{\tilde{\mathrm{W}}}}$, and in the bottom panel the spectral function of the $f$ electrons $-2 \im G_{f_{\bm{k}}}(T, \omega)_{\tilde{\mathrm{W}}}$, where $\varepsilon^c_{\bm{k}} \coloneqq \kappa$, are plotted for different centre-of-mass times, denoted in \cref{fig:gd-time-trace}.} \label{fig:kondo-collapse-k} \end{figure} \subsection{Photon re-emission} In a reflection geometry, such as in the THz experimental setting of~\cite{Wetli_2018}, an incident THz pulse arrives, interacts with the system, is reflected, and propagates away, carrying some information about the system (\cref{fig:time-reflectivity}). There is an \textit{instantaneous} reflection originating from stimulated intra-band excitations within the conduction bands, which leave the heavy quasiparticles intact -- and are even decoupled from the model~\eqref{eq:full-model}, which only contains the conduction electrons which interact with the $f$ electrons. Additionally, a \textit{delayed} -- echo-like -- reflection response is observed, originating from inter-band transitions between hybridised conduction and $f$ electrons and understood as stemming from the recovery of the Kondo singlets after their destruction by the THz pulse. Establishing a faithful correspondence between theory and such experiments is somewhat complicated: unlike the latter, which can only measure reflection and not incidence intensity, the former can only compute local temporal \textit{distortions} of the incident intensity. Geometrical and dispersive aspects of the THz pulse cannot be fully considered given the local description of light interacting with matter, namely the light-matter Hamiltonian~\eqref{eq:full-model} and the associated equations of motion of the non-interacting part of pulse~\eqref{eq:photon-eom}. Nonetheless, despite resulting in a small distortion due to the smallness of $g_0$, the intensity of the renormalised incident pulse shows interesting dependence on the system parameters, corroborating the experimental observations. For short times $T \ll 1/{g_0}$ upon the pulse interacting with the system, the incident pulse mode is not distorted, as spontaneous emission either into the pulse mode or into the environment is negligible. At later times, emission by one-photon transitions drives the system back to its original state through radiative recombination, renormalising the temporal mode of the incident pulse, which develops a macroscopic \textit{delayed} secondary pulse. This long sought-after experimental signature encodes and carries information about the ground state and possibly physics hidden within the underlying system (\cref{sec:thz}). The dependence of the delayed pulse on several system parameters strongly indicates its relation with low-energy Kondo physics, decisively explaining in great part the time-resolved THz experimental observations of~\cite{Wetli_2018}. In the following analyses, the renormalised photon intensity -- obtained through the lesser component of~\eqref{eq:nca-a-projected} -- is shifted by its background value $G^<_a(T, \tau)_\mathrm{W} \to G^<_a(T, \tau=0)_\mathrm{W} - \lim_{\Gamma \to \infty} G^<_a(\Gamma, \tau=0)_\mathrm{W}$. Since the electronic system is always coupled to the electromagnetic field, the renormalisation of the latter always results in a finite but experimentally unmeasurable background photon number. Similarly, the normalised intensity of the delayed pulse of the inset plots is obtained via $-\im \frac{G^<_a(T, \tau=0)_\mathrm{W}}{\lim_{\Gamma \to \infty} G^<_a(\Gamma, \tau=0)_\mathrm{W}} - 1$, and measures how much brighter the delayed pulse is in comparison with the background noise. The incident pulse intensity $-\im G^<_{a_0}(T, \tau=0)_\mathrm{W}$, shaded in grey, is visible in all plots, but due to the large zoom of the vertical axis, its Gaussian profile~\eqref{eq:ga-lesser} is not visible. Finally, a red-shaded area indicates where some results may become inaccurate due to the lack of points in the $\tau$ direction to resolve sharp spectral features. \paragraph{Delayed pulse vs pulse intensity (\cref{fig:reflectivity-5})} Experiments with light interacting with matter are typically performed in two different regimes: the low-fluence regime, where the system is perturbed as gently as possible to minimise heating effects and the high-fluence regime, where the system is perturbed as strongly as required, possibly non-perturbatively, to induce phase transitions or create non-thermally-accessible metastable states~\cite{Basov_2011}. Experiments such as~\cite{Wetli_2018} are within the former regime, with a low photon flux hitting the heavy-fermion compounds. However, there is a long-standing question on whether the observed experimental delayed pulse is not related to Kondo physics but instead to superradiant decay. The archetype example of superradiance is when a dense ensemble of incoherently excited two-level systems lock their dipoles in phase and develop a macroscopic dipole proportional to the number of inverted atoms $N$~\cite{Cong2016}. The dipole decays at an accelerated rate, emitting a delayed pulse with intensity $I(t) \approx \gamma \del{\frac{N}{2}}^2 \operatorname{sech}^2 \sbr{\gamma \frac{N}{2} (t - t_\mathrm{D})}$, for some decay rate $\gamma$ and delay time $t_\mathrm{D}$~\cite{Benedict2018}. The intensity of the emitted light is characterised by being proportional to the \textit{square} of the number of excited atoms $N$ and its duration inversely proportional to $N$. These characteristics are in strong opposition with the \textit{linear} dependence of the delayed pulse on the intensity of the incident pulse, controlled by $\bar{n}_a$. This parameter should also increase the number of excited electrons in the system, even though a heavy-fermion lattice is undoubtedly more complicated than an ensemble of weakly interacting two-level systems. An estimate of the \textit{fraction} of excited atoms is given by \begin{equation} N \approx \max_t \sbr{G^>_{c_\sigma}(t, t) + G^>_{f_\sigma}(t, t)}~, \end{equation} which measures the charge carriers' maximum depopulation (or hole number) upon interaction with the external pulse. Here, the delayed pulse's intensity also depends linearly on $N$. Finally, the duration of the delayed pulse does not depend on the inverse of $N$, another definitive indication that the obtained delayed pulse is fundamentally different from superradiant emission. \begin{figure}[!htb] \centering \includesvg[width=\figwidth]{figs/reflectivity-5} \caption{Time trace of the renormalised photon pulse as a function of the incident pulse's maximum intensity $\bar{n}_a$ and $f$ electron energy $\varepsilon^0_f$. The insets show a linear dependence of the maximum intensity of the delayed pulse on the incidence intensity $\bar{n}_a$ and on the fraction of excited electrons $N$.} \label{fig:reflectivity-5} \end{figure} \paragraph{Delayed pulse vs pulse frequency (\cref{fig:reflectivity-3})} Another key dependence of the delayed pulse is on the incident pulse central frequency $\omega_0$. For a large frequency $\omega_0$, the pulse distortion appears to be a longer -- or more extended -- exponential-like decay rather than a secondary pulse. Upon lowering $\omega_0$ to a physical regime\footnote{Despite being possible to generate results for larger $\omega_0$, the simplifying assumptions on the dipole matrix element (\cref{sec:alm-light}), such as zero-momentum transfer on light-mediated electronic transitions, become unfounded. For a conduction bandwidth $D \sim 4 \upsilon \sim 4$ eV, a photon frequency of $\omega_0 = 0.1 \upsilon \sim 25$ THz is already roughly $100$ times larger than in the experiments.} more closely related to the experimental settings~\cite{Wetli_2018}, the delayed pulse, with a slightly longer decay time, is recovered. The reason why this happens is encoded in the inset of \cref{fig:reflectivity-3}, which shows how the retarded component of the t-matrix $T_a$ of the photon pulse looks like in thermal equilibrium (or, equivalently for the purpose of this analysis, when Kondo coherence is fully developed). The t-matrix $T_a$ comprises a fast decaying rapid oscillation (related to the broad single-particle peak centred at $\varepsilon^0_f$) and a slow decaying slow oscillation (associated with the sharp Kondo peak centred almost at $0$). The renormalisation of the photon pulse, obtained through time convolutions~\eqref{eq:nca-a-projected}, results in vastly different outcomes depending on the oscillations of the photon Green functions. For a large central frequency $\omega_0$, the slow decaying part of $T_a$ roughly averages to zero and hence temporal information about Kondo coherence is lost, and only single-particle information is retained. Conversely, for a small $\omega_0$, the slow decaying part of $T_a$ is retained, resulting in an echo-like distortion of the incident temporal pulse mode. This shows how Kondo coherence is required for a secondary echo-like pulse to form, as opposed to an exponential-like decay, characteristic of single-particle relaxations. \begin{figure}[!htb] \centering \includesvg[width=\figwidth]{figs/reflectivity-3} \caption{Time trace of the renormalised photon pulse as a function of the incident pulse's central frequency $\omega_0$. The inset shows the time trace in the relative time $\tau$ direction of the retarded components of the incident photon's Green functions and (rescaled) t-matrix $T_a$ in thermal equilibrium, hence centre-of-mass time $T$ independent.} \label{fig:reflectivity-3} \end{figure} % \begin{figure}[!htb] % \centering % \includesvg[width=\figwidth]{figs/reflectivity-3} % \caption{Time trace of the renormalised photon pulse as a function of the incident pulse's central frequency $\omega_0$.} % \label{fig:reflectivity-3} % \end{figure} % \begin{figure}[!htb] % \centering % \includesvg[width=\figwidth]{figs/t-matrix} % \caption{Time trace in the relative-time $\tau$ direction of the incident photon's Green functions and t-matrix retarded components in thermal equilibrium, hence centre-of-mass time $T$ independent.} % \label{fig:reflectivity-3-2} % \end{figure} % \paragraph{Delayed pulse vs pulse bandwidth} \cref{fig:reflectivity-4} Despite a system described by smaller $J_K$ being more susceptible to photoexcitation. % \begin{figure}[!htb] % \centering % \includesvg[width=\figwidth]{figs/reflectivity-4} % \caption{Dependence of the delayed pulse on the pulse bandwidth.} % \label{fig:reflectivity-4} % \end{figure} \paragraph{Delayed pulse vs bath temperature (\cref{fig:reflectivity-1})} Ultimately, the delayed pulse is strongly temperature dependent, vanishing at high temperatures and growing as the bath temperature is lowered -- both in qualitative agreement with experiments~\cite{Wetli_2018}. At high temperatures, the electronic system is described by single-particle physics, with the system being composed of localised $f$ particles interacting weakly with delocalised conduction electrons. Here, there are no heavy quasiparticles since thermal fluctuations wash out many-body quantum effects, and light-induced effects have a negligible impact on the system. However, as the bath temperature approaches and is set below $T^*_K$, there is a pronounced renormalisation of the incident pulse, resulting in a secondary, \textit{delayed} pulse. The amplification of the delayed pulse's intensity with the lowering of the bath temperature is another indication of its relation with low-energy physics, appearing to plateau as $\beta\to\infty$, as the Kondo effect saturates (much lower temperatures cannot be reached numerically within DMFT + NCA and these Anderson lattice model parameters). \begin{figure}[!htb] \centering \includesvg[width=\figwidth]{figs/reflectivity-1} \caption{Time trace of the renormalised photon pulse as a function of the bath temperature $\beta$. The inset shows a logarithmic dependence of the maximum intensity of the delayed pulse on $\beta$. The Kondo timescale $\tau^*_K$ is defined as $\tau^*_K = \frac{\hbar}{k_\mathrm{B} T^*_K}$.} \label{fig:reflectivity-1} \end{figure} \paragraph{Delayed pulse vs bath coupling (\cref{fig:reflectivity-2})} A larger bath coupling strength $\alpha$ induces more decoherence in conduction electrons, which washes out some of the Kondo coherence, directly influencing the delayed pulse, which becomes less intense and with a shorter tail. Furthermore, $\alpha$ sets the timescale for thermalisation, controlling how long it takes for the system to build back its in-equilibrium coherence, estimated by the Kondo peak relative height. A system with a lower $T^*_K$, related to a larger $\varepsilon_0^f$, is more stable to the interaction with the incident pulse, noting a considerably weaker collapse of the Kondo peak height and a less intense delayed pulse. This indicates that the delayed pulse is related not only to the Kondo coherence but also to how much of it is destroyed. Despite the height roughly measuring how much Kondo coherence there is, it never vanishes entirely (cf.~\cref{fig:gd-time-trace}) as it is just the value of the spectral function at $\omega = 0$ and is not an isolated measurement of the spectral weight of Kondo quasiparticles. \begin{figure}[!htb] \centering \includesvg[width=\figwidth]{figs/reflectivity-2} \caption{Time trace of the renormalised photon pulse as a function of the bath coupling strength $\alpha$ and $f$ electron energy $\varepsilon^0_f$. The time trace of the relative Kondo peak height (unity when the system is in thermal equilibrium with the bath) is juxtaposed. Its drop-off at large times $T$ is related to the lack of points in the relative-time direction to properly resolve the sharp spectral features of Kondo observables (cf.~\cref{fig:gd-time-trace}).} \label{fig:reflectivity-2} \end{figure}