\chapter{Heavy-fermion systems driven by a terahertz light pulse} \label{sec:thz} Heavy-fermion compounds boast an intricate phase diagram (\cref{fig:qcp}) due to the interplay of magnetic- and electric-ordering quantum fluctuations arising from a high concentration of localised valence electrons interacting with mobile conduction electrons. For example, in the heavy-fermion compound CeCu$_{6-x}$Au$_{x}$, there is a competition controlled by chemical doping between a heavy Fermi liquid and an anti-ferromagnetic phase, with a quantum critical point at $x=0.1$ doping fraction~\cite{Schr_der_2000}. The CeCu$_{6}$ compound is in a heavy Fermi liquid phase, and doping with Au atoms distorts the lattice and induces a reduction of the exchange coupling strength between the $f$ and conduction electrons responsible for the Kondo effect, which tends to \textit{locally} screen of the spin of valence $f$-electrons and is behind the heavy Fermi liquid phase (\cref{sec:kondo-pheno}). The lattice distortion favours the \textit{long-range} Ruderman–Kittel–Kasuya–Yosida (RKKY) interaction between neighbouring $f$-electrons, which tends to align anti-ferromagnetically the valence electrons in the lattice. Because of these quantum phases' inherent complexity and antagonist character, there is still an ongoing debate about the system's quasiparticles near or at the quantum critical point, where the phases can co-exist. In the work~\cite{Wetli_2018}, CeCu$_{6-x}$Au$_{x}$ was irradiated by an ultrafast\footnote{The ultrafast timescales are between the femtosecond ($10^{-15}$ s) and the nanosecond ($10^{-9}$ s), with one picosecond ($10^{-12}$ s) being to one second as one second is to approximately 32000 years.} pulse of terahertz (THz) radiation. THz time-domain spectroscopy is one of the few experimental techniques that can temporally resolve the ultrafast dynamics of materials, which can be used to inspect and characterise their low-energy excitations. Unlike photo-emission experiments, which use intense ultraviolet radiation to eject photo-excited electrons, ultrafast THz spectroscopy only slightly stirs up the electrons with non-ionising radiation and evades the pitfalls of the former such as lattice distortion and heating, which would heavily disturb the delicate low-energy physics of these compounds. A pulse of linearly polarised THz light with a frequency range of $0.1$–$3$ THz photoexcites a CeCu$_{6-x}$Au$_{x}$ sample, which initiates a dynamical response. Due to the penetration depth of THz fields in conductive mediums being only a few nanometers long, the reflected (instead of the transmitted) electric field is measured. For the anti-ferromagnetic compound CeCu$_{5}$Au$_{1}$, as well as for a test Pt mirror~\cite{Wetli_2017}, only an instantaneous reflection of the pulse is observed (\cref{fig:time-reflectivity}). Such immediate response is, however, expected for any metallic compound and is simply the reflectivity caused by conduction electrons. However, for the heavy Fermi liquid CeCu$_{6}$, and the quantum critical CeCu$_{5.9}$Au$_{0.1}$ compounds, a delayed reflected pulse -- or echo -- appears (\cref{fig:time-reflectivity}). \begin{figure}[!htb] \centering \includegraphics[width=0.8\figwidth]{diagrams/sketches/phase-diagram} \caption{Phase diagram of heavy-fermion systems. } \label{fig:qcp} \end{figure} The delayed pulse does not show an exponential decay, as expected for weakly interacting or single-particle relaxation. Instead, it is a compact pulse, with a pronounced \enquote{dark time} of $\tau_\mathrm{echo} \approx 6.2$ ps in CeCu$_{6}$ and $\tau_\mathrm{echo} \approx 5.8$ ps in CeCu$_{5.9}$Au$_{0.1}$ after the instantaneous reflection. Furthermore, coherence time in metals is typically within femtosecond timescales~\cite{Knoesel_1998}, which hints at the relation of the echo with low-energy excitations with longer characteristic timescales. Namely, the delay time agrees\footnote{\label{fn:time-energy-uncer} The finite lifetime of Kondo quasiparticles can be related to the width of their excitation spectrum $k_\mathrm{B} T_K$ (\cref{sec:kondo-pheno}) via the time-energy uncertainty relation~\cite{Briggs_2008}. Its nature is very different from Heisenberg's uncertainty relations due to time not being an operator. However, despite its correctness being somewhat disputed, it is of great heuristic value. In this case, the uncertainty in the energy of Kondo quasiparticles leads to an uncertainty in the lifetime of its quasiparticles, related to the time for photon emission.} with the reported~\cite{Klein_2008} compound's Kondo lattice-coherence temperature $\tau_\mathrm{echo} \sim \frac{h}{k_\mathrm{B} T^*_K} \approx 8$ ps, where $h$ is the Planck's constant and $k_\mathrm{B}$ the Boltzmann constant, and is also strongly temperature dependent, vanishing at high temperatures, also in remarkable agreement with Kondo physics. Hence, a direct link was established between the origin of the echo and the Kondo quasiparticles characteristic of the heavy Fermi liquid phase of the heavy-fermion compound. The generation of the echo pulse was qualitatively described via a rate equation~\cite{Wetli_2018} which proposed it to be a response from the reconstruction of the Kondo lattice coherence after its destruction by the incident pulse. The THz radiation induces dipole interband transitions between the heavy $4f$ and the light $5d$ bands, which break the Kondo singlets. Due to an occupation-dependent transition rate specific to heavy Fermi liquid systems, the electrons' recombination can only occur after enough time has passed for the build-up of the Kondo lattice coherence, leading to a delayed response. The importance of this class of experiments lies in directly addressing questions such as whether $T^*_K$ vanishes or remains finite at the quantum critical point of CeCu$_{6-x}$Au$_x$~\cite{Schr_der_2000} and hence whether heavy quasiparticles survive or not criticality. \begin{figure}[!htb] \centering \includegraphics[width=\figwidth]{figs/wetli-time-reflectivity.pdf} \caption{Time-resolved reflectivity of CeCu$_{6-x}$Au$_{x}$, colour matched with the phases of~\cref{fig:qcp}. The instantaneous reflection occurs at sampling time $0$, and the time-delayed reflected pulse is highlighted in the insets. Reproduced from~\cite{Wetli_2018}, with permission.} \label{fig:time-reflectivity} \end{figure} \section{The Anderson lattice model} \label{sec:anderson-lattice-model} %Mixed-valence Ce compounds such as CeCu$_6$ are characterised by having valence electrons in the $4f$ shell of Ce. The seven orbitals of the $4f$ shell -- each doubly degenerate due to the electrons' spin -- are split in energy by spin-orbit coupling and crystal field effects and, when embedded in a lattice, create narrow valence bands due to the spatial confinement of the orbitals, with only a few bands left at the vicinity of the Fermi level in CeCu$_6$. In accordance with the valence of Ce in CeCu$_6$~\cite{Pal2019}, the $4f$ shell is considered to contain just one valence electron and, due to strong Coulomb repulsion -- weakly screened on atomic length scales -- only the lowest-lying doublet (\cref{fig:mutiplets}) is significantly occupied, with the physics stemming from considering the other doublets deemed unimportant at low temperatures. \begin{figure}[!htb] \centering \includegraphics[width=0.5\figwidth]{diagrams/sketches/multiplets} \caption{Schematic of $f$-orbital splitting in CeCu$_6$~\cite{Matsumoto_2020}. Spin-orbit coupling splits the $4f$ orbitals of Ce in $j=5/2$ and $j=7/2$ multiplets, separated by $\Delta \approx 250$ meV~\cite{Ehm_2007}. Crystal field effects further split the $j=5/2$ multiplet into three doublets separated by $\delta_1=7$ meV and $\delta_2=15$ meV.} \label{fig:mutiplets} \end{figure} The Coulomb interaction strength $U$ in $f$-valence compounds is the largest energy scale, typically comparable with the conduction bandwidth~\cite{Gunnarsson_1985}. As a result, the empty $f^0$ and singly-occupied $f^1$ state of the doublet are close in energy and much lower than the doubly-occupied $f^2$, which can be \textit{projected out} in low-energy models. However, this is not universal, and for systems with larger energy separations between the $f^0$ and $f^1$ configuration, a $U\to\infty$ limit may not even be suitable for qualitative studies. Although initially introduced to describe the effects of a low concentration of magnetic impurities in metals, it is believed that the $U\to\infty$ limit of the Anderson lattice model~\cite{Anderson_1961} captures the essential low-energy physics of heavy-fermion compounds~\cite{Millis1987} \begin{equation} \label{eq:ALM} \hat H_{\mathrm{ALM}} = -\sum_{\ev{i,j} \sigma} t^c_{ij}\hat c^\dagger_{i\sigma}\hat c_{j\sigma} -\sum_{\ev{i,j} \sigma} t^f_{ij}\ket{\sigma}^f_i \bra{\sigma}^f_j + \sum_{i \sigma} \varepsilon_0^f \ket{\sigma}^f_i \bra{\sigma}^f_i + V_0 \sum_{i \sigma} \del{\hat c^\dagger_{i \sigma} \ket{0}^f_i\bra{\sigma}^f_i + \hc}~. \end{equation} This model describes the hybridisation of a band of non-interacting conduction electrons $c$ (related to the extended $s$, $p$ or $d$ orbitals of the compound) with localised interacting electrons $f$. The $f$ electrons are considered \textit{localised} due $\varepsilon^f_0 < 0$ and the spatial confinement of the $f$ orbitals, resulting in a hopping strength $|t^f| \ll |t^c|$. The Anderson lattice model parameters of CeCu$_6$ were estimated~\cite{Matsumoto_2020} to be $\varepsilon^f_0 = -1.61$ eV, $V_0 = 0.41$ eV and $U \sim 5$ eV, similar to the conduction bandwidth, which supports the doubly-occupied $f$-states to be projected out of the theory. The states $\ket{0}^f_i$ and $\ket{\sigma}^f_i$ represent, respectively, the empty $f^0$ and singly occupied $f^1$ configurations of an $f$-electron with spin $\sigma$ on site the $i$. The Hamiltonian is equivalently described in terms of creation/annihilation operators $\hat f$ acting on a \textit{restricted} Hilbert space, subject to the operator constraint \begin{equation} \sum_{\sigma} \hat f^\dagger_{i\sigma} \hat f_{i\sigma} \le 1~. \end{equation} This inequality, challenging to implement due to the $c$-$f$ hybridisation, can be turned into an equality by the addition of a new auxiliary particle (\cref{sec:aux-particles}), a boson field $b$, with the projection operators of the Hamiltonian mapped to \begin{equation} \ket{0}^f_i\bra{\sigma}^f_i \to \hat b^\dagger_i \hat f_{\sigma i}~, \qquad \ket{\sigma}^f_i\bra{\sigma}^f_i \to \hat f^\dagger_{i \sigma} \hat f_{i \sigma}~. \end{equation} Due to projection operators, the strong Coulomb repulsion character was built-in from the beginning, and thanks to the auxiliary-particle mapping, all particle operators satisfy canonical commutation relations and the hybridisation is promoted to a proper interaction, amenable to perturbative methods. The resulting Hamiltonian reads \begin{equation} \label{eq:ALM-ap} H_{\mathrm{ALM}} = -\sum_{\ev{i,j} \sigma} t^c_{ij}\hat c^\dagger_{i\sigma}\hat c_{j\sigma} -\sum_{\ev{i,j} \sigma} 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} + V_0 \sum_{i \sigma} \del{\hat c^\dagger_{i \sigma} \hat b^\dagger_i \hat f_{i \sigma} + \hc}~, \end{equation} where the operator constraint \begin{equation} \label{eq:q-constraint} \hat{Q}_i = \hat b_i^\dagger \hat b_i + \sum_\sigma \hat f^\dagger_{i\sigma} \hat f_{i\sigma} = \hat{\mathds{1}}~, \end{equation} must be enforced at all lattice sites $i$. \subsection{Phenomenology of Anderson-impurity models} \label{sec:kondo-pheno} For a \textit{single}-impurity Anderson model, where there is just one $f$ shell (denoted as the \enquote{impurity}) in a sea of conduction $c$ electrons, the model is characterised by the Kondo $|\varepsilon^f_0| \gg \Delta$ and the mixed-valence $\Delta \gtrsim |\varepsilon^f_0|$ regimes. In both regimes, the physics is dominated by the Coulomb repulsion, whose strength $U$ is much larger than the hybridisation strength $W_0 = V_0^2 \rho^c_0(\varepsilon_F)$, where $\rho^c_0(\varepsilon_F)$ is the bare density of states of $c$-electrons at the Fermi energy. But while in the former, the $f$-shell average occupation is one, in the latter, the proximity of the $f$-shell bare energy to the Fermi energy leads to a mixed-valence regime, with an average occupation of less than one (i.e., possible loss of the \enquote{local moment}). \subsubsection{The Kondo effect} Focusing on the Kondo regime, the strong Coulomb repulsion favours the $f$ shell to be occupied by either a spin-up or spin-down electron, and a \textit{sequential} tunnelling process between $c$ and $f$ electrons is energetically preferred. The Kondo temperature $T_K$ \begin{equation} \label{eq:kondo-temp} T_K = D \del{\frac{N W_0}{D}}^{1/N} \exp\del{\frac{\varepsilon^f_0}{N W_0}}~, \end{equation} where $D$ is the half-bandwidth of $c$ electrons and $N$ is the degeneracy of the local moment~\cite{Kim_1990}, sets the energy scale for which the tunneling becomes \textit{resonant}. As the temperature is lowered $T \lesssim T_K$, scattered electrons scatter again coherently with the impurity and become \textit{correlated}. Quantum mechanically, an electron can hop out of the impurity, within a short timescale $\hbar / |\varepsilon_0^f|$ (\cref{fn:time-energy-uncer}) into the surface of the Fermi sphere, which populates the impurity level with another electron, possibly with a different spin~\cite{Kouwenhoven_2001}. The sequential spin-flip scattering between impurity electrons and excitations in the Fermi sea gives rise to a \textit{many-body} resonance, which exhibits a sharp peak of width $\sim k_\mathrm{B} T_K$ in the impurity electrons' spectral function at the Fermi energy. The resonance is known as the Abrikosov-Suhl or Kondo resonance and is the smoking gun of the Kondo effect. % As the temperature approaches zero, the spin fluctuations are frozen, and a cloud of conduction electrons completely screens the $f$-shell electron. At high temperatures, $k_\mathrm{B} T \gg U$, decoherence through thermal fluctuations washes out many-body coherent effects, and the physics is dominated by single-particle dynamics. In a lattice of $f$ and $c$- electrons, resonant scattering at each lattice site will generate a Kondo lattice-coherent resonance with width\footnote{Unlike the single-impurity Kondo temperature $k_\mathrm{B} T_K$, which is understood as the energy at which perturbative renormalisation group breaks down~\cite{Coleman_2015_hf}, there is not an agreed formula for the Kondo (lattice-)coherence temperature $T^*_K$.} $\sim k_\mathrm{B} T^*_K$ at the Fermi energy of the on-site spectral function of each $f$-electron. As these are embedded in a lattice, a coherent $f$-band is formed, with significant spectral density at the Fermi energy because of the resonance. This increases the number of mobile electrons in the system, and due to the local character of $f$-electrons, the band is mostly flat, resulting in a large effective mass of the charge carriers at low temperatures. These \textit{heavy} itinerant electrons lead to an observable expansion of the Fermi volume, as the Fermi surface must expand to accommodate the additional number of indistinguishable electrons in the Fermi sea. This is, however, not general, as reported in~\cite{Kummer_2015}. Despite the RKKY interaction not being addressed in this thesis, the relation of its energy scale $E_\mathrm{RKKY} \sim J_K^2 \rho_0^c(\varepsilon_F)$~\cite{Coleman_2015_hf} with the Kondo temperature can roughly estimate the phase of the heavy-fermion compound. For $k_\mathrm{B} T^*_K \gg E_\mathrm{RKKY}$ the system is in a heavy Fermi liquid phase and decreasing the exchange coupling $J_K \sim -2 V^2\frac{U}{(U + \varepsilon_0^f)\varepsilon_0^f}$~\cite{Frithjof_2012} between $c$ and $f$ electrons can lead to $E_\mathrm{RKKY} \gg k_\mathrm{B} T^*_K$ and hence an anti-ferromagnetic phase. \section{Light-matter interaction} The literature has typically treated light-matter interactions in non-equilibrium many-body systems via the Peierls substitution, which dresses the electronic hopping matrix elements $t_{ij}$ as \begin{equation} t_{ij} \to t_{ij} e^{i \int_{\bm{R}_i}^{\bm{R}_j} \dif \bm{x} \cdot \bm{A}\del{\bm{x}}}~, \end{equation} where $\bm{A}$ is the electromagnetic (EM) vector potential. Although a compact and elegant formulation -- and even finding some uses in strong light-matter coupling scenarios, this formulation has some limitations. First, the light-matter coupling strength is fixed by the strength of the field and does not depend on material properties~\cite{Dmytruk_2021}. Second, it only couples light to non-local hopping elements and neglects local hoppings, such as light-induced orbital transitions. Moreover, $\bm{A}$ is typically considered a classical field, resulting in a semi-classical theory of light-matter interactions. Even though such semi-classical approximations can be appropriate for the study of some non-linear optics or laser theory~\cite{Loudon2000}, a \textit{quantum} theory of light is required for classes of problem dealing with spontaneous/stimulated emission or coherent emission, such as super-radiance -- with similar quantum effects expected from the interaction of an external pulse of THz radiation. \subsection{Quantisation of the electromagnetic field} % The Maxwell action of \textit{classical} electromagnetism is given by % \begin{equation} % S_\mathrm{Maxwell} = -\frac{1}{4} \int \dif^4 x\,F_{\mu \nu}F^{\mu \nu} - \int \dif^4x\, A_\mu J^\mu \qquad F_{\mu \nu} = \partial_\mu A_\nu - \partial_\nu A_\mu~, % \end{equation} % where $J^\mu$ is some external current depending on fields different than $A_\mu$, and the sum is considered to follow the Minkowski space sign convention $(+\,-\,-\,-)$. The Euler-Lagrange equations of motion % \begin{equation} % \fd{}{A_\nu}S_\mathrm{Maxwell} = \partial_\mu F^{\mu\nu} + J^\nu = 0 % \end{equation} % correspond to the inhomogeneous Maxwell equations, whereas the homogeneous ones follow directly from the Bianchi identity % \begin{equation} % \partial_\gamma F_{\mu \nu} + \partial_\mu F_{\nu \gamma} + \partial_\nu F_{\gamma \mu} = 0~, % \end{equation} % with % \begin{equation} % \frac{1}{c}E^{i} = -F^{0 i} % \qquad % \epsilon^{ijk}B_k = -F^{ij}~, % \end{equation} % where the Latin indices denote space indices. The theory described the fields $A_\mu$ and action $S_\mathrm{Maxwell}$ is invariant under the transformation % \begin{equation} % A_\mu \to A_\mu + \delta_{\mu} f~, % \end{equation} % which is known as \textit{gauge} invariance. Unlike, e.g., global symmetries, gauge symmetries are not physical and are simply a redundancy of the formalism chosen to describe some theory~\cite{Schwartz_2013}. By \textit{choosing} an appropriate gauge, this invariance of the theory can be exploited to simplify the equations for the fields $E$ and $B$ and the canonical quantisation of the theory. The Coulomb gauge % \begin{equation} % \bm{\nabla} \cdot \bm{A}(\bm{x},t) = 0 % \end{equation} % is chosen, i.e., if $A^\mu$ does not fulfil the Coulomb gauge, a gauge transformation with $\nabla^2 f = - \bm{\nabla}\cdot\bm{A}$ ensures the Coulomb gauge always holds. % The main hurdle in quantising the EM field is finding the canonical conjugate of the field $A_\mu$. Namely, the conjugate fields of $A_\mu$ % \begin{equation} % \Pi_\mu \coloneqq \pd{L}{(\partial_0 A^\mu)} = -F^{0 \mu} % \end{equation} The canonical \textit{quantisation} of a classical theory promotes the classical Poisson brackets, related to the canonical variables, to a commutator\footnote{In Maxwell's theory, the canonical variable is the vector potential $\bm{A}$ and its conjugate variable is the transverse electric field $E^\perp \coloneqq - \pd{\bm{A}}{t}$. Such relations closely resemble the momentum being the time-derivative of the position (up to a constant) of classical mechanics.}. Regrettably, the condensed-matter literature awkwardly overlooks many important theoretical details of this procedure when quantising the EM field, such as the emergence of the Coulomb interaction from coupling the EM field to matter currents. The actual quantisation is also typically framed as an inexplicable replacement of scalar coefficients by quantum operators. Running the risk of falling down the rabbit hole on a by-now textbook problem, refer to~\cite{Weinberg_1995_ed} for a complete and rigorous quantisation of the EM field. In a free region of space\footnote{Note $\bm{A}$ is given by the sum of the external $\bm{A}_\mathrm{ext}$ and the internal $\bm{A}_\mathrm{int}$ vector potential which is self-consistently generated by matter being composed of charged moving particles. It has been shown~\cite{Skolimowski_2020} that neglecting $\bm{A}_\mathrm{int}$ can be inaccurate in metals, namely for light frequencies of $\bm{A}_\mathrm{ext}$ typically in the \textit{mid} to \textit{near}-infrared region (roughly $30$-$300$ THz). However, given the light-matter coupling of interest being in the \textit{far}-infrared region -- i.e., roughly two orders of magnitude lower in frequency, only $\bm{A} = \bm{A}_\mathrm{ext}$ will be considered. }, i.e., in the absence of current densities, the EM vector potential $\bm{A}$ in the Coulomb gauge $\bm{\nabla} \cdot \bm{A}(\bm{x},t) = 0$ is described by the wave-equation \begin{equation} \label{eq:em-wave-eq} -\nabla^2 \bm{A}(\bm{x},t) + \frac{1}{c^2} \pd[2]{}{t}\bm{A}(\bm{x},t) = \bm{0}~, \end{equation} with the total energy of the system given by \begin{equation} \label{eq:em-field-energy} H_\mathrm{em} = \frac{1}{2}\int\dif \bm{x}\,\del{\varepsilon_0 \bm{E}(\bm{x}, t)^2 + \frac{1}{\mu_0} \bm{B}(\bm{x}, t)^2}~, \end{equation} where $\bm{E}(\bm{x}, t) = -\pd{}{t}\bm{A}(\bm{x}, t)$ and $\bm{B}(\bm{x}, t) = \bm{\nabla} \times \bm{A}(\bm{x}, t)$. Considering the field to be contained in a finite volume $V$ and imposing periodic boundary conditions, the EM field can be \textit{box quantised}~\cite{Garrison_2008}. The vector potential is expanded as a sum of the discrete modes and separated into two complex terms \begin{equation} \label{eq:em-wave-ansatz} \bm{A}(\bm{x},t) = \sum_{\bm{k}s}\bm{e}_{\bm{k}s}A_{\bm{k}s}(\bm{x},t) = \sum_{\bm{k}s}\bm{e}_{\bm{k}s}\sbr{A_{\bm{k}s}(t)\frac{e^{i \bm{k} \cdot \bm{x}}}{\sqrt{V}} + \hc}~, \end{equation} where $s$ denotes the polarisation and $\bm{e}_{\bm{k}s}$ are unit polarisation vectors which satisfy the orthogonality condition $\bm{e}_{\bm{k}s} \cdot \bm{e}_{\bm{k}s'}=\delta_{ss'}$ and the Coulomb gauge condition $\bm{k} \cdot \bm{e}_{\bm{k}s} = 0$ in reciprocal space. Re-scaling the variable $A_{\bm{k}s}(t) = \sqrt{\frac{\hbar \omega_{\bm{k}}}{2 \varepsilon_0}} a_{\bm{k}s}(t)$, the solution of~\eqref{eq:em-wave-eq} \begin{equation} \del{\pd[2]{}{t} + \omega_{\bm{k}}^2} a_{\bm{k}s}(t) = 0~, \end{equation} with $\omega_{\bm{k}} = c\envert{\bm{k}}$, results~\cite{Garrison_2008} in an expression for the EM vector potential \begin{equation} \label{eq:classical-em-field} \bm{A}(\bm{x}, t) = \sum_{\bm{k}s} \sqrt{\frac{\hbar}{2 \omega_{\bm{k}} \varepsilon_0 V}} \sbr{e^{i\del{\bm{k}\cdot\bm{x} - \omega_{\bm{k}} t}} \bm{e}_{\bm{k}s}a_{\bm{k}s} + \hc} \coloneqq \sum_{\bm{k} s} \sbr{\tilde{\bm{A}}_{\bm{k} s}(\bm{x},t) a_{\bm{k} s} + \hc}~, \end{equation} identical to a collection of quantised simple harmonic oscillators. % \footnote{One should note, however, that $a_{\bm{k}s}(t)$ should not be seen as the wave-function of a photon~\cite{CohenTannoudji1997} as its Fourier transform is not the photon wave-function in real-space -- which cannot be formulated for massless vector fields~\cite{Weinberg_1995_ed}. %-- and it does not have a Schr\"odinger-form in the presence of sources. % The Green function $\Delta(\bm{x},\bm{y})$ of a massive vector field has singularities for $m \to 0$ as seen in $\Delta_{\mu\nu}(\bm{x},\bm{y}) = \int \frac{\dif^4 q}{(2\pi)^4} e^{i \bm{q} \cdot (\bm{x} - \bm{y})} \frac{\eta_{\mu \nu} + q_\mu q_\nu/m^2}{q^2 + m^2 - i \epsilon}$.} The total energy~\eqref{eq:em-field-energy} then reads \begin{equation} H_\mathrm{em} = \sum_{\bm{k}s} \frac{\hbar \omega_{\bm{k}}}{2}\del{a^*_{\bm{k}s} a_{\bm{k}s} + a_{\bm{k}s} a^*_{\bm{k}s}}~, \end{equation} which is mathematically equivalent to the energy of a set of harmonic oscillators with frequency $\omega_{\bm{k}}$. Given the formal equivalence to oscillators, the quantisation of the EM field can heuristically be accomplished by promoting $a_{\bm{k}s}^*$ and $a_{\bm{k}}$ to mutually adjoint operators satisfying bosonic commutation relations \begin{equation} \sbr{\hat a_{\bm{k}s}, \hat a_{\bm{k'}s'}} = \sbr{\hat a^\dagger_{\bm{k}s}, \hat a^\dagger_{\bm{k'}s'}} = 0~,\qquad \sbr{\hat a_{\bm{k}s}, \hat a^\dagger_{\bm{k'}s'}} = \delta_{\bm{k}\bm{k'}}\delta_{ss'}~. \end{equation} These are the creation and annihilation operators of photons, the quanta (or excited states) of the quantised EM field, with total energy \begin{equation} \hat H_\mathrm{em} = \sum_{\bm{k}s} \hbar \omega_{\bm{k}} \del{\hat a^\dagger_{\bm{k}s} \hat a_{\bm{k}s} + \frac{1}{2}}~. \end{equation} \subsection{Anderson lattice model coupled to terahertz light} \label{sec:alm-light} \subsubsection{Light-matter coupling (dipole gauge)} The interaction between the EM and matter fields can be introduced via the \textit{minimal coupling}, where the transformation of the momentum operator \begin{equation} -i \hbar \bm{\nabla} \to -i \hbar \bm{\nabla} + e \bm{A} \end{equation} ensures Lorenz gauge invariance of the light-matter theory~\cite{Weinberg_1995_ed}. The projection of the continuum theory of electronic matter fields $\Psi(\bm{x})$ interacting with the EM field \begin{equation} \label{eq:hamiltonian-light-matter} \hat H = \hat H_\mathrm{em} + \int \dif\bm{x}\, \hat \Psi^\dagger(\bm{x}) \sbr{\frac{(-i \hbar \bm{\nabla} + e \bm{A})^2}{2m} + V(\bm{x})} \hat \Psi(\bm{x}) + \ldots~, \end{equation} where into a low-energy electronic basis $\hat \Psi^\dagger(\bm{x}) = \sum_{i \mu} \psi_{i \mu}(\bm{x}) \hat c_{i \mu}$, where $\psi_{i \mu}(\bm{x})$ denotes electronic Wannier functions and $\hat c_{i \mu}$ the annihilation operator of electrons with orbital $\mu$ localised at site $i$, hides subtle issues regarding the preservation of gauge invariance~\cite{Li_2020}. It is advised~\cite{Dmytruk_2021} first to project the matter fields and then add the minimal coupling through the unitary transformation $\hat H \to \hat U^\dagger \hat H \hat U$, where $\hat U = \exp\sbr{i \sum_{\bm{k}s} (\hat a_{\bm{k}s} + \hat a_{\bm{k}s}^\dagger) \int \dif\bm{x}\, \hat \Psi^\dagger(\bm{x}) \chi_{\bm{k}s}(\bm{x}) \hat \Psi(\bm{x})}$ and $\bm{\nabla} \chi_{\bm{k}s}(\bm{x}) = e \tilde{\bm{A}}_{\bm{k}s}(\bm{x})$. As a result, the EM field couples linearly to matter \begin{equation} \hat H = \hat H_\mathrm{em} + i \sum_{\bm{k}s} \sum_{ij} \sum_{\mu \nu} \hbar \omega_{\bm{k}} \chi_{\bm{k}s}^{\mu \nu, ij} \del{\hat a_{\bm{k}s} - \hat a^\dagger_{\bm{k}s}} \hat c^\dagger_{i \mu} \hat c_{i \nu} + \hat H_\mathrm{matter}~, \end{equation} where $\chi_{\bm{k}s}^{\mu \nu, ij} = \int \dif\bm{x}\, \psi^*_{i \mu}(\bm{x}) \chi_{\bm{k}s}(\bm{x}) \psi_{j \nu}(\bm{x})$, and can induce and intra- or inter-site photon-mediated hybridisation/hopping terms as well, to a smaller degree, shift the orbital energies. In addition, the EM field generates a quartic, self-interaction term between matter fields~\cite{Dmytruk_2021}, which was dropped due to being relevant only for strong light-matter coupling regimes. \subsubsection{Dipole selection rules} Assuming that the wavelength of the EM field is much larger than the atom size -- atomic radii are of the order of $10^{-10}$m, whereas THz radiation is of the order $10^{-4}$m -- to a good approximation, the EM field does not change over the atomic scale. As such, $\chi_{\bm{k} s}(\bm{x}) \approx e \tilde{\bm{A}}_{\bm{k} s}(\bm{x}) \cdot \bm{x}$ and the light-matter matrix-element reads \begin{equation} \label{eq:light-matter-coupling} \chi_{\bm{k}s}^{\mu \nu, ij} = e \sqrt{\frac{\hbar}{2 \omega_{\bm{k}} \varepsilon_0 V}} \bm{e}_{\bm{k}s} \cdot \int \dif\bm{x}\, \psi^*_{i \mu}(\bm{x})\, \bm{x}\, \psi_{j \nu}(\bm{x}) \end{equation} Due to the position operator in the integral kernel, for same-site coupling $i=j$, to a good approximation, light couples only matter in orbitals with different parity, in agreement with the selection rule $\Delta L = 1$ for light-induced transitions in the \enquote{dipole approximation}. Considering the $d$ conduction bands and $4f$ bands of \cref{sec:anderson-lattice-model}, the coupling of the EM field induces on-site transitions between the two orbitals. Neglecting effects of light-induced inter-site hopping, the coupling of the EM field to the Anderson lattice model~\eqref{eq:ALM-ap} reads \begin{equation} \label{eq:AndersonLightModel} \hat H = \sum_{\bm{k}s} \hbar \omega_{\bm{k}} \hat a^\dagger_{\bm{k}s} \hat a_{\bm{k}s} + \hbar g_0 \sum_{\bm{k} s} i \del{\hat a_{\bm{k} s} - \hat a^\dagger_{\bm{k} s}} \sum_{i \sigma} \del{\hat c^\dagger_{i \sigma} \hat b^\dagger_i \hat f_{i \sigma} + \hc} + \hat H_{\mathrm{ALM}}~, \end{equation} where an approximation similar to Weisskopf-Wigner's~\cite{Scully_1997}, where the light-matter coupling strength is assumed constant for the frequencies of interest $g^{\mu\nu,ij}_{\bm{k} s} \coloneqq \omega_{\bm{k} s} \chi^{\mu\nu,ij}_{\bm{k} s} \stackrel{!}{=} g_0 \delta_{d, f}\delta_{ij}$ is considered. This is ultimately related to the light-matter interaction, to a good approximation, not transfering momentum between electrons, as the momentum carried by a THz photon $q \sim 10^{-1}$ m$^{-1}$ is much smaller than the typical Fermi momentum $k_F \sim 10^6$ m$^{-1}$ -- for the interaction expressed in momentum basis $\hat H_\mathrm{int} \sim i \sum_{\bm{k} \bm{q}\sigma} \del{\hat a_{\bm{q}} - \hat a^\dagger_{-\bm{q}}} \del{\hat c^\dagger_{\bm{k}+\bm{q}\sigma} \hat f_{\bm{k}\sigma} + \hc} \sim i \sum_{\bm{q}} \del{\hat a_{\bm{q}} - \hat a^\dagger_{-\bm{q}}} \sum_{\bm{k}\sigma }\del{\hat c^\dagger_{\bm{k}\sigma} \hat f_{\bm{k}\sigma}}$. \section{Dissipative dynamics} \label{sec:heat-baths} A difficulty encountered in the theoretical description of driven non-equilibrium systems is that the work done by external fields is transformed into heat through quantum scattering mechanisms~\cite{Amaricci_2012}. Beyond the low-energy toy models used to describe strongly interacting systems, real materials interact further with \textit{dissipative}\footnote{For a proper understanding of the meaning of \enquote{dissipative} channels, one has to introduce \textit{open} quantum systems and their dynamics, which is beyond the intended scope of this thesis. In modelling an open system, the system of interest $S$ is coupled to a \textit{reservoir} $R$, an environment with infinite degrees of freedom. Despite the system plus reservoir defining a \textit{closed} system -- with unitary dynamics generated by their joint Hamiltonian, it can be assumed that the system $S$ will have a negligible influence on the reservoir $R$, which can be traced out~\cite{Breuer_2007, Gardiner_2004, Wiseman_2009}. The elimination of these degrees of freedom results in an \textit{open} description of the system $S$, with irreversible dynamics, due to information loss -- or \textit{decoherence} -- to the environment.} channels such as background photons, phonons or other electronic bands, which can remove the excess energy deposited by external fields. Adding an environment can be critical in controlling the population of high-energy states or maintaining the system's internal energy conserved on average. This is especially important in \textit{continuously} driven systems, where a balance between the action of the driving field and the dissipation of high-energy electrons is necessary to reach some kind of stationary state. For a \textit{short} drive such as an ultrafast pulse, a finite amount of energy is deposited in the system, and stationarity can be reached without coupling to an environment. However, depending on the intensity of the injected pulse, the resulting heating of the system may be unwanted if thermal fluctuations entirely wash out low-energy quantum effects. A reservoir in thermal equilibrium -- or heat bath -- is typically modelled~\cite{Amaricci_2012, Werner_2012} as a set of independent fermionic or bosonic modes described by a time-dependent Gibbs ensemble~\eqref{eq:gibbs-ensemble} and coupled through exchange terms to the system's electronic degrees of freedom. In a scenario where the environment is a heat bath, and the system's Hamiltonian becomes static, its stationary state is likely to be described by a Gibbs ensemble at the same temperature as the reservoir~\cite{Breuer_2007}, a process known as \textit{thermalisation}. The timescale at which a system thermalises depends on the possible channels of thermalisation, which not only depend on the system-bath coupling strength but also on momentum, energy or particle conservation constraints specific to the material. Beyond the canonical way of treating open quantum systems through master equations\footnote{\label{fn:markov-master}Under a \textit{Markovian} approximation of the reservoir $R$, with short correlation times, and the assumption of it not being significantly affected by the interaction with the system $S$ (for which the factorisation of the total density matrix $\hat \rho(t) \stackrel{!}{=} \hat \rho_S(t) \otimes \hat \rho_R$ can hold), the reservoir degrees of freedom are traced out~\cite{Breuer_2007, Gardiner_2004, Wiseman_2009} resulting in the quantum master equation \begin{equation*} \partial_t \hat \rho_S = \hat{\mathcal{L}}\, \hat \rho_S = -i \sbr{\hat H, \hat \rho_S} + \mathcal{D}[\mathcal{L}] = -i \sbr{\hat H, \hat \rho_S} + \sum_\alpha \kappa_\alpha \del{\hat L_\alpha \hat \rho_S \hat L_\alpha^\dagger - \frac{1}{2}\cbr{\hat L_\alpha^\dagger \hat L_\alpha, \hat \rho_S}}~. \end{equation*} While the first part of the \textit{Liouvillian} $\hat{\mathcal{L}}$ describes the unitary/coherent dynamics generated by Hamiltonian $\hat H$ -- as encountered in the von Neumann equation, the second part $\mathcal{D}[\hat{\mathcal{L}}]$ describes the \textit{dissipative} dynamics that result from the connection to the environment. The jump operator $\hat L_\alpha$ is the remnant system operator which was coupled to the environment and $\kappa_\alpha$ encodes an effective coupling strength to the environment. Note that a series of other approximations were employed in obtaining the $\kappa_\alpha$, which generally can be time-dependent. The transformation from the quantum master equation to a path integral formulation follows a similar procedure to \cref{sec:path-integral} and can be found in detail in~\cite{Sieberer_2016}. In short, the dissipator $\mathcal{D}[\hat{\mathcal{L}}]$ enters the Keldysh action~\eqref{eq:keldysh-action} as \begin{equation*} \sum_\alpha \kappa_\alpha \cbr{L_{\alpha}(t_-)L^*_{\alpha}(t_+) - \frac{1}{2}\sbr{L^*_{\alpha}(t_-)L_{\alpha}(t_-) + L^*_{\alpha}(t_+)L_{\alpha}(t_+)}}~, \end{equation*} and couples the forward and backward branches of the contour, ultimately breaking the cancellation of the forward and backward action of the contour under unitary time evolution.}, \textit{irreversible} dynamics can be obtained by integrating out the reservoir and replacing all reservoir-related 2-point functions in~\eqref{eq:gamma2-functional} by \textit{non-interacting} 2-point functions. Reservoirs are, by construction, not renormalised by interactions with the system, and such a procedure will enforce the self-energy of reservoir 2-point functions to vanish, effectively blocking those degrees of freedom from being affected by the interaction with the system. In opposition, the system is coupled to and renormalised by the reservoir. Hence, due to the asymmetrical nature of the system-reservoir coupling, information from the system will be lost to the reservoir, generating irreversible dynamics. \subsection{Fermionic heat bath} \label{sec:fermionic-bath} The B\"uttiker model~\cite{Murakami_2017} couples an identical fermionic reservoir to each lattice site \begin{equation} \label{eq:fermionic-bath} \hat H_{\mathrm{bath}} = \sum_\ell \sum_{i \sigma} \varepsilon^a_{\ell} \hat a^\dagger_{\ell i \sigma} \hat a_{\ell i \sigma} + \alpha \sum_\ell \sum_{i \sigma} \del{\hat c^\dagger_{i \sigma} \hat a_{\ell i \sigma} + \hc}~, \end{equation} where $\hat a^\dagger_{\ell i \sigma}$ is a reservoir electron creation operator of a fermion with spin $\sigma$ at site $i$ with energy $\varepsilon_{\ell}$. The heat bath is, by definition, site and spin diagonal, and upon integration, generates the additional \textit{non-interacting} effective action term \begin{equation} \label{eq:fermionic-bath-action} i S_{\mathrm{eff}}^{\mathrm{bath}} = \int_\gamma \dif z \dif z'\, \sum_{i \sigma} c^*_{i \sigma}(z) \sbr{\alpha^2 \Delta_{a}(z, z')} c_{i \sigma}(z')~, \end{equation} which is effectively an additional hybridisation of conduction electrons. The lack of momentum conservation resulting from such coupling is crucial to allow the electrons' momenta to relax from the excitation by external fields. The heat bath 2-point function reads \begin{equation} \Delta_{a}(z, z') = \int_{-\infty}^{+\infty} \dif \varepsilon\, \rho^a_{0}(\varepsilon) \sbr{-i \del{\Theta_\gamma(z, z') - f(\varepsilon)}e^{-i \varepsilon (z - z')}}~, %\delta_{ij} \delta_{\sigma\sigma'}~, \end{equation} where $\rho^a_{0}(\varepsilon) = \sum_\ell \delta(\varepsilon - \varepsilon^a_\ell)$ is the bath density of states and $f$ is the Fermi-Dirac distribution. It is important that the environment is infinitely large and involves a continuum of frequencies for the bath Green functions to decay in time (cf.~\cref{sec:dmft-aux-particles}) and properly generate irreversible dynamics~\cite{Breuer_2007}. For a small system-bath coupling strength $\alpha$, the effect of the bath on the original conduction electron bandwidth is negligible and the details of $\rho^a_{0}(\varepsilon)$ unimportant. As such, it is typically taken as a flat density of states or being the same as the conduction electrons' bare one. Due to the particle-exchange character of the coupling, the reservoir can also tune the total electronic particle number of the system. \subsection{Bosonic heat bath} Another way to incorporate dissipation into the electronic system is via baths of bosonic particles -- typically phonons or photons. This kind of coupling allows energy exchange with the environment but, unlike the former, microscopically conserves electronic particle number due to the conservation of electronic charge. The resulting system-bath coupling is hence an interacting vertex and can introduce more interesting thermalisation transients~\cite{Peronaci_2020} than their fermionic counterpart. As such, the inclusion of bosonic bath channels requires perturbative expansions, and hence it will be assumed that the system-bath coupling strength $\alpha$ is small\footnote{Refer to~\cite{Chen_2016} for strong system-bath couplings.}. \paragraph{Ohmic bath} The Holstein model adds a minimal coupling between phonon modes and the electronic degrees of freedom and is used to describe the effects of vibrations in such systems~\cite{Chen_2016} \begin{equation} \hat H_{\mathrm{bath}} = \sum_\ell \sum_{i \sigma} \hbar \omega^a_{\ell} \hat a^\dagger_{\ell i} \hat a_{\ell i} + \alpha \sum_\ell \sum_{i \sigma} \del{\hat c^\dagger_{i \sigma} \hat c_{i \sigma} - \hat{\mathds{1}}}\del{\hat a^\dagger_{\ell i} + \hat a_{\ell i}}~, \end{equation} where $\hat a^\dagger_{\ell i}$ is a reservoir phonon creation operator of a particle with frequency $\omega_{\ell}$ at site $i$ and $\alpha$ is the coupling strength between the phonon modes and the conduction electrons. By integrating out the phonon degrees of freedom, a dissipative channel enabling energy and momentum relaxation arises through a new \textit{interacting} effective action term \begin{equation} i S_{\mathrm{eff}}^{\mathrm{bath}} = \int_\gamma \dif z \dif z' \sum_{i \sigma} \sbr{c_{i \sigma}^*(z) c_{i \sigma}(z) - 1} i \alpha^2 \sbr{\Pi_{a}(z, z') + \Pi_{a}(z', z)} \sbr{c_{i \sigma}^*(z') c_{i \sigma}(z') - 1}~, \end{equation} for which the lowest order loop expansion results in a self-energy contribution to the conduction electrons~\cite{Eckstein_2013, Grandi_2021}. The heat bath 2-point function reads \begin{equation} \label{eq:bosonic-thermal-gf} \Pi_{a}(z, z') = \int_{-\infty}^{+\infty} \dif \varepsilon\, \rho^a_0(\varepsilon) \sbr{-i \del{\Theta_\gamma(z, z') + b(\varepsilon)}e^{-i \varepsilon (z - z')}}~, %\delta_{ij}~, \end{equation} where \begin{equation} \label{eq:ohmic-dos} \rho^a_{0}(\varepsilon) = \sum_\ell \delta(\varepsilon - \omega^a_\ell) = \varepsilon \frac{\Lambda^2}{\varepsilon^2 + \Lambda^2} \Theta(\varepsilon) \end{equation} is an ohmic density of states, linear at low frequencies and with a soft cut-off $\Lambda$ and $b$ is the Bose-Einstein distribution function. While in normal conductors, it is known that the thermalisation time is generally within dozens/hundreds of femtoseconds, it can go up to 100 picoseconds in heavy-fermion compounds. It has been hypothesized~\cite{Demsar_2006} that this is due to the suppression of electron-phonon scattering in the vicinity of the Fermi energy. Due to the flat bands near the Fermi energy, the Fermi velocity can be slower than the sound velocity, which in turn, suppresses the phase space available for the elastic scattering between electrons and phonons, dramatically increasing the thermalisation time via phonon scattering. \subsection{Photonic heat bath} \label{sec:photonic-bath} Although the hypothesised suppression of phononic scattering in heavy-fermion compounds, another possible channel of bosonic thermalisation is through spontaneous emission. Namely, in heavy-fermion systems, electrons in the conduction band can fill holes in the valence band via the spontaneous emission of a photon, a process known as \textit{radiative recombination}. This exchange arises due to the EM field coupling electronic states of different orbitals and, due to the small momentum carried by photons, transitions effectively not transferring momenta (\cref{sec:alm-light}). \subsubsection{Continuous-mode light} In quantum optics, the EM field is usually confined within a cavity, for which \textit{box quantisation} results in a complete set of discrete modes. However, for optical experiments such as THz spectroscopy, there is no identifiable cavity, but rather light flowing from a source through some interaction region to a set of detectors. In the absence of a cavity, physical quantities such as the light-matter coupling strength~\eqref{eq:light-matter-coupling} should hence not be determined by the quantisation volume. Moreover, real experiments are more accurately described by Gaussian beams of light instead of plane waves, and more quantitative calculations require a vastly different quantisation procedure~\cite{Deutsch_1993}. To quantise the EM field in \textit{free space}~\cite{Blow_1990}, the quantisation volume goes to infinity $V\to\infty$, such that the sum over the discrete modes \begin{equation} \sum_{\bm{k}} \to \frac{1}{\Delta k}\int\dif\bm{k}~\, \end{equation} is converted to an integral, where $\Delta k = \frac{(2\pi)^3}{V}$, and the discrete-mode creation and annihilation operators are related to \textit{continuous-mode} ones through \begin{equation} \hat a_{\bm{k}} \to \sqrt{\Delta k} \hat a(\bm{k}) \qquad \hat a^\dagger_{\bm{k}} \to \sqrt{\Delta k} \hat a^\dagger(\bm{k})~, \end{equation} satisfying the continuous-mode commutation relation \begin{equation} \sbr{\hat a(\bm{k}), \hat a^\dagger(\bm{k}')} = \delta(\bm{k} - \bm{k}')~. \end{equation} The continuous-mode vector potential~\eqref{eq:classical-em-field} reads \begin{equation} \hat{\bm{A}}(\bm{x}, t) = \int \frac{\dif\bm{k}}{(2\pi)^{3/2}}\sum_{s} \sqrt{\frac{\hbar}{2 \omega(\bm{k}) \varepsilon_0}} \sbr{e^{i\del{\bm{k}\cdot\bm{x} - \omega(\bm{k}) t}} \bm{e}_{s}(\bm{k})\hat a_{s}(\bm{k}) + \hc}~, \end{equation} and partitioning the integral into sections of the solid angle $\int \dif\bm{k} \to \int \dif k\,k^2\sum_m \int \dif\Omega_m$, with area $\Delta \Omega_m \equiv \int_{\Omega_m} \dif \Omega$, results a decomposition of one-dimensional modes at $\bm{x}=0$ in the limit $\Delta_m\to0$~\cite{Ko_2022} \begin{equation} \label{eq:em-field-steradian} \hat{\bm{A}}(t) = \sum_{m s} \int_0^{\infty} \dif \omega \sqrt{\frac{\hbar \omega \Delta\Omega_m}{16 \pi^3 \varepsilon_0 c^3}} \sbr{e^{-i \omega t} \bm{e}_{ms} \hat a_{ms}(\omega) + \hc}~, \end{equation} where $\hat a_{ms}(\omega)$ is an annihilation operator of a mode with frequency $\omega$, solid angle section $m$ and polarisation $s$. The bare Hamiltonian of the continuous-mode EM field reads \begin{equation} \label{eq:free-em-hamiltonian} \hat H_\mathrm{em} = \sum_{ms} \int_{0}^{\infty} \dif\omega\, \hbar\omega \hat a^\dagger_{ms}(\omega) \hat a_{ms}(\omega)~, \end{equation} where the vacuum energy was removed. The EM field vacuum (a state devoid of photons) contains an infinite amount of energy $E_\mathrm{vac} = \frac{1}{2}\sum_{ms} \int_{0}^{\infty} \dif\omega \hbar \omega$ and is a result of vacuum fluctuations. However, this superficial divergence in free space is trivially removed as it simply defines the reference from which all energies must be calculated. \subsubsection{Photonic heat bath} The quantum nature of the EM field in free space generates an infinite number of modes that serve as a reservoir to which atoms can radiate. This kind of irreversible disappearance of light quanta cannot occur in a box since light cannot be lost by bouncing back (either physically or theoretically through periodic boundary conditions). In free space, \eqref{eq:free-em-hamiltonian} forms the bath into which the system~\eqref{eq:ALM-ap} can radiate to \begin{equation} \label{eq:photonic-bath} \hat H_{\mathrm{bath}} = \hat H_\mathrm{em} + \hbar g_0 \sqrt{\eta} \sum_{m s} \int \dif\omega\, i \del{\hat a_{m s}(\omega) - \hat a^\dagger_{m s}(\omega)} \sum_{i \sigma} \del{\hat c^\dagger_{i \sigma} \hat b^\dagger_i \hat f_{i \sigma} + \hc}~, \end{equation} where $\eta$ is a dimensionless parameter characterising the strength of the coupling to the EM field environment. Note that the prefactors of $g_0$ in~\eqref{eq:photonic-bath} are slightly different from the ones in~\eqref{eq:AndersonLightModel}. Integrating out the photonic degrees of freedom generates the new \textit{interacting} effective action term \begin{equation} \label{eq:photonic-bath-action} 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} \hbar^2 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}$}~, \end{equation} where $\Pi_a(z, z')$ is described by~\eqref{eq:bosonic-thermal-gf}. Despite $\Lambda\to\infty$ in the EM field vacuum, a soft cut-off to the density of states is kept because the light-matter coupling was assumed to be constant in frequency, which must be remedied by a cut-off in the coupling to high-energy vacuum modes. \section{Travelling pulse of light} \label{sec:travelling-pulse} When a travelling pulse of quantum radiation interacts with non-linear quantum systems in free space, its photon number content or temporal mode (\enquote{wavepacket shape}) may change due to quantum absorptive or dispersive effects. The main difficulty of treating such systems is that due to spontaneous emission, not only can the structure of the temporal mode change, the emission is, in general, not restricted to that specific mode~\cite{Silberfarb_2003}. The theoretical frameworks that treat the interaction of such travelling pulses of quantum radiation with quantum systems date back to the influential input-output theory~\cite{Gardiner_1985}. Here, the \textit{output} radiation (asymptotically free after interaction) can be related to the \textit{input} radiation (asymptotically free before interaction) through the simple formula \begin{equation} \hat a_\mathrm{out}(t) = \hat a_\mathrm{in}(t) + g\,\hat c(t)~, \end{equation} where $g$ denotes an effective coupling between some system ($c$) and radiation ($a$). However\footnote{Recent extensions~\cite{Kiilerich_2019} to this theory can describe an incident mode of radiation acting on some non-linear interacting system. The key idea is that a cavity driven by white noise (i.e., coupled to a Markovian environment) can be engineered to generate coloured noise~\cite{Gough2012}, i.e. any arbitrary input radiation can be generated by choosing appropriate jump operators and manipulating the pump/loss rates (\cref{fn:markov-master}). A pulse of radiation interacting with a cavity is then treated as an open (cascaded~\cite{Combes2017}) quantum system, where an upstream pseudo-cavity leaks the \textit{input} radiation that drives the cavity-system of interest, with the \textit{output} radiation mode being acquired by some downstream pseudo-cavity. While this technique should, in principle, be able to model the interaction of a travelling pulse with an Anderson lattice model, the two pseudo-cavities mentioned above can only describe interaction with a single mode of the incident and outgoing radiation, with the rest of the modes being treated as loss. Hence, it can be used to examine the quantum contents of a specific output temporal mode after interaction with the system of interest. While the theory can be extended~\cite{Kiilerich_2020} to multi-input and output modes, it becomes prohibitively expensive as the number of modes grows. % Hence, it is best suited for systems in actual cavities, with a well-defined mode instead of free space. % It would, nonethless, be interesting to investigate the model~\eqref{eq:ALM-ap} driven by an input pulse with shape $xi$ with the theory presented in~\cite{Kiilerich_2020}. Here, where the quantum master equation (\cref{fn:markov-master}) of the system plus radiation pulse is described by the Hamiltonian and % \begin{equation*} % \hat H = ~, % \end{equation*} % and time-dependent jump operators % \begin{equation*} % \hat L(t) = g_\xi^*(t) \hat a_\xi + g_0 % \end{equation*} }, this requires a series of considerations~\cite{Gardiner_1985} incompatible with the system~\eqref{eq:ALM-ap}, namely the rotating wave approximation -- in continuum systems, there are no characteristic slow and fast frequencies. Since the non-equilibrium field theory was framed concerning 2-point functions, related with ensemble averages, and not specific quantum states, the modelling of a travelling pulse in free space needs only to be resolved at the level of Green functions and can hence be reverse-engineered through the investigations of~\cref{sec:neqft-interpretability}. % For pulse modes distributed over the Gaussian density of states % \begin{equation} % \rho_\mathrm{pulse}(\omega) = \sqrt{\frac{2\pi}{\Omega_0^2}} \exp{\sbr{- \frac{\del{\omega - \omega_0}^2}{2 \Omega_0^2}}}~, % \end{equation} % where $\Omega_0$ denotes the frequency width and $\omega_0$ the central, carrier frequency, the forward EM field reads % \begin{equation} % \hat E^+(\bm{x}, t) \approx \int\frac{\dif \omega}{2\pi} \rho_\mathrm{pulse}(\omega) \sum_{s} \sqrt{\frac{\hbar\omega}{2 \varepsilon_0 V}} \bm{e}_{s}(\omega) a_{s}(\omega) e^{i\del{\bm{k}\cdot\bm{x} - \omega t}} \approx \sqrt{\frac{\hbar \omega_0}{2 \varepsilon_0 V}} e^{-i \omega_0 t} \xi(T - T_0)~, % \end{equation} Assuming a one-dimensional Gaussian pulse to be travelling along some direction, with a carrier frequency $\omega_0$ and Gaussian envelope~\cite{Wang2011} $\xi_{\Omega_0}(t)$ centred at and $T_0$ \begin{equation} \xi_{\Omega_0}(t) = \exp\del{-\frac{\Omega_0^2}{4}t^2}~, \end{equation} where $\Omega_0$ denotes the frequency width of the pulse, the electric field pulse at $\bm{x}=0$ (cf.~\eqref{eq:em-field-steradian}) reads \begin{equation} \label{eq:quantum-em-field} \hat{\bm{E}}(t) \approx \sum_{s} \sqrt{\frac{\hbar\omega_0^3 \Delta\Omega_\mathrm{P}}{16 \pi^3 \varepsilon_0 c^3}} e^{-i \omega_0 t} \xi_{\Omega_0}(t - T_0) \bm{e}_s(\omega_0) \hat a_s(\omega_0) + \hc~, \end{equation} assuming enveloping bandwidths much smaller than the carrier frequency and $\Delta\Omega_\mathrm{P}$ the solid angle covered by the pulse. % The annihilation operator of the photon pulse mode % \begin{equation} % \hat a_\xi = \int \frac{\dif \omega}{2\pi} \sbr{\int \dif t\, e^{+i \omega t} \del{\xi_{\Omega_0}(t - T_0) e^{-i \omega_0 t}}}\hat a(\omega)~, % \end{equation} % satisfies canonical commutation relations $\sbr{\hat a_\xi, \hat a^\dagger_\xi} \stackrel{!}{=} 1$, assuming enveloping bandwidths much smaller than the carrier frequency. The intensity at $\bm{x}=0$, determined by the Poynting operator, for linearly polarised light~\cite{Loudon2000} \begin{equation} \hat I(t) = 2 \varepsilon_0 c \hat{\bm{E}}^{-}(t) \cdot \hat{\bm{E}}^{+}(t)~, \end{equation} where $\hat{\bm{E}}^{\pm}(t)$ denotes the forward and backward components of EM field of~\eqref{eq:quantum-em-field}, respectively, is proportional to the number of photons \begin{equation} \hat I(t) = \frac{\hbar \omega_0}{\tilde{A}} \hat a^\dagger(\omega_0) \hat a(\omega_0) \del{\xi_{\Omega_0}(T - T_0) e^{-i \omega_0 t}}^* \del{\xi_{\Omega_0}(T - T_0) e^{-i \omega_0 t}}~, \end{equation} where $\tilde{A} = \frac{4 \pi^3 c^2}{\Delta\Omega_P \omega_0^2}$. A generalisation of the Poynting operator to two times results in an expression for the lesser Green function of the pulse $i \left. G^<_{a_0}(t, t')\right|_{t'\to t} \sim \frac{\tilde{A}}{\hbar \omega_0}\ev{\hat I(t)}$. In Wigner coordinates, it reads \begin{equation} \label{eq:ga-lesser} G^<_{a_0}(T, \tau)_{\mathrm{W}} = -i \bar{n}_a \xi_{\sqrt{2}\Omega_0}(T - T_0) e^{-i \omega_0 \tau} \xi_{\Omega_0 / \sqrt{2}}(\tau)~, \end{equation} where $\bar{n}_a$ is the maximum number of photons in the Gaussian mode, $\xi_{\sqrt{2}\Omega_0}(T - T_0)$ controls the modulation of the occupation and $e^{-i \omega_0 \tau} \xi_{\Omega_0 / \sqrt{2}}(\tau)$ controls the spectral content of the Green function -- its Fourier transform is a normalised Gaussian distribution with width $\Omega_0 / \sqrt{2}$ centred at frequency $\omega_0$\footnote{The spectral content of a photonic pulse being a Gaussian distribution contains the problem that for $\omega_0 < \Omega_0$, there is a significant portion of the photonic field with \textit{unphysical} zero and negative energies. To eliminate this spectral region, the spectral information can instead be encoded by the half-Fourier transform \begin{equation} e^{-i \omega_0 \tau} \xi_{\Omega_0 / \sqrt{2}}(\tau) \to \int_{0}^{+\infty} \frac{\dif \omega}{2\pi} e^{-i \omega \tau} \, \omega \int_{-\infty}^{+\infty} \dif \tau'\, e^{+i \omega \tau'} \sbr{e^{-i \omega_0 \tau'} \xi_{\Omega_0 / \sqrt{2}}(\tau')}~, \end{equation} which ensures the photons always have positive energies.}. Note that this implies that the pulse's distribution function is flat in frequency, which differs significantly from the ones found in thermal equilibrium, such as the Bose-Einstein distribution. While it would be possible to control the occupation by a time-dependent chemical potential and a Bose-Einstein distribution, this was not considered as that would mean that the field occupation has different spectral decompositions at different times. Since the pulse field can be considered non-interacting, the greater Green function \begin{equation} G^>_{a_0}(T, \tau)_{\mathrm{W}} = -i \del{1 + \bar{n}_a \xi_{\sqrt{2}\Omega_0}(T - T_0)} e^{-i \omega_0 \tau} \xi_{\Omega_0 / \sqrt{2}}(\tau)~, \end{equation} ensures the spectral content of the pulse is constant in centre-of-mass time. The resulting equations of motion read\footnote{These equations of motion have a striking similarity to those found in open quantum systems described by a quantum master equation (\cref{fn:markov-master}). Due to the difficulty of resolving the joint photon-system interactions, the standard technique is to integrate the photonic input field and study the resulting effective cavity. For a single-particle pumping $\hat L = \hat a^\dagger$ and loss $\hat L = \hat a$, with respective effective pump $\kappa_p$ and loss rate $\kappa_l$, the equations of motion for a non-interacting Hamiltonian $\hat H = \omega_0 \hat a^\dagger \hat a$ read \begin{equation*} \label{eq:gf-eom-open} \begin{split} i \partial_t G^<(t,t') &= \sbr{\omega_0 - \frac{i}{2}\del{\kappa_l(t) + \kappa_p(t)}}G^<(t,t') + i \kappa_p(t) \sbr{\Theta(t-t') G^<(t,t') + \Theta(t'-t) G^>(t,t')} \\ i \partial_t G^>(t,t') &= \sbr{\omega_0 + \frac{i}{2}\del{\kappa_l(t) + \kappa_p(t)}}G^>(t,t') - i \kappa_l(t) \sbr{\Theta(t-t') G^>(t,t') + \Theta(t'-t) G^<(t,t')}~. \end{split} \end{equation*} While it is not possible to compare them directly with~\eqref{eq:photon-eom}, it can be argued that the imaginary terms of the right-hand-side are similar to time-dependent loss/pump rates -- however, unlike the rates found in Markovian dynamics~\cite{Bacon_2001}, these are not necessarily positive-definite.} \begin{equation} \label{eq:photon-eom} \begin{split} i \partial_t G^<_{a_0}(t, t') &= % i \del{\frac{1}{2} \partial_T + \partial_\tau} G^<_{a_0}(T, \tau)_\mathrm{W} = \del{\omega_0 -i \frac{\Omega_0^2}{4} (t-t')} G^<_{a_0}(t, t') - i \frac{\Omega_0^2}{2} \del{\frac{t+t'}{2} - T_0} G^<_{a_0}(t, t') \\ i \partial_t G^>_{a_0}(t, t') &= % i \del{\frac{1}{2} \partial_T + \partial_\tau} G^>_{a_0}(T, \tau)_\mathrm{W} = \del{\omega_0 -i \frac{\Omega_0^2}{4} (t-t')} G^>_{a_0}(t, t') - i \frac{\Omega_0^2}{2} \del{\frac{t+t'}{2} - T_0} G^<_{a_0}(t, t') \end{split} \end{equation} Note that the travelling pulse of light is only resolved at the level of 2-point functions, and expressions for higher-point functions would require further analysis. The main challenge arising from a THz pulse interacting with a heavy-fermion system is the no separation of the timescales. The bandwidth of the pulse sets the duration of the interaction $\Omega_0 \not\ll k_\mathrm{B} T^*_K$, which is in the same scale as the Kondo coherence temperature. This results in proper non-equilibrium dynamics, which prevents the problem from being simplified, for example, by being cast into a quasi-static form (\cref{sec:wigner-basis}).