\chapter{Non-equilibrium saddle-point theory} \label{sec:thz-mf} Saddle-point approximations are ubiquitous in condensed-matter theory. At heart, these theories approximate the quantum averaging of the partition function~\eqref{eq:path-integral-Z} by a single \textit{classical} field contribution. That is, in the classical limit $\hbar \to 0$, the integral of the partition function \begin{equation} Z = \int D\sbr{\phi} e^{i S\sbr{\phi}/\hbar}~, \end{equation} is dominated by the extremum of $S\sbr{\phi}$, with the amplitudes from the other quantum paths destructively interfering with each other~\cite{Kleinert2009}. This results in Euler-Lagrange equations of motion for the now \enquote{classical} $\phi_0$ field, which satisfies the saddle-point condition \begin{equation} \left.\fd{S[\phi]}{\phi}\right|_{\phi=\phi_0} = 0~. \end{equation} Such theories have successfully described many complex phenomena, such as superconductivity~\cite{Sun_2020} or Bose-Einstein condensates~\cite{Lappe_2018}. The caveat is that the approximation is only valid if the amplitude of the quantum fluctuations $\delta \phi \coloneqq \phi - \phi_0$ around the classical field $\phi_0$ is small. In general, the saddle-point condition identically satisfies an associated mean-field theory, where a mean-field $\bar{\phi}$ extremises the free energy $F[\phi] = -\frac{1}{\beta} \log Z[\phi]$. Equivalently, the mean-field theory is only valid if the free energy associated with the fluctuations is much smaller than the mean-field free energy. The Ginzburg criterion dictates the validity of the mean-field approximation in equilibrium systems: fluctuations only become relevant when their length scale is longer than the system's correlation length. In non-equilibrium systems, however, much less is understood. Saddle-point theories were among the earliest successes in describing Kondo regimes via auxiliary-particle techniques (\cref{sec:aux-particles}). Seemingly counterintuitive, given the strongly-interacting nature of Kondo physics, there is a direct correspondence between the Kondo temperature and the fields' saddle points at zero temperature. This fortuitous coincidence makes saddle-point theory a powerful starting point when studying this class of problems as it offers a window into the underlying physics at a fraction of the cost of resolving the properly interacting\footnote{Note that even at the saddle-point, auxiliary-particle theories are, in general, \textit{strongly interacting} due to the constraint. What is meant is that the quantum fluctuations are neglected.} theory. However, it should be noted that auxiliary-particle saddle-point approximations suffer from severe fundamental problems, such as featuring a spurious phase transition (vanishing saddle-point solution) at finite temperatures~\cite{Kroha2004} or displaying a Fermi liquid phase in multi-channel (non-Fermi liquid) models. Nevertheless, for an Anderson lattice model in the $U\to\infty$ regime in thermal equilibrium, the saddle-point approximation and associated fluctuations are well-understood~\cite{Millis1987, Fr_sard_2011}. For Anderson-impurity models in non-equilibrium regimes, some saddle-point formulations have been proposed. However, these have been applied only to quasi-equilibrium problems, such as quantum-dot steady-states~\cite{Wu2008, Ratiani_2009}, high-frequency Floquet drive~\cite{Takasan_2017} and almost-adiabatic~\cite{Romeo_2010} drive. A proper solution in time -- which requires the problem to be formulated and solved in the Keldysh contour -- is still lacking. This is most likely due to the complicated form of the resulting equations of motion -- a set of differential-algebraic equations. Driven by the non-separable timescales that dictate the time-evolution of the system (\cref{sec:travelling-pulse}), a general non-equilibrium saddle-point theory of auxiliary particles for a time-resolved pulse of light interacting with an Anderson lattice is developed. \section{Non-equilibrium effective action} The partition function associated with the pulse of light interacting with an Anderson lattice has several field degrees of freedom that must be summed over. Ideally, the more these degrees of freedom are integrated, the more quantum effects are contained in the resulting effective action. Concerning the Hamiltonian of interest~\eqref{eq:full-model}, one could integrate the auxiliary bosons and derive an effective theory with boson-mediated electron interactions. However, this theory would require further transformations, such as the Hubbard-Stratonovich, to decouple the new effective interacting terms between matter fields. This is not only more laborious, but it also dismisses the great appeal of treating the auxiliary bosons and constraint at the saddle-point: their amplitude is directly proportional~\cite{Fr_sard_2011} to the Kondo temperature (\cref{sec:kondo-pheno}) of the system \begin{equation} \label{eq:kondo-temp-sp} \begin{split} V_0 \rho_0^c(\varepsilon_F) \envert{b_0}^2 &\approx T_K \qquad \textrm{for the single-impurity Anderson model}, \\ \lambda_0 - \varepsilon^f_0 &\sim T^*_K \qquad \textrm{for the Anderson lattice model}. \end{split} \end{equation} Fulfilling the constraint~\eqref{eq:q-constraint} at every lattice site $j$ can conveniently be imposed via the path-integral formalism, with the introduction of several Lagrange multipliers $\lambda_j$. In essence, a \textit{projection}~\cite{Fr_sard_2011} of the partition function onto the $U \to \infty$ subspace\footnote{Note that the action of the operator $\mathds{P}$ amounts to the Abrikosov projection (\cref{sec:ap-zeta-scaling}). By a change of variables to the fugacity, where $\zeta = e^{-i \lambda}$, the projector reads \begin{equation*} \mathds{P} \coloneqq \oint \frac{\dif \zeta}{2\pi i} \frac{1}{\zeta^2} f(\zeta) = \lim_{\zeta\to0} \od{}{\zeta}f(\zeta) = \lim_{\zeta\to0} \hat{Q}\zeta^{\hat{Q} - \hat{\mathds{1}}}~, \end{equation*} with $f(\zeta) = \zeta^{\hat Q}$. The projection of the partition function $\mathds{P} Z \equiv \lim_{\zeta \to 0} \hat{Q} \zeta^{\hat Q - \hat{\mathds{1}}} Z_G$ is hence entirely equivalent to~\eqref{eq:projection-Q}, modulo a scaling factor.} \begin{equation} \label{eq:projected-Z} Z_{\mathrm{physical}} = \del{\prod_j \mathds{P}_j} Z~, \end{equation} via some projector $\mathds{P}_j$ \begin{equation} \mathds{P}_j = \int_{-\pi}^{\pi}\frac{\dif\lambda_j}{2\pi} e^{-i \lambda_j \del{\hat{Q}_j - \hat{\mathds{1}}}}~, \end{equation} is required. Moreover, the introduction of the auxiliary particles adds a gauge degree of freedom \begin{equation} f_{j \sigma}(z) \to e^{i \phi_j(z)} f_{j \sigma}(z) \qquad b_{j}(z) \to e^{i \phi_j(z)} b_{j}(z)~, \end{equation} which is known~\cite{Arrigoni_1994} to produce divergences in the auxiliary-boson Green function. The radial gauge \begin{equation} \lambda_j \to \lambda_j - i \partial_z \phi_j(z) \eqqcolon \lambda_j(z)~, \end{equation} promotes $\lambda$ to a field while leaving the action invariant under gauge transformations \begin{equation} S[\cdot, \lambda_j(z)] \to S[\cdot, \lambda_j(z) - i \partial_z \phi_j] = S[\cdot, \lambda(z)] + \int_\gamma \dif z\, i \partial_z \del{\frac{1}{2} \envert{b_j(t)}^2 + \phi_j(z)}~, \end{equation} as the last term is a total derivative term and vanishes. The invariance under gauge transformations is especially important~\cite{Kroha2004} since, according to Elitzur's theorem, a non-zero saddle-point solution of the $b$ fields is only possible for gauge-invariant theories. In anticipation of the saddle-point approximation and considering the system to be infinite and translational invariant, the constraint field $\lambda$, as well as $b$, are site-independent, i.e., $\lambda_i = \lambda$ and $b_i = b$. The resulting projected partition function~\eqref{eq:projected-Z}, from now referred as $Z$, reads \begin{equation} Z = \int D\sbr{b^*, b} D\sbr{c^*, c} D\sbr{f^*, f} D\sbr{a^*, a} D\sbr{\lambda} e^{i S\sbr{b^*,b,c^*,c,f^*,f,a^*,a,\lambda}}~, \end{equation} where the Keldysh action $S\sbr{b^*,b,c^*,c,f^*,f,a^*,a,\lambda}$ is given by \begin{equation} \begin{split} S\sbr{b^*,b,c^*,c,f^*,f,a^*,a,\lambda} &= \int_\gamma \dif z\, \Bigg\{ \mathcal{N} b^*(z) \sbr{i\partial_z - \lambda(t)}b(z) + a^*(z) G_a^{{-1}}(z,z) a(z) \\ &+ \sum_{\bm{k}\sigma} \begin{bmatrix}c^*_{\bm{k}\sigma}(z) &f^*_{\bm{k}\sigma}(z)\end{bmatrix} \bm{G}^{-1}_{\bm{k}\sigma, \bm{k}\sigma}(z,z) \begin{bmatrix} c_{\bm{k}\sigma}(z) \\ f_{\bm{k}\sigma}(z) \end{bmatrix} + \mathcal{N} \lambda(z) \Bigg\}~, \end{split} \end{equation} where $\mathcal{N} \eqqcolon \sum_j 1$ denotes the number of lattice sites, $G_a^{^{-1}}(z,z')$ the inverse Green function of the photonic fields, and the fermionic inverse Green function $\bm{G}^{-1}_{\bm{k}\sigma,\bm{k}'\sigma'}(z,z')$ is given by \begin{equation} \bm{G}^{-1}_{\bm{k}\sigma,\bm{k}'\sigma'}(z,z') = \resizebox{0.8\linewidth}{!}{$\sbr{ \underbrace{\begin{bmatrix} i \partial_z - \varepsilon^c_{\bm{k}\sigma} &-V_0 b^*(z) \\ -V_0 b(z) &i \partial_z - \varepsilon^f_0 - \lambda(z) \end{bmatrix}}_{\bm{G}^{-1}_0} - i g_0 \underbrace{(a(z) - a^*(z)) \begin{bmatrix} 0 &b^*(z) \\ b(z) &0 \end{bmatrix}}_{\bm{A}} } \delta_\gamma(z,z') \delta_{\bm{k} \bm{k}'} \delta_{\sigma \sigma'}$}~. \end{equation} Note that while $\bm{G}^{-1}$ is block-diagonal in $\bm{k}$ and $\sigma$ space, the time differential operator $\partial_z$ is formally non-diagonal (\cref{sec:2pi}), as the continuous-time notation is but an abbreviation of the discrete path integral. For this reason, $\bm{G}^{-1}$ and its inverse are labelled by one momentum and spin index but by two time indices. This has important implications in the evaluation of the saddle-point equations, where, due to the definition of the $\delta$-distribution, $\bm{G}(z,z) = \bm{G}(z,z+0^+)$ which evaluates to $\bm{G}^<(t,t)$ for any $z\in\gamma$ in the Keldysh contour $\gamma$. By integrating out the Grassmann fields, the effective action $S_\mathrm{eff}$ reads \begin{equation} \begin{split} S_\mathrm{eff}\sbr{b^*,b,a^*,a,\lambda} &= -i \tr \log \sbr{- i \bm{G}^{-1}} \\ &+ \int_\gamma \dif z\, \cbr{\mathcal{N} b^*(z) \sbr{i\partial_z - \lambda(z)}b(z) + a^*(z) G_a^{{-1}}(z,z) a(z) + \mathcal{N} \lambda(z)}~. \end{split} \end{equation} Note that the first term traces over \textit{all} degrees of freedom of $\bm{G}$, including time. Furthermore, unless explicitly denoted by its discrete or continuous components, any boldface symbol should be considered to contain all respective degrees of freedom. For an invertible $\bm{G}^{-1}_0$, the logarithm can be expanded as \begin{equation} \tr \log -i \bm{G}^{-1} = \tr \log \sbr{-i \bm{G}_0^{-1} \del{\mathds{1} - i g_0 \bm{G}_0 \bm{A}}} = \tr \log -i \bm{G}_0^{-1} - \sum^{+\infty}_{n=1} \frac{\del{i g_0}^n}{n} \tr \del{\bm{G}_0 \bm{A} }^n~, \end{equation} with \begin{equation} \begin{split} \tr \del{\bm{G}_0 \bm{A}} &= \sum_{\bm{k}\bm{k}'} \sum_{\sigma \sigma'} \int_\gamma \dif z \dif z'\, \tr \cbr{\bm{G}_{0_{\bm{k}\sigma, \bm{k}'\sigma'}}(z,z') \bm{A}_{\bm{k}'\sigma', \bm{k}\sigma}(z',z)} % \\ % &= \sum_{\bm{k}\sigma} \int_\gamma \dif z \tr % \cbr{ % \begin{bmatrix} % G^{\del{cc}}_{0_{\bm{k}\sigma}}\del{z,z} & G^{\del{cf}}_{0_{\bm{k}\sigma}}\del{z,z} \\ % G^{\del{fc}}_{0_{\bm{k}\sigma}}\del{z,z} & G^{\del{ff}}_{0_{\bm{k}\sigma}}\del{z,z} % \end{bmatrix} % \begin{bmatrix} % 0 & b^*\del{z} \\ % b\del{z} & 0 % \end{bmatrix} % } \sbr{a^*\del{z} + a\del{z}} \\ &= \int_\gamma \dif z \underbrace{\sbr{\sum_{\bm{k}\sigma}\del{G^{\del{cf}}_{0_{\bm{k}\sigma}}(z,z) b(z) + G^{\del{fc}}_{0_{\bm{k}\sigma}}(z,z) b^*(z)}}}_{\alpha(z)} \sbr{a(z) - a(z)^*}~. \end{split} \end{equation} % \begin{equation} % \begin{split} % \tr \del{\bm{G}_0 \bm{A} \bm{G}_0 \bm{A}} = % &\int_\gamma \dif z \underbrace{\sum_{\bm{k}\sigma} % \sbr{b(z)^2 G^{\del{fc}}_{0_{\bm{k}\sigma}}(z)^2 + % b^*(z)^2 G^{\del{cf}}_{0_{\bm{k}\sigma}}(z)^2 + % \envert{b(z)}^2 G^{\del{cc}}_{0_{\bm{k}\sigma}}(z) G^{\del{ff}}_{0_{\bm{k}\sigma}}(z)}}_{\beta(z)} \\ % &\times \del{a(z)^* + a(z)}^2 % \end{split} % % \tr \del{\bm{G}_0 A \bm{G}_0 A} = \sum_{\bm{k}\sigma} % % \int_\gamma \dif z \tr_\mathcal{F} \cbr{ % % \begin{bmatrix} % % G^{\del{cc}}_{0_{\bm{k}\sigma}}\del{z} &G^{\del{cf}}_{0_{\bm{k}\sigma}}\del{z} \\ % % G^{\del{fc}}_{0_{\bm{k}\sigma}}\del{z} &G^{\del{ff}}_{0_{\bm{k}\sigma}}\del{z} % % \end{bmatrix} % % \begin{bmatrix} % % 0 &b\del{z} (a^*\del{z} + a\del{z}) \\ % % b^*\del{z} (a^*\del{z} + a\del{z}) &0 % % \end{bmatrix} % % \begin{bmatrix} % % G^{\del{cc}}_{0_{\bm{k}\sigma}}\del{z} &G^{\del{cf}}_{0_{\bm{k}\sigma}}\del{z} \\ % % G^{\del{fc}}_{0_{\bm{k}\sigma}}\del{z} &G^{\del{ff}}_{0_{\bm{k}\sigma}}\del{z} % % \end{bmatrix} % % \begin{bmatrix} % % 0 &b\del{z} (a^*\del{z} + a\del{z}) \\ % % b^*\del{z} (a^*\del{z} + a\del{z}) &0 % % \end{bmatrix} % % } % \end{equation} For a light-matter coupling far away from the ultrastrong regime $g_0 < 1$ and a small photon number $\ev{\hat a^\dagger(z)\hat a(z)} \lesssim 1$, it is justified to truncate the expansion at linear order. The photon fields are integrated by \textit{completing the square} \begin{equation} \begin{split} Z &= e^{\tr \log -i \bm{G}_0^{-1}} \int D\sbr{b^*, b} D\sbr{\lambda}e^{i \int_\gamma \dif z\, \cbr{\mathcal{N} b^*(z) \sbr{i\partial_z - \lambda(z)}b(z) + \mathcal{N} \lambda(z) }} \\ &\times \int D\sbr{a^*, a} e^{i \int_\gamma \dif z\, \cbr{ a^*(z) G_a^{{-1}}(z,z) a(z) + g_0 \,\alpha(z) i \sbr{a(z) - a(z)^*} + \mathcal{O}\del{g_0^2}}} \\ &\simeq \int D\sbr{b^*, b} D\sbr{\lambda} e^{i S_\mathrm{eff}\sbr{b^*,b,\lambda}}~, \end{split} \end{equation} where \begin{equation} \begin{split} S_\mathrm{eff}\sbr{b^*,b,\lambda} &= - i \tr \log -i \bm{G}_0^{-1} + i \tr \log -i \bm{G}_a^{{-1}} \\ &+ \int_\gamma \dif z\ \cbr{\mathcal{N} b^*(z) \sbr{i\partial_z - \lambda(z)}b(z) + \mathcal{N} \lambda(z)} - \int_\gamma \dif z\, \dif z'\, g_0^2\,\alpha(z) G_a(z,z') \alpha(z')~. \end{split} \end{equation} Note that the truncation of the $\log$ term at linear order neglects the renormalisation of the photon pulse arising from higher order coupling terms, e.g., \begin{equation} \begin{split} \tr \del{\bm{G}_0 \bm{A} \bm{G}_0 \bm{A}} &= \int_\gamma \dif z \dif z'\, \sum_{\bm{k}\sigma} \sbr{a^*(z) + a(z)}\sbr{a^*(z') + a(z')} \\ &\times\Bigg[ b(z) G^{\del{fc}}_{0_{\bm{k}\sigma}}\del{z,z'} G^{\del{fc}}_{0_{\bm{k}\sigma}}\del{z',z} b(z') + b(z) G^{\del{ff}}_{0_{\bm{k}\sigma}}\del{z,z'} G^{\del{cc}}_{0_{\bm{k}\sigma}}\del{z',z} b^*(z') \\ &+ b^*(z) G^{\del{cf}}_{0_{\bm{k}\sigma}}\del{z,z'} G^{\del{cf}}_{0_{\bm{k}\sigma}}\del{z',z} b^*(z') + b^*(z) G^{\del{cc}}_{0_{\bm{k}\sigma}}\del{z,z'} G^{\del{ff}}_{0_{\bm{k}\sigma}}\del{z',z} b(z') \Bigg]~. \end{split} \end{equation} Despite neglecting the build-up of coherence in the incident pulse, including such a term in the equations of motion is of considerable technical challenge and ultimately should not change the results much in this specific problem. \section{Non-equilibrium saddle-point equations} Apart from the truncation of the photon interactions, the partition function remains exact and accounts for most of the quantum effects. The partition function is now approximated by the saddle-point approximation \begin{equation} \label{eq:saddle-point-system} Z = \int D\sbr{b^*, b}D\sbr{\lambda} e^{i S_{\mathrm{eff}}\sbr{b^*,b,\lambda}} \approx e^{i S_\mathrm{eff}\sbr{b_0^*,b_0,\lambda_0}}~, \end{equation} where \begin{equation} \left.\fd{S_{\mathrm{eff}}\sbr{b^*,b,\lambda}}{\del{b^*(z), b(z), \lambda(z)}}\right|_{b^*(z) = b^*_0,\, b(z) = b_0,\, \lambda(z)=\lambda_0} = 0~. \end{equation} In non-equilibrium, the saddle-point condition generates a set of self-consistent equations describing the saddle-points' time evolution. Noting that $b^*_0$ and $b$ are related by complex conjugation, the saddle-point equations read \begin{subequations} \label{eq:saddle-point-equations} \begin{align} i \partial_t b_0(t) &= \lambda_0(t) b_0(t) + i \frac{1}{\mathcal{N}} \fd{\tr \log -i \bm{G}_0^{-1}}{b_0^*(t)} + g_0^2 \fd{\alpha(t)}{b_0^*(t)} \int_\gamma \dif z\,\sbr{G_a(t,z) + G_a(z,t)}\alpha(z)~, \\ 1 &= \envert{b_0(t)}^2 + i \frac{1}{\mathcal{N}} \fd{\tr \log -i \bm{G}_0^{-1}}{\lambda_0(t)} + g_0^2 \fd{\alpha(t)}{\lambda_0(t)} \int_\gamma \dif z\,\sbr{G_a(t,z) + G_a(z,t)}\alpha(z)~. \end{align} \end{subequations} Unlike the saddle-point equation for $b_0(t)$, the saddle-point equation for $\lambda_0(t)$ is not a differential equation but an algebraic one. The set of saddle-point equations is hence a differential-algebraic system of equations, where the equation for $\lambda_0(t)$ \textit{constrains} the possible trajectories of $b_0(t)$. Moreover, since the Hamiltonian is equally defined on both Keldysh contours, the fields follow the same equations of motion on each contour branch. This highlights the \textit{classical} nature of saddle-point equations, where the action on the forward path is cancelled by the one of the backward path~\cite{Kamenev2011_bosons}. Differentiating through the $\log$ term \begin{equation} \begin{split} \fd{\tr \log -i \bm{G}_0^{-1}}{x(t)} = \tr \bm{G}_0 \fd{ \bm{G}^{-1}_0}{x(t)} % &= \sum_{\bm{k}\bm{k}'} \sum_{\sigma \sigma'} \int_{\gamma} \dif z \dif z'\,\tr \bm{G}_{0_{\bm{k}\sigma,\bm{k}'\sigma'}}(z,z') \fd{\bm{G}^{-1}_{\bm{k}'\sigma',\bm{k}\sigma}(z',z)}{x(z)}\fd{x(z)}{x(t)} % \\ % &= = \sum_{\bm{k}\sigma} \tr \bm{G}_{0_{\bm{k}\sigma}}(t,t) \fd{\bm{G}^{-1}_{\bm{k}\sigma}(t,t)}{x(t)}~, \end{split} \end{equation} yields \begin{equation} \begin{split} \fd{\tr \log -i \bm{G}_0^{-1}}{b^*(t)} = % \sum_{\bm{k}\sigma} \,\tr \bm{G}_{0_{\bm{k}\sigma}}(t,t) \begin{bmatrix} % 0 &-V_0\\ % 0 &-\varepsilon^f_{\bm{k}\sigma} b(t) % \end{bmatrix} % = % -\sum_{\bm{k}\sigma} \sbr{V_0 G^{\del{fc}}_{0_{\bm{k}\sigma}}(t,t) +\varepsilon^f_{\bm{k}\sigma} b(t) G^{\del{ff}}_{0_{\bm{k}\sigma}}(t,t)}~, -\sum_{\bm{k}\sigma} V_0 G^{\del{fc}}_{0_{\bm{k}\sigma}}(t,t)~, \qquad \fd{\tr \log -i \bm{G}_0^{-1}}{\lambda(t)} = % \sum_{\bm{k}\sigma} \,\tr \bm{G}_{0_{\bm{k}\sigma}}(t,t) \begin{bmatrix} % 0 &0\\ % 0 &-1 % \end{bmatrix} % = -\sum_{\bm{k}\sigma} G^{\del{ff}}_{0_{\bm{k}\sigma}}(t,t)~. \end{split} \end{equation} Differentiating through the $\alpha$ term, % \begin{equation} % \begin{split} % \fd{}{x(t)} \int_{\gamma} \dif z\, g_0^2 G_a(z,z) \alpha(z)^2 % &= \int_{\gamma} \dif z\, 2 g_0^2 G_a(z,z) \alpha(z) \fd{\alpha(z)}{x(z)} \fd{x(z)}{x(t)} % % \\ % = 2 g_0^2 G_a(t,t) \alpha(t) \fd{\alpha(t)}{x(t)}~, % \end{split} % \end{equation} given $\fd{}{x}\bm{U} = \bm{U} \fd{}{x}\del{-\bm{U}^{-1}} \bm{U}$, yields \begin{subequations} \begin{align} \fd{\alpha(t)}{b^*(t)} % &= \sum_{\bm{k}\sigma}\Bigg[G^{\del{fc}}_{0_{\bm{k}\sigma}}(t,t) + \del{V_0 G^{\del{cc}}_{0_{\bm{k}\sigma}}(t,t) % + \varepsilon^f_{\bm{k}\sigma} b(t) G^{\del{cf}}_{0_{\bm{k}\sigma}}(t,t)}G^{\del{ff}}_{0_{\bm{k}\sigma}}(t,t) b(t) % &\phantom{= \sum_{\bm{k}\sigma}} + \del{V_0 G^{\del{fc}}_{0_{\bm{k}\sigma}}(t,t) + \varepsilon^f_{\bm{k}\sigma} b(t) G^{\del{ff}}_{0_{\bm{k}\sigma}}(t,t) }G^{\del{fc}}_{0_{\bm{k}\sigma}}(t,t) b^*(t)\Bigg]~, &= \sum_{\bm{k}\sigma}\sbr{G^{\del{fc}}_{0_{\bm{k}\sigma}}(t,t) + V_0 \del{G^{\del{cc}}_{0_{\bm{k}\sigma}}(t,t) G^{\del{ff}}_{0_{\bm{k}\sigma}}(t,t) b(t) + G^{\del{fc}}_{0_{\bm{k}\sigma}}(t,t) G^{\del{fc}}_{0_{\bm{k}\sigma}}(t,t) b^*(t)}}~, \\ \fd{\alpha(t)}{\lambda(t)} &= \sum_{\bm{k}\sigma} \sbr{G^{\del{cf}}_{0_{\bm{k}\sigma}}(t,t)G^{\del{ff}}_{0_{\bm{k}\sigma}}(t,t)b(t) + G^{\del{ff}}_{0_{\bm{k}\sigma}}(t,t)G^{\del{fc}}_{0_{\bm{k}\sigma}}(t,t)b^*(t)}~. \end{align} \end{subequations} Furthermore the approximation \begin{equation} \label{eq:sp-pulse-approximation} \int_\gamma \dif z\,\sbr{G_a(t,z) + G_a(z,t)}\alpha(z) \approx \sbr{G_a(t,t) + G_a(t,t)}\alpha(t)~, \end{equation} is employed, which holds for short pulses, turning the non-equilibrium saddle-point equations~\eqref{eq:saddle-point-equations} of motion into a Markovian set of equations. \section{Initial conditions} The system is in thermal equilibrium before the pulse's arrival. Due to~\eqref{eq:saddle-point-system} being an effective non-interacting problem and the saddle-points fields 1-point functions, unlike in interacting problems of 2-point functions, the problem can be separated into two disjoint parts. First, solving the saddle-point equations on the Matsubara branch of the Konstantinov-Perel' contour (\cref{fig:konstantinov-contour}). Second, solving the time-dependent saddle-point equations on the Keldysh branches of the Konstantinov-Perel' contour using the thermal solutions as initial conditions. This is also highlighted by the purely Markovian form of the saddle-point set of equations~\eqref{eq:saddle-point-system}. Considering the light pulse to vanish in the Matsubara branch and the fields to be static in imaginary time, the saddle-point equations for the system in thermal equilibrium read \begin{subequations} \begin{align} 0 &= \lambda_0 b_0 - \frac{1}{\beta} \frac{1}{\mathcal{N}} \fd{\tr \log \bm{G}_0^{-1}}{{b_0^*}}~, \\ 1 &= \envert{b_0}^2 - \frac{1}{\beta} \frac{1}{\mathcal{N}} \fd{\tr \log \bm{G}_0^{-1}}{\lambda_0}~, \end{align} \end{subequations} where $\bm{G}_0^{-1} = \mathds{1}\partial_\tau + \bm{H}$. Since $\bm{G}_0^{-1}$ is block-diagonal in imaginary time, \begin{equation} \tr \log \bm{G}_0^{-1} = \log \det \bm{G}_0^{-1} = \sum_{\bm{k}\sigma} \log \det \del{\mathds{1}\partial_\tau + \bm{H}_{\bm{k}\sigma}} = \sum_{\bm{k}\sigma} \sum_j \log \det \del{\partial_\tau + \omega_j}~, \end{equation} where $\omega_{\bm{k}\sigma}^j$ are the eigenvalues of $\bm{H}_{\bm{k}\sigma}$. For Grassmann fields, it can be shown~\cite{Kamenev2011_fermions} that $\det \del{\partial_\tau + \omega_{\bm{k}\sigma}^j} = 1 + \rho(\omega_{\bm{k}\sigma}^j)$, where $\rho(\omega)$ is the Boltzmann factor $e^{-\beta \omega}$. This non-trivial result can be attributed to the boundary terms from constructing the path integral (\cref{sec:path-integral}). In the limit of vanishing temperature \begin{equation} \lim_{\beta \to \infty} -\frac{1}{\beta}\log\del{1 + e^{-\beta \omega}} = \Theta(-\omega) \omega~, \end{equation} the saddle-point equations read \begin{subequations} \begin{align} 0 &= \lambda_0 b_0 + \sum_{\bm{k}\sigma} \sum_{\omega_{\bm{k}\sigma} < 0} \fd{}{{b_0^*}} \omega_{\bm{k}\sigma}^j~, \\ 1 &= \envert{b_0}^2 + \sum_{\bm{k} \sigma} \sum_{\omega_{\bm{k}\sigma} < 0} \fd{}{\lambda_0} \omega_{\bm{k}\sigma}^j~. \end{align} \end{subequations} This unusual form of the saddle-point equations is due to the functional derivatives being taken through a $\log \det$ instead of the typical $\tr \log$, which would result in a set of saddle-point equations dependent on 2-point averages instead of derivatives of eigenvalues of the Hamiltonian. The power hidden in the \enquote{unsual} derivatives is that the saddle-point equations can be obtained through automatic differentiation by specifying a single Hamiltonian instead of laborious and error-prone expansions of 2-point averages. Furthermore, considering a $d$-dimensions lattice volume $\mathcal{V}$, the sum over the lattice momenta takes the form \begin{equation} \lim_{\mathcal{N} \to \infty} \frac{1}{\mathcal{N}} \sum_{\bm{k}} = \lim_{\mathcal{N} \to \infty} \frac{1}{\mathcal{N} \Delta k} \sum_{\bm{k}} \Delta k = \lim_{\mathcal{N},\mathcal{V} \to \infty} \frac{\frac{\mathcal{V}}{(2\pi)^d}}{\mathcal{N}} \sum_{\bm{k}} \Delta k = \frac{\mathcal{V}}{\mathcal{N}} \int \frac{\dif\bm{k}}{(2\pi)^d} = \int_{\mathrm{FBZ}} \frac{\dif\bm{k}}{(2\pi)^d}~, \end{equation} where in the last equality the momenta were rescaled as $a^d \bm{k} \to \bm{k}$, where $a$ is the lattice constant and $\frac{\mathcal{V}}{\mathcal{N} a^d} \stackrel{!}{=} 1$. The last integral is taken over the first Brillouin zone with a unit lattice constant. \section{Non-equilibrium saddle-point solutions} Since the fluctuations of the boson field are expected to occur at frequencies associated with $\varepsilon^f_0$~\cite{Coleman_2015}, the constraint field is shifted as $\lambda_0(t) \to \lambda_0(t) - \varepsilon^f_0$ such that $\lambda_0(t)$ is small and directly corresponds to the Kondo coherence temperature~\eqref{eq:kondo-temp-sp}. Furthermore, the saddle-point approximation freezes out the $b$ field, which encodes the state of an \textit{empty} $f$-site. As a result, the $f$-electron degrees of freedom left in the system are of a singly-occupied $f$-site, and it is expected that $f$ describes similar dynamics to the Kondo Hamiltonian~\cite{Schrieffer_1966}, that is, spin-fluctuations at around the Fermi energy. Moreover, due to the approximation~\eqref{eq:sp-pulse-approximation}, the system of equations cannot \enquote{see} the carrier frequency of the photon pulse and the system is driven only by the intensity of the photon pulse. For these reasons, interpreting Kondo saddle-point results is quite challenging and serves only as a starting point for inspecting the physics of the problem. In order to solve the non-equilibrium saddle-point equations~\eqref{eq:saddle-point-equations}, solutions to the non-interacting Green functions of matter fields can be obtained by \begin{equation} \bm{G}_0(t,t') = -i \, \mathds{T}\sbr{e^{-i \int_{t_0}^t \dif \bar{t}\,{\bm{H}}(\bar{t})}} \sbr{\Theta_\gamma(t,t')\mathds{1} + \xi \bm{n}} \bar{\mathds{T}}\sbr{e^{+i \int_{t_0}^{t'} \dif \bar{t}\,{\bm{H}}(\bar{t})}}~, \end{equation} where $\bm{n}$ is the initial occupation of the fields, at time $t=t_0$. However, at equal times, it suffices to calculate \begin{equation} i \partial_t \bm{G}_0(t,t) = \sbr{\bm{H}(t), \bm{G}_0(t,t)}~, \end{equation} subject to $\bm{G}_0(t_0,t_0') = \sbr{\Theta_\gamma(t_0,t_0')\mathds{1} + \xi \bm{n}}$. For the following analysis, the system parameters are $\varepsilon_0^f = -0.35$, $g_0 = 0.045$ and $V_0=0.3$, in units of the conduction electron hopping $\upsilon$, in a Bethe lattice with infinite connectivity (\cref{sec:dmft}), for congruency with later results. The set of differential-algebraic equations~\eqref{eq:saddle-point-equations} is very unstable and solved via an \textit{implicit Euler} scheme~\cite{Rackauckas2017}, with $\texttt{rtol} \sim 10^{-5}$ and $\texttt{atol} \sim 10^{-5}$. Until the system is perturbed by the external pulse, the saddle-point fields do not change in time since the system is in thermal equilibrium. Note that \textit{thermal equilibrium} at the level of 1-point fields is a very loose term since the saddle-point fields only capture field averages and hence any $n$-point correlation that could encode thermal distributions -- and fulfil, e.g. fluctuation-dissipation relations (\cref{sec:fdr}) -- are not available. Note, however, that thermal distributions can be encoded in the initial occupation of the matter fields $\bm{G}_0(t_0,t_0)$ but are \textit{fixed} by the initial conditions since their time-evolution is dictated by a \textit{non-interacting} Hamiltonian. In~\cref{fig:saddle-points-hybridisation}, the photon-pulse intensity is kept at a reasonable value $n_a=1.0$, which is expected to drive the system out of equilibrium, however, without collapsing the Kondo ground-state. Here, it is evident how a larger Kondo coherence temperature, set by the hybridisation strength $V_0$, generates faster oscillations in the relaxation of $\lambda_0$ towards its thermal-equilibrium value. Despite the magnitude of the perturbation to the saddle-point fields being similar whether the pulse is short (\cref{fig:saddle-points-hybridisation-1}) or long (\cref{fig:saddle-points-hybridisation-2}), two slightly different relaxation scenarios arise. For large Kondo temperatures, under a short pulse, the system is driven out of equilibrium and relaxes back to its ground state with intrinsic fast oscillations. However, for pulses longer than the period of the oscillations, the pulse drive is effectively adiabatic, resulting in a smooth, oscillation-free relaxation. For systems with small Kondo temperatures, both short and long pulses have shorter durations than the oscillation period, resulting in a long and slow oscillatory relaxation to $\lambda_0(t_0)$. \begin{figure}[!htb] \centering \begin{subfigure}[b]{0.49\linewidth} \centering \includesvg[width=\linewidth]{figs/sp0} \caption{A gaussian pulse with width $\Omega_0 = 0.5$ and pulse maximum at $t\,\upsilon= 20$, shaded in grey.} \label{fig:saddle-points-hybridisation-1} \end{subfigure} \hfill \begin{subfigure}[b]{0.49\linewidth} \centering \includesvg[width=\linewidth]{figs/sp1} \caption{A gaussian pulse with width $\Omega_0 = 0.1$ and pulse maximum at $t\,\upsilon = 100$, shaded in grey.} \label{fig:saddle-points-hybridisation-2} \end{subfigure} \caption{Time dependence of $\lambda_0(t)$ for different values of the hybridisation strength $V_0$, with maximum pulse intensity $n_a=1$.} \label{fig:saddle-points-hybridisation} \end{figure} In~\cref{fig:saddle-points-pulse}, it is revealed that the system responds rather strongly to the pulse's intensity. As expected at the saddle-point level, the recovery of the $\lambda_0$ field is modulated by a single frequency, independent of the pulse's intensity. For values of intermediate pulse intensity $n_a \gtrsim 0.5$, the value of $\lambda_0(t)$ changes significantly, and for $n_a = 8.0$, the value of $\lambda_0(t)$ increases to almost double its equilibrium value. Moreover, for long pulses and a large photon-pulse intensity (\cref{fig:saddle-points-pulse-2}) the system relaxes to a different value of $\lambda_0$. This is indeed surprising since the saddle-point equations~\eqref{eq:kondo-temp-sp} in thermal equilibrium are thought to have a unique solution within the physical regime ($b_0^2 \le 1$ and $\lambda_0 \sim 0$), which should imply that $\lim_{t\to\infty} \lambda_0(t) = \lambda_0(t_0)$. Despite not appearing to be a numerical problem, as running the differential-equation solver with different tolerances yields similar results, this cannot be ruled out since the tolerances cannot be set to lower values without incurring numerical instabilities. In this regime, however, the Kondo ground state is expected to collapse due to a strong interaction with the photon pulse. Since that is not verified, in the form of no solution for $\lambda_0(t)$, there could be a dynamical breakdown of the saddle-point approximation. That is, for time-dependent problems, the assumption that the fluctuations of the $b$ fields are small most likely does not hold, and further investigation into criteria for the validity of non-equilibrium saddle points is warranted. \begin{figure}[!htb] \centering \begin{subfigure}[b]{0.49\linewidth} \centering \includesvg[width=\linewidth]{figs/sp2} \caption{A gaussian pulse with width $\Omega_0 = 0.5$ and pulse maximum at $t\,\upsilon = 20$, shaded in grey.} \label{fig:saddle-points-pulse-1} \end{subfigure} \hfill \begin{subfigure}[b]{0.49\linewidth} \centering \includesvg[width=\linewidth]{figs/sp3} \caption{A gaussian pulse with width $\Omega_0 = 0.1$ and pulse maximum at $t\,\upsilon = 100$, shaded in grey.} \label{fig:saddle-points-pulse-2} \end{subfigure} \caption{Time dependence of $\lambda_0(t)$ for different values of the maximum photon intensity $n_a$.} \label{fig:saddle-points-pulse} \end{figure}