\chapter{Auxiliary particles in non-equilibrium}
\label{sec:aux-particles}

The introduction of auxiliary particles dates back to the foundational work of Abrikosov~\cite{Abrikosov_1965}, which introduced a faithful representation of spin-$1/2$ operators in terms of auxiliary fermionic operators. Unlike the former, which obey the SU$(2)$ algebra and are incompatible with Wick's theorem, the latter obey canonical commutation relations and fully support standard diagrammatic techniques and saddle-point theories. Especially relevant to this thesis are the posterior works~\cite{Read_1983, Coleman_1984} which introduced auxiliary bosonic particles to study \textit{impurity} fermionic particles $f$ subject to infinitely strong Coulomb repulsion. In such systems, the Coulomb interaction term can be dropped out of the Hamiltonian provided that the operator constraint
\begin{equation}
    \sum_\sigma \hat f^\dagger_\sigma \hat f_\sigma \le \hat{\mathds{1}}~,
\end{equation}
holds, where $\sigma$ denotes the local degrees of freedom of $f$, such as spin. Such non-holonomic constraints are of difficult analytic implementation but can be turned into holonomic constraints via the introduction of an auxiliary bosonic particle
\begin{equation}
\label{eq:auxiliary-particle-map}
    \hat f_\sigma^\dagger \to \hat f_\sigma^\dagger \hat b~,
\end{equation}
while still retaining the original fermionic character
\begin{equation}
    \cbr{\hat f_\sigma, \hat f^\dagger_{\sigma'}} \to \cbr{\hat b^\dagger \hat f_\sigma, \hat f^\dagger_{\sigma'} \hat b} \stackrel{!}{=} \delta_{\sigma\sigma'}~,
\end{equation}
if the operators are subject to the operator constraint
\begin{equation}
\label{eq:auxiliary-particle-constraint}
    \hat Q \coloneqq \hat b^\dagger \hat b + \sum_\sigma \hat f^\dagger_\sigma \hat f_\sigma = \hat{\mathds{1}}~.
\end{equation}
Auxiliary-particle mappings such as~\eqref{eq:auxiliary-particle-map} are generally not unique -- in fact, for certain problems, there may be an infinite number of possible equivalent representations. However, this equivalence between representations can break down when approximations -- such as diagrammatic truncations -- are employed~\cite{Bickers_1987}. Since virtually any analytic calculation requires approximations, it is important that the auxiliary-particle mapping has a somewhat physical foundation. Only then can it be possible to infer whether the resulting approximation captures the relevant physical processes and is valid. These questions will be later explored regarding the physics resulting from saddle-point (\cref{sec:thz-mf}) or diagrammatic truncations (\cref{sec:thz-nca}) of theories relying on auxiliary particles.

Even though the use of auxiliary particles in non-equilibrium settings dates back to works in the early 1990s~\cite{Langreth_1991, Wingreen_1994}, a proper extension to all contour-ordered Green function components was yet lacking. Despite modern attempts~\cite{Eckstein_2010} at extending the formalism to the Konstantinov-Perel' contour (\cref{fig:konstantinov-contour}), the auxiliary-particle formalism can be made much more transparent and complete, especially regarding thermal Green functions, which are necessary as initial conditions in the study of systems driven out of thermal equilibrium (\cref{sec:neqft-components}).

\section{Abrikosov's projection}
\label{sec:ap-zeta-scaling}
The addition of auxiliary particles enlarges the original -- or physical -- Hilbert space with unphysical degrees of freedom which must be \textit{projected} out. This is carried out by enforcing the constraint~\eqref{eq:auxiliary-particle-constraint} via a projection technique~\cite{Zawadowski_1969}, which generalized Abrikosov's~\cite{Abrikosov_1965} original proposal. Consider the grand-canonical ensemble related to the auxiliary-particle charge $Q$ and associated chemical potential $\lambda$ or fugacity $\zeta = e^{-\beta\lambda}$. By construction -- or virtue of the Hamiltonian, each auxiliary-particle Fock space is disjoint, $\sbr{\hat{Q}, \hat{H}} = 0$, and hence the grand-canonical partition $Z_\zeta$ can be written as a sum over \textit{canonical} partition functions $Z_Q \coloneqq \tr_{\mathcal{H_{Q}}} \hat{\rho}_0$,
\begin{equation}
Z_\zeta = \tr\sbr{\zeta^{\hat{Q}-\hat{\mathds{1}}}\,\hat{\rho}_0} = 
\sum_{Q=0}^{\infty}\, \zeta^{Q-1}\,\tr_{\mathcal{H}_Q} \hat{\rho}_0 = 
\zeta^{-1} \sbr{Z_0 + \zeta Z_1 + \mathcal{O}(\zeta^2)}~,
\end{equation}
where $\mathcal{H}_Q$ denotes the Hilbert space of the system in the auxiliary-particle Fock sector with charge $Q$. The grand-canonical ensemble average of some operator $\hat{O}$
\begin{equation}
\label{eq:grand-canonical-ensemble-average}
\ev{\hat{O}}_\zeta
% = 
% \frac{\tr\sbr{\zeta^{\hat{Q}}\,\hat{\rho}_0\,\hat{O}}}{Z_{\mathcal{G}}(\zeta)} 
\coloneqq
\frac{\sum_{Q=0}^{\infty}\, \zeta^{Q-1}\,\tr_{\,\mathcal{H}_Q}\sbr{\hat{\rho}_0\,\hat{O}}}{Z_\zeta}~,
\end{equation}
can be related to the \textit{physical} (constrained to the $\mathcal{H}_{Q=1}$ sector) canonical ensemble average via the fugacity limit
\begin{equation}
\label{eq:projection-Q}
\begin{split}
\ev{\hat{O}}_\mathrm{physical}
\coloneqq
\frac{\tr_{\mathcal{H}_1}\sbr{\hat{\rho}_0\, \hat{O}}}{\tr_{\mathcal{H}_1}\hat{\rho}_0}
% \\
&\stackrel{\mathrm{!}}{\equiv} 
\lim_{\zeta\to0} \frac{\ev{\hat{O}\, \hat{Q}}_\zeta}{\ev{\hat{Q}}_\zeta} 
= 
\lim_{\zeta\to0} \frac{ \cancelto{0}{\tr_{\mathcal{H}_0}\sbr{\hat{\rho}_0\, \hat{O}\,\hat{Q}}} + \zeta \tr_{\mathcal{H}_1}\sbr{\hat{\rho}_0\, \hat{O}\, \hat{Q}} + \mathcal{O}(\zeta^2)}{\cancelto{0}{\tr_{\mathcal{H}_0}\sbr{\hat{\rho}_0\, \hat{Q}}} + \zeta \tr_{\mathcal{H}_1}\sbr{\hat{\rho}_0\, \hat{Q}} + \mathcal{O}(\zeta^2)}~.
% =\lim_{\zeta\to0} \frac{\tr_{\mathcal{H}_1}\sbr{\hat{\rho}_0\, \hat{O}} + \mathcal{O}(\zeta)}{\tr_{\mathcal{H}_1}\hat{\rho}_0 + \mathcal{O}(\zeta)}~.
\end{split}
\end{equation}
This equation is at the heart of the formalism and relates a canonical ensemble average in a constrained Hilbert space with an unconstrained grand canonical one, taken in some particular limit. All its constituents play an essential role -- for example, $\ev{Q}_\zeta$ is not \textit{just} another constant and will play a critical role in determining leading-order diagrams.

Note that an impurity operator $\hat O_\mathrm{imp}$ acting on the physical Hilbert space must be transformed to an operator acting on the enlarged Hilbert space, e.g., the prototypical $\hat O_\mathrm{imp} = \hat f^\dagger \hat f \equiv \hat f^\dagger \hat b \hat b^\dagger \hat f$. Such \textit{normal-ordered} auxiliary-particle operators annihilate the $Q=0$ sector by construction due to the presence of an auxiliary-particle annihilation operator on the right. Because of this, $\hat Q$ can be dropped out of the numerator -- since its purpose was to annihilate $\mathcal{H}_{Q=0}$, and the projection~\eqref{eq:projection-Q} reduces to
\begin{equation}
\label{eq:projection}
    \ev{\hat{O}_\mathrm{imp}}_\mathrm{physical} \equiv \lim_{\zeta\to0} \frac{\ev{\hat{O}_\mathrm{imp}}_\zeta}{\ev{\hat{Q}}_\zeta}~.
\end{equation}
However, the same does \textit{not} apply to non-impurity operators, and the projection~\eqref{eq:projection-Q} must be considered. Now that it is known how to orchestrate some grand-canonical ensemble averages to obtain a constrained canonical ensemble average, all remaining is to determine the expressions or equations of motion for said grand-canonical ensemble averages.

\section{\texorpdfstring{$\zeta$}{Zeta}-scaling of contour-ordered Green functions}
\label{sec:project-gf}

The projection~\eqref{eq:projection-Q} solely requires grand-canonical ensemble averages in the $\zeta\to0$ limit. For this reason, any expression that can be expanded in powers of $\zeta$ will be taken only to the leading order. The contour-ordered Green function~\eqref{eq:gf-time-decomposition} of auxiliary particles
\begin{equation}
\label{eq:gf-time-decomposition-projected}
\begin{split}
    G_\zeta(z, z') &= \Theta_\gamma(z, z') G_\zeta^>(z, z') + \Theta_\gamma(z', z) G_\zeta^<(z, z')
    \\
    &=\Theta_\gamma(z, z') \resizebox{0.33\linewidth}{!}{$\frac{\tr_{\mathcal{H}_0}\sbr{\hat{\rho}_0\, (-i)\hat f(z) \hat f^\dagger(z')} + \mathcal{O}(\zeta)}{Z_\zeta}$} + \Theta_\gamma(z', z) \resizebox{0.33\linewidth}{!}{$\frac{\zeta\tr_{\mathcal{H}_1}\sbr{\hat{\rho}_0\, (-i \xi)\hat f^\dagger(z')\hat f(z)} + \mathcal{O}(\zeta^2)}{Z_\zeta}$}
    \\
    &=\Theta_\gamma(z, z') \sbr{G_{\zeta,\mathcal{H}_0}^>(z, z') +\mathcal{O}(\zeta)} + \Theta_\gamma(z', z) \sbr{\zeta G_{\zeta,\mathcal{H}_1}^<(z, z') + \mathcal{O}(\zeta^2)}~,
\end{split}
\end{equation}
now displays a key feature of the grand-canonical formalism: the greater and lesser components have \textit{different} $\zeta$-scaling. In the limit $\zeta\to0$, the greater function is only traced in the $Q=0$ sector, and as a result, it does not carry information about particle occupation and hence only spectral information. Similarly, spectral functions such as the retarded Green function reduce to
\begin{equation}
\label{eq:retarded-decomp}
    \lim_{\zeta\to0}G^R_\zeta(t,t') = \lim_{\zeta\to0}\Theta(t-t') \sbr{G^>_\zeta(t,t') - G^<_\zeta(t,t')} = \lim_{\zeta\to0}\Theta(t-t') G^>_\zeta(t,t')~.
\end{equation}
Note that this unusual result is not some esoteric property of auxiliary particles but simply a formal \textit{consequence} of taking the grand-canonical ensemble average in the $\zeta\to0$ limit, a requirement of the projection~\eqref{eq:projection-Q} to the physical Hilbert space. Note, however, that non-auxiliary-particle Green functions always scale as $\mathcal{O}(1)$ and hence do not fulfil such relations.

\subsection{Kadanoff-Baym equations}
\label{sec:project-kb}

The rather unusual $\zeta\to0$ limit results in a modified set of Langreth rules~\eqref{eq:countour-conv} that directly enter the Kadanoff-Baym equations of motion~\eqref{eq:kb} for $\zeta$-scaled contour-ordered Green functions of auxiliary particles. This follows from keeping only the leading $\zeta$ terms of the $\zeta$-decomposition of contour-ordered Green functions~\eqref{eq:gf-time-decomposition-projected} and self-energies
\begin{equation}
\begin{split}
\int_\gamma \dif \bar{z} \,\Sigma_\zeta(z, \bar{z}) G_\zeta(\bar{z}, z') &= 
\Theta_\gamma(z,z')\underbrace{\int_{z'}^z \dif \bar{z}\,\Sigma_\zeta^>(z,\bar{z})G_\zeta^>(\bar{z},z')}_{\textrm{at least }\mathcal{O}(1)} +
\Theta_\gamma(z',z)\underbrace{\int_{z}^{z'} \dif \bar{z}\,\Sigma^<_\zeta(z,\bar{z})G^<_\zeta(\bar{z},z')}_{\textrm{at least }\mathcal{O}(\zeta^2)}
\\
&+
\underbrace{\int_{t_0}^{\min_\gamma(z,z')} \dif \bar{z}\,\Sigma_\zeta^>(z,\bar{z})G_\zeta^<(\bar{z},z')
+\int_{\max_\gamma(z,z')}^{t_0-i\beta} \dif \bar{z}\,\Sigma^<_\zeta(z,\bar{z})G^>_\zeta(\bar{z},z')}_{\textrm{at least }\mathcal{O}(\zeta)}~.
\end{split}
\end{equation}
The following analyses will consider that the auxiliary-particle self-energies have at least the \textit{same} $\zeta$-scaling as the contour-ordered Green functions. Note, however, that the $\zeta$-scaling of the self-energy of non-auxiliary particles is model-dependent and cannot be inferred a priori.

For a Green function with real-time arguments $(z,z') \to (t, t')$, the equations of motion reduce to
\begin{subequations}
\label{eq:kb-projected}
\begin{align}
\lim_{\zeta\to0}\sbr{i \partial_t - h_0(t)} G^>_\zeta(t, t') &=
\lim_{\zeta\to0}\cbr{\int_{t'}^t \dif \bar{t}\, 
\Sigma^>_\zeta(t, \bar{t}) G^>_\zeta(\bar{t}, t') + \mathcal{O}(\zeta)}
\\
\label{eq:kb-projected-lesser}
\begin{split}
\lim_{\zeta\to0}\sbr{i \partial_t - h_0(t)} G^<_\zeta(t, t') &= 
\lim_{\zeta\to0} \Bigg\{\int_{t_0}^{t} \dif \bar{t}\, 
\Sigma^>_\zeta(t, \bar{t}) G^<_\zeta(\bar{t}, t')
-\int_{t_0}^{t'} \dif \bar{t}\,
\Sigma^<_\zeta(t, \bar{t}) G^>_\zeta(\bar{t}, t')
\\
&\phantom{=}
+ \int_{t_0}^{t_0 - i \beta} \dif \bar{\tau}\, \Sigma^\rceil_\zeta(t, \bar{\tau}) G^\lceil_\zeta(\bar{\tau}, t')
+ \mathcal{O}(\zeta^2)\Bigg\}~.
\end{split}
\end{align}
\end{subequations}

\subsection{Gaussian initial conditions}

For an arbitrary non-thermal initial density matrix, the thermal-branch integral of~\eqref{eq:kb-projected-lesser} is zero -- since, by definition, it only arises from the extension of the Keldysh contour to imaginary times for a thermal initial density matrix (\cref{sec:kp-contour}). Despite the $\zeta$-term appearing to be temperature dependent -- which is a useful construction for \cref{sec:lambda-scaling}, in the absence of a definition of temperature, as in arbitrary density matrices, this detail can be ignored as long as the proper $\zeta$ limit is considered. For a Gaussian initial density matrix $\hat \rho_0 = \kappa^{\hat Q}$, the boundary term in the path integral quantum field theory (\cref{sec:path-integral}) that encodes the initial distribution is instead the expectation value of $\lim_{\zeta\to0} \zeta^{\hat Q} \kappa^{\hat Q}$ -- due to a different definition of the ensemble average (cf.~\eqref{eq:grand-canonical-ensemble-average}), and the initial conditions read~\cite{Kamenev2011_bosons, Kamenev2011_fermions, lappe2021non}
\begin{subequations}
\begin{align}
    \lim_{\zeta\to0} G_\zeta^>(t_0, t_0) &= \lim_{\zeta\to0} \sbr{-i \del{1 + \xi \frac{\zeta \kappa}{1 - \zeta \kappa}}} = -i
    \\
    \lim_{\zeta\to0} G_\zeta^<(t_0, t_0) &= \lim_{\zeta\to0} \sbr{- i \xi \frac{\zeta \kappa}{1 - \zeta \kappa}}~.
\end{align}
\end{subequations}

\section{\texorpdfstring{$\lambda$}{Lambda}-scaling of thermal Green functions}
\label{sec:lambda-scaling}
The application of the decomposition~\eqref{eq:gf-time-decomposition-projected} equation of motion of the Matsubara Green function $\mathcal{G}_\zeta(\tau, \tau')$ truncates the thermal-branch integral in the $\zeta\to0$ limit~\cite{Eckstein_2010}. However, due to the convenient analytic properties of Matsubara frequencies, it is worthwhile to keep the full integral. For that, instead of treating the grand-canonical ensemble average as a $\zeta$-expansion~\eqref{eq:grand-canonical-ensemble-average}, consider the equivalent $\lambda$-expansion for an initial thermal density matrix
\begin{equation}
    \ev{\hat O}_\zeta \equiv \ev{\hat O}_\lambda \coloneqq \frac{\tr \sbr{e^{-\beta \del{\hat{\mathcal{H}} + \lambda \del{\hat Q - \hat{\mathds{1}}}}}\, \hat O}}{\tr \sbr{e^{-\beta \del{\hat{\mathcal{H}} + \lambda \del{\hat Q - \hat{\mathds{1}}}}}}}~.
\end{equation}
The auxiliary chemical potential enters directly in the non-interacting part of the auxiliary-particle Hamiltonian (cf.~\eqref{eq:hamiltonian-kp-contour}), and the right-hand-side of the equation of motion follows standard Langreth rules~\eqref{eq:countour-conv}
\begin{equation}
    \lim_{\lambda\to\infty} \sbr{-\partial_{\tau} - h_0 - \lambda}\mathcal{G}_\lambda(\tau, \tau') = i\delta(\tau-\tau') - i \lim_{\lambda\to\infty} \cbr{\int_{0}^{\beta}\dif \bar{\tau}\, \Sigma_\lambda(\tau, \bar{\tau}) \mathcal{G}_\lambda(\bar{\tau}, \tau')}~.
\end{equation}
Note that $\Sigma_\lambda$ here denotes the Matsubara self-energy $\Sigma_\lambda(\tau, \tau') \coloneqq \Sigma_\lambda(t_0 - i \tau, t_0 - i \tau')$. Leveraging the translational symmetry $\mathcal{G}_\lambda(\tau, \tau') \equiv \mathcal{G}_\lambda(\tau - \tau')$ and the (anti)periodicity $\mathcal{G}_\lambda(\tau + \beta, \tau') = \xi\, \mathcal{G}_\lambda(\tau, \tau')$ of the Matsubara functions -- which follows directly from their definition~\eqref{eq:contour-gf}, the imaginary-time integrals can be replaced by products in Matsubara frequencies $i \omega_n$
\begin{equation}
\label{eq:fourier-matsubara}
% \begin{split}
\mathcal{G}_\lambda(\tau) = \frac{i}{\beta}\sum_{i \omega_n} e^{-i \omega_n \tau} \mathcal{G}_\lambda(i \omega_n)
\qquad
\mathcal{G}_\lambda(i \omega_n) =
-i \int_0^\beta \dif \tau e^{i \omega_n \tau} \mathcal{G}_\lambda(\tau)~,
% \end{split}
\end{equation}
where $i \omega_n = (2n)\pi/\beta$ or $i \omega_n = (2n + 1)\pi/\beta$ for bosonic or fermionic Green functions, respectively. Thus,
\begin{equation}
\label{eq:matsubara-gf-aux-particles}
    \lim_{\zeta\to0} \mathcal{G}_\zeta(i \omega_n) \equiv \lim_{\lambda\to\infty} \mathcal{G}_\lambda(i \omega_n) = \lim_{\lambda\to\infty} \sbr{i \omega_n - h_0 - \lambda - \Sigma_\lambda(i \omega_n)}^{-1}~.
\end{equation}
The mixed Green functions have similar properties in the imaginary time and hence can also be expressed in terms of Matsubara frequencies
\begin{equation}
\begin{split}
\lim_{\lambda\to\infty}\sbr{i \partial_t - h(t)} G^\rceil_\lambda(t, i \omega_n) &= 
\resizebox{0.65\linewidth}{!}{$\lim_{\lambda\to\infty} \cbr{\int_{t_0}^{t} \dif \bar{t}\, 
\sbr{\Sigma^>_\lambda(t, \bar{t}) - \Sigma^<_\lambda(t, \bar{t})} G^\rceil_\lambda(\bar{t}, i \omega_n)
 + \Sigma^\rceil_\lambda(t, i \omega_n) \mathcal{G}_\lambda(i \omega_n)}$}~,
\end{split}
\end{equation}
with an analogous equation of motion for $G^\lceil(i \omega_n, t)$.


\subsection{Analytic continuation}
\label{sec:analytic-continuation}

A series of formal analytic relations between Green functions, namely the analytic continuation~\cite{Baym_1961_2} of complex time arguments, such as the Matsubara time $\tau$ and frequencies $i\omega_n$ to real-time or frequency, can be obtained via the Lehmann representation, which is an expansion of the ensemble average over the complete set of eigenstates of the Hamiltonian. Consider the Lehmann representation of the mixed Green function~\eqref{eq:contour-gf}
\begin{equation}
\label{eq:mixed_lehmann}
\begin{split}
G_\lambda^\lceil(i \omega_n, t) =  &-i \int_0^\beta \dif\tau\, e^{+i \omega_n \tau} G^\lceil_\lambda(\tau, t)
\\
% = &-i \int_0^\beta \dif\tau e^{+i \omega_n \tau}
% \cbr{-i \tr \sbr{e^{-\beta H} e^{H \tau} d e^{-H \tau} d^\dagger\del{t}}}
% \\
= &-i \int_0^\beta \dif\tau\, e^{+i \omega_n \tau}
\sbr{- i \xi \frac{1}{Z} \sum_{n n'} \braket{n|e^{-(\hat{\mathcal{H}} + \lambda \hat Q) (\beta - \tau)} \hat d^\dagger |n'} \braket{n' | e^{-(\hat{\mathcal{H}} + \lambda \hat Q) \tau} \hat d(t) | n}}
\\
= &-i \frac{1}{Z} \sum_{n n'} \braket{n|\hat d^\dagger|n'} \braket{n'|\hat d(t)|n}
e^{-\beta E^\lambda_n} \sbr{-i \xi \int_0^\beta \dif\tau\, e^{(+i \omega_n + E^\lambda_n - E^\lambda_{n'})\tau}}
\\
% = &- \sum_{n n'} \frac{\braket{n|d|n'} \braket{n'|d^\dagger(t)|n}}{i \omega_n + E_n - E_{n'}}
% \sbr{e^{i \omega_n \beta} e^{-\beta E_{n'}} - e^{-\beta E_n}}
% \\
= &- \frac{1}{Z} \sum_{n n'} \frac{\braket{n'|\hat d^\dagger|n} \braket{n|\hat d(t)|n'}}{i \omega_n + E^\lambda_n - E^\lambda_{n'}}
\sbr{\xi e^{-\beta E^\lambda_{n'}} - e^{-\beta E^\lambda_n}}~.
\end{split}
\end{equation}
Likewise, it can be shown that $G_\lambda^\rceil(i \omega_n, t) = G_\lambda^\lceil\del{t, (i \omega_n)^*}^\dagger$ and hence only one of the two mixed Green functions needs to be considered. But more surprisingly, it also follows from the Lehmann representation that $G_\lambda^\rceil(i \omega_n, t_0) = \mathcal{G}_\lambda(i \omega_n)$. Similarly to the Matsubara Green function $\mathcal{G}_\lambda(z)$, the mixed Green functions are also analytical everywhere in the complex plane, except on the real axis, where they have a branch cut and $\mathcal{G}_\lambda(i\omega_n \to \omega \mp i \eta) \stackrel{!}{=} G_\lambda^{A/R}(\omega)$~\cite{Rammer_2007_gf}. So although appearing in form as greater/lesser Green functions, the analytically-continued mixed Green functions are instead related to the advanced/retarded Green functions. These would be equivalent to Wigner-Ville-transforming the ones presented in \cref{sec:neqft-gf}, assuming thermal equilibrium conditions.

\subsection{Dyson and Kadanoff-Baym equations}
\label{sec:project-dyson}

\subsubsection{Dyson equations}

A solution of the Green functions in \textit{real-frequency} can be obtained directly from the analytic continuation of the Matsubara Dyson equation~\eqref{eq:matsubara-gf-aux-particles}
\begin{equation}
\label{eq:dyson-real-freq-aux}
    \lim_{\lambda\to\infty} {G}_\lambda^{A/R}(\omega + \lambda) = \lim_{\lambda\to\infty} \sbr{\omega - h_0 - \Sigma^{A/R}_\lambda(\omega + \lambda)}^{-1}~,
\end{equation}
for which $G^A_\lambda(\omega)^\dagger = G^R_\lambda(\omega)$. The real-frequency greater/lesser Green functions are connected through the fluctuation-dissipation relation
\begin{equation}
\label{eq:fluctuation-dissipation-gtrless}
    \lim_{\lambda\to\infty} G_\lambda^<(\omega + \lambda) = \lim_{\lambda\to\infty} \xi\, e^{-\beta(\omega + \lambda)} G_\lambda^>(\omega + \lambda)~.
\end{equation}
However, due to the $\lambda$ limit, the exponential term usually diverges \textit{numerically} and, unlike standard thermal equilibrium techniques, both Green function components should be calculated independently. Whereas the real-frequency greater component in the $\lambda$-limit can be directly obtain through~\eqref{eq:retarded-decomp} and~\eqref{eq:dyson-real-freq-aux}, the definition of the lesser component $G_\lambda^<(\omega+\lambda)$ follows from the property $a - b = a (b^{-1} - a^{-1}) b$ and~\eqref{eq:fluctuation-dissipation-gtrless}
\begin{equation}
\label{eq:lesser-real-freq-aux}
\begin{split}
    G_\lambda^<(\omega+\lambda) &\coloneqq {G}^R_\lambda(\omega + \lambda) \sbr{\xi\, e^{-\beta (\omega+\lambda)} \del{\Sigma^R_\lambda(\omega +\lambda) - \Sigma^A_\lambda(\omega+\lambda)}} {G}^A_\lambda(\omega+\lambda)
    \\
    &={G}^R_\lambda(\omega + \lambda) \Sigma^<_\lambda(\omega + \lambda) {G}_\lambda^A(\omega + \lambda)~.
\end{split}
\end{equation}

\subsubsection{Kadanoff-Baym equations}

Similarly, the real-frequency mixed and complementary lesser Green function can be defined as
\begin{subequations}
\begin{align}
     G_\lambda^\rceil(t, \omega) &\coloneqq G_\lambda^\rceil(t, i \omega_n \to \omega + i \eta)
     \\
     G_\lambda^<(t, \omega) &\coloneqq \xi\, e^{-\beta (\omega + \lambda)} \sbr{G_\lambda^\rceil(t, \omega) - G_\lambda^\rceil(t, \omega)^\dagger}~,
\end{align}
\end{subequations}
and follow the equations of motion
\begin{subequations}
\begin{align}
\lim_{\lambda\to\infty}\sbr{i \partial_t - h(t)} G^\rceil_\lambda(t, \omega) &= 
\lim_{\substack{\zeta\to0 \\ \lambda\to\infty}} \cbr{\int_{t_0}^{t} \dif \bar{t}\, 
\Sigma^>_\zeta(t, \bar{t}) G^\rceil_\lambda(\bar{t}, \omega)
 + \Sigma^\rceil_\lambda(t, \omega) G^\rceil_\lambda(t_0, \omega)
+ \mathcal{O}(\zeta)}
\\
\begin{split}
\lim_{\lambda\to\infty}\sbr{i \partial_t - h(t)} G^<_\lambda(t, \omega) &= 
\lim_{\substack{\zeta\to0 \\ \lambda\to\infty}} \Bigg\{\int_{t_0}^{t} \dif \bar{t}\, 
\Sigma^>_\zeta(t, \bar{t}) G^<_\lambda(\bar{t}, \omega)
\\
&+ \Sigma^<_\lambda(t, \omega) G^\rceil_\lambda(t_0, \omega)^\dagger
+ \Sigma^\rceil_\lambda(t, \omega) G^<_\lambda(t_0, \omega) + \mathcal{O}(\zeta^2)\Bigg\}~.
\end{split}
\end{align}
\end{subequations}

Finally, the thermal collision integral appearing in the equation of motion~\eqref{eq:kb-projected-lesser} of the lesser Green function $G_\zeta(t,t')$ can be calculated as
\begin{equation}
\begin{split}
    &\lim_{\zeta\to0}\int_{t_0}^{t_0 - i \beta} \dif \bar{\tau}\, \Sigma^\rceil_\zeta(t, \bar{\tau}) G^\lceil_\zeta(\bar{\tau}, t') = \lim_{\lambda\to\infty}\frac{i}{\beta} \sum_{i \omega_n} \Sigma_\lambda^\rceil(t, i \omega_n) G_\lambda^\lceil(i \omega_n, t) e^{-i \omega_n 0^+}
    \\
    = &\lim_{\lambda\to\infty} \int_{-\infty}^{+\infty} \frac{\dif \omega}{2\pi} \sbr{\Sigma_\lambda^<(t, \omega + \lambda) G^\rceil(t', \omega + \lambda)^\dagger + \Sigma^\rceil_\lambda(t, \omega + \lambda) G^<(t', \omega + \lambda)}~.
\end{split}
\end{equation}

\subsection{Thermal initial conditions}

Since the $\lambda$-scaling is formulated through a \textit{thermal} density matrix, it is only valid for thermal-related observables. This is the case of thermal initial conditions~\eqref{eq:gf-initial-conditions} for the Kelysh greater/lesser components of the Konstantinov-Perel' contour (\cref{fig:konstantinov-contour}), which read
\begin{subequations}
\begin{align}
\begin{split}
\lim_{\zeta\to0}G_\zeta^>(t_0, t_0) &=  \lim_{\lambda\to\infty}\mathcal{G}_\lambda(0^+, 0)
= \lim_{\lambda\to\infty}\frac{i}{\beta} \sum_n e^{-i \omega_n 0^+} \mathcal{G}_\lambda(i \omega_n)
\\
&= \lim_{\lambda\to\infty} \int_{-\infty}^{+\infty} \frac{\dif \omega}{2\pi} \sbr{1 + \xi\, n(\omega)} \sbr{{G}^R_\lambda(\omega) - {G}^A_\lambda(\omega)}
\\
&= \lim_{\lambda\to\infty} \int_{-\infty}^{+\infty} \frac{\dif \omega}{2\pi} \sbr{{G}^R_\lambda(\omega+\lambda) - {G}^A_\lambda(\omega+\lambda)} + \lim_{\zeta\to0}\mathcal{O}(\zeta)
\end{split}
\\
\begin{split}    
\lim_{\zeta\to0}G_\zeta^<(t_0, t_0) &= \lim_{\lambda\to\infty}\mathcal{G}_\lambda(0, 0^+)
= \lim_{\lambda\to\infty}\frac{i}{\beta} \sum_n e^{+i \omega_n 0^+} \mathcal{G}_\lambda(i \omega_n)
\\
&= \lim_{\lambda\to\infty} \int_{-\infty}^{+\infty} \frac{\dif \omega}{2\pi} \xi\,n(\omega) \sbr{{G}^R_\lambda(\omega) - {G}^A_\lambda(\omega)}
\\
&= \lim_{\lambda\to\infty} \int_{-\infty}^{+\infty} \frac{\dif \omega}{2\pi} \xi\,e^{-\beta (\omega+\lambda)} \sbr{{G}^R_\lambda(\omega +\lambda) - {G}^A_\lambda(\omega+\lambda)} + \lim_{\zeta\to0}\mathcal{O}(\zeta^2)~,
\end{split}
\end{align}
\end{subequations}
where
\begin{equation}
    \lim_{\lambda\to\infty}n\del{\omega + \lambda} = \lim_{\lambda\to\infty} e^{-\beta (\omega + \lambda)} \frac{1}{1 - \xi\, e^{-\beta(\omega + \lambda)}} = \lim_{\zeta\to0} \sbr{\zeta\, e^{-\beta \omega} + \mathcal{O}(\zeta^2)}~.
\end{equation}
Note that the $\zeta$-scaling is not \textit{explicitly} encoded in the real-frequency Green functions, but on how the non-interacting energy of the auxiliary particles being at infinity effectively modifies the distribution functions. 

The initial conditions for the mixed Green functions follow directly from the boundary condition~\eqref{eq:gf-initial-conditions} at $t=t_0$
\begin{subequations}
\begin{align}
    G^\rceil_\lambda(t_0, \omega) &= G^R_\lambda(\omega)~,
    \\
    G^<_\lambda(t_0, \omega) &= G^<_\lambda(\omega)~.
\end{align}
\end{subequations}