\chapter{Non-equilibrium quantum field theory} \label{sec:neqft} The time-evolution operator $\hat{U}$ \begin{equation} \label{eq:time-evolution-op} \hat{U}(t, t') = \mathds{T}\sbr{e^{-i \int_{t'}^t \dif \bar{t}\,\hat{H}(\bar{t})}} \qquad \mathrm{for}\; t > t'~, \end{equation} where $\hat{H}(t)$ is a (time-dependent) Hamiltonian and $\mathds{T}$ the time-ordering operator, generates the time evolution of wavefunctions $\ket{\Psi(t)} = \hat{U}(t, t')\ket{\Psi(t')}$ or density matrices $\hat{\rho}(t) = \hat{U}(t, t') \hat{\rho}(t') \hat{U}(t', t)$, for all times $t$. Similarly to the wave function of a quantum system, a density matrix scales, at worst, exponentially with system size and, despite describing a quantum ensemble of particles exactly, can be impractical for the calculation of observables. This difficulty is addressed by non-equilibrium quantum field theory, which can directly reformulate the quantum problem and its time evolution in terms of the observables of interest, typically objects of much lower dimensionality. \section{The Schwinger contour} The time-ordering operator $\mathds{T}$ is not an operator in the quantum-mechanical, conventional sense -- associated with an observable -- but rather establishes a rule on how to arrange products of operators which depend on time: given a time-grid $t_1 < t_2 < \ldots < t_n$ \begin{equation} \label{eq:time-ordered-op} \mathds{T}\sbr{\hat{H}(t_{P(\ell)})\ldots\hat{H}(t_{P(2)})\hat{H}(t_{P(1)})} = \hat{H}(t_\ell)\ldots\hat{H}(t_{2})\hat{H}(t_{1})~, \end{equation} for all permutations $P$ of $\cbr{1, 2, \ldots, \ell}$. For commuting Hamiltonian operators $\sbr{\hat{H}(t_i), \hat{H}(t_j)} = 0$, $\mathds{T}$ has no influence on~\eqref{eq:time-evolution-op} and the integral can be evaluated directly. However, in general, Hamiltonian operators at different times do not commute, and $\hat{U}(t, t')$ is the continuous limit of an infinite product of locally constant Hamiltonian operators \begin{equation} \label{eq:trotterization} \hat{U}(t, t') \equiv \lim_{N\to\infty} e^{-i \hat{H}\del{t} \delta t} e^{-i \hat{H}\del{t - \delta t} \delta t} \ldots e^{-i \hat{H}\del{t - \del{N-1}\delta t} \delta t}~, \end{equation} where $\delta t = (t - t')/(N-1)$ is some infinitesimal time-step. From the group property that \begin{equation} \label{eq:time-evolution-group} \hat{U}(t, t') \hat{U}(t', t) = \hat{\mathds{1}}~, \end{equation} it can be inferred through a decomposition akin to~\eqref{eq:trotterization} that \begin{equation} \hat{U}(t, t') = \bar{\mathds{T}}\sbr{e^{+i \int_{t}^{t'} \dif \bar{t}\,\hat{H}(\bar{t})}} \qquad \mathrm{for}\, t < t'~, \end{equation} where $\bar{\mathds{T}}$ is the anti-time-ordering operator, ordering the operators in the reverse order of~\eqref{eq:time-ordered-op}. A time-dependent observable $O(t)$ is obtained by calculating the ensemble average of its associated operator $\hat{O}(t)$ in the Heisenberg picture, \begin{equation} \label{eq:ensemble-average-def} O(t) = \ev*{\hat{O}(t)} \coloneqq \frac{\tr \cbr{\hat{\rho}(t_0)\hat{O}(t)}}{\tr \cbr{\hat{\rho}(t_0)}} = \frac{\tr \cbr{\hat{\rho}(t_0)\hat{U}(t_0, t) \hat{O}(t) \hat{U}(t, t_0)}}{\tr \cbr{\hat{\rho}(t_0)}}~, \end{equation} where $\hat{\rho}(t_0)$ is the initial density matrix -- describing an arbitrary interacting or non-interacting state at $t=t_0$, with the trace taken over the Hilbert space on which the Hamiltonian acts. Note that -- despite the abuse of notation -- the operator $\hat O$ retains a time argument, which solely records the instance of time that the (constant-in-time) Schr\"odinger operator should act on. By defining an oriented time path $\gamma = (t_0, t) \oplus (t, t_0)$, the ensemble average can directly be written as a path-ordered product of operators on $\gamma$, \begin{equation} \ev*{\hat{O}(t)} = \frac{\tr \cbr{\hat{\rho}(t_0)\mathds{T}_\gamma \sbr{e^{-i \int_\gamma \dif \bar{t}\,\hat{H}(\bar{t})} \hat{O}(t)}}}{\tr \cbr{\hat{\rho}(t_0)}}~, \end{equation} where $\mathds{T}_\gamma$ is a time-ordering operator, rearranging the operators by their \textit{chronological} order on $\gamma$, i.e., preceding operators on the contour on the right. Finally, due to~\eqref{eq:time-evolution-group}, the path in the denominator can be extended to $\gamma$ yielding \begin{equation} \label{eq:ensemble-average} \ev*{\hat{O}(t)} = \frac{\tr \cbr{\hat{\rho}(t_0)\mathds{T}_\gamma \sbr{e^{-i \int_\gamma \dif \bar{t}\,\hat{H}(\bar{t})} \hat{O}(t)}}}{\tr \cbr{\hat{\rho}(t_0)\mathds{T}_\gamma \sbr{e^{-i \int_\gamma \dif \bar{t}\,\hat{H}(\bar{t})}}}}~. \end{equation} The dynamics of the system are then entirely determined by its initial state together with the Hamiltonian, and the calculation of an observable is realised by evolving the initial state forward and then backwards in time, resulting in the ubiquitous closed-time Schwinger~\cite{Schwinger_1961} contour. It is possible to eliminate the backward branch of the contour~\cite{Kamenev2011_intro} in a system in equilibrium, as the ground state in the distant future is known. However, this does not hold in non-equilibrium systems since there is no guarantee that the system can return to its ground state. \begin{figure}[!htb] \centering \includegraphics[width=\figwidth]{diagrams/sketches/schwinger-contour} \caption{The Schwinger contour. Note that the forward (denoted by subscript $(-)$) and backward (denoted by a subscript $(+)$) branches are only displaced for graphical purposes. In addition, generally $\hat{O}(t_+) $ can be different from $\hat{O}(t_-)$ for any $t \in \gamma$.} \label{fig:schwinger-contour} \end{figure} \section{Interacting initial conditions} \label{sec:contours} Unlike in a many-body interacting system at large times $t \gg t_0$, where it generally holds that it no longer contains information related to its initial state, i.e., the initial \textit{correlations} of the system have decayed with time, the transient dynamics during early times $t \gtrsim t_0$ are primarily dominated by the system's initial state at time $t_0$, encoded by the initial density matrix $\hat \rho(t_0)$. This time $t_0$ typically serves as an actual boundary between the past -- formally a preparatory stage for an interacting state -- and the onset of a qualitatively different process -- for example, a sudden change in the environment or coupling to an external field. A fundamental problem of starting from an interacting state -- containing non-Gaussian correlations -- is that it is incompatible with the Wick decomposition~\cite{Danielewicz_1984}: if the density matrix is not separable into a product of Gaussian operators, the expectation value of a product of operators does not factorise into a product of expectation values of pairs of operators. \subsection{The Keldysh contour} This problem was at first\footnote{Refer to~\cite{_pi_ka_2014} and references therein for a historical overview of the theoretical development concerning the inclusion of interacting initial conditions.} resolved by adiabatically switching on the interactions from a non-interacting state at a very distant past~\cite{Keldysh_1964}. The interacting density matrix $\hat \rho(t_0)$ is generated through the adiabatic switch \begin{equation} \hat \rho(t_0) = \hat U(t_0, -\infty) \hat \rho(t_{-\infty}) \hat U(-\infty, t_0)~, \end{equation} where $\hat \rho(t_{-\infty})$ is a non-interacting density matrix, and $\hat U$ evolves with the Hamiltonian \begin{equation} \hat H(t) = \begin{cases} \hat h_0 + e^{-\eta \envert{t - t_0}} \hat h_{\mathrm{int}}~, &t t_0$, following coherently after the (historical) state preparation. %This interpretation is equivalent to Keldysh's for a preparatory stage starting from a non-interacting state at $-\infty$. \begin{figure}[!htb] \centering \begin{subfigure}[c]{\figwidth} \includegraphics[width=\linewidth]{diagrams/sketches/keldysh-contour} \caption{Original contour~\cite{Keldysh_1964}} \label{fig:keldysh-contour} \end{subfigure} \centering \begin{subfigure}[c]{\figwidth} \includegraphics[width=\linewidth]{diagrams/sketches/velicky-contour} \caption{Time-partitioned contour~\cite{Velick__2010}} \label{fig:velicky-contour} \end{subfigure} \caption{The Keldysh contour is also often called the Schwinger-Keldysh contour. The shaded region denotes the preparatory stage, starting in some distant past.} \label{fig:keldysh-parent} \end{figure} \subsection{The Konstantinov-Perel' contour} \label{sec:kp-contour} For an interacting thermal initial state, a far more practical -- and arguably less artificial -- approach is preparing the initial state through time evolution in the complex-time plane~\cite{Danielewicz_1984, Wagner_1991}. For an interacting system in thermodynamic equilibrium with a reservoir at inverse temperature $\beta$ and chemical potential $\mu$, the density matrix is given by~\cite{Schwabl2006} \begin{equation} \label{eq:gibbs-ensemble} \hat{\rho}(t_0) = \frac{e^{-\beta (\hat{\mathcal{H}} - \mu \hat N)}}{\tr e^{-\beta (\hat{\mathcal{H}} - \mu \hat N)}} = \frac{e^{-i \int_{t_0}^{t_0 - i \beta} \dif \bar{z}\, (\hat{\mathcal{H}} - \mu \hat N)}}{\tr e^{-i \int_{t_0}^{t_0 - i \beta} \dif \bar{z}\, (\hat{\mathcal{H}} - \mu \hat N)}}~, \end{equation} where $\hat{\mathcal{H}}$ is the Hamiltonian for the system in equilibrium and $\hat N$ is the particle-number operator. Denoting the path $(t_0, t_0 - i \beta) \coloneqq \gamma^M$, the ensemble average~\eqref{eq:ensemble-average} can be written as \begin{equation} \ev*{\hat{O}(z)} = \frac{\tr \cbr{ e^{-i \int_{\gamma^M} \dif \bar{z}\, (\hat{\mathcal{H}} - \mu \hat{N})} \mathds{T}_\gamma \sbr{e^{-i \int_\gamma \dif \bar{z}\,\hat{H}(\bar{z})} \hat{O}(z)}}}{\tr \cbr{e^{-i \int_{\gamma^M} \dif \bar{z}\, (\hat{\mathcal{H}} - \mu \hat{N})}\mathds{T}_\gamma \sbr{e^{-i \int_\gamma \dif \bar{t}\,\hat{H}(\bar{t})} }}}~. \end{equation} The superscript $M$ stands for Matsubara, a formalism~\cite{Matsubara_1955} for systems in thermal equilibrium where the density matrix is treated as the (imaginary-)time-evolution operator. Defining a new oriented time path $\gamma = \gamma \oplus \gamma^M$, the cyclic property of the trace allows the integral over $\gamma^M$ to be brought inside $\mathds{T}_\gamma$, \begin{equation} \label{eq:operator-average} \ev*{\hat{O}(z)} = \frac{\tr \cbr{ \mathds{T}_\gamma \sbr{e^{-i \int_\gamma \dif \bar{z}\,\hat{H}(\bar{z})} \hat{O}(z)}}}{\tr \cbr{\mathds{T}_\gamma \sbr{e^{-i \int_\gamma \dif \bar{z}\, \hat H(\bar{z})}}}}, \end{equation} where \begin{equation} \label{eq:hamiltonian-kp-contour} \hat{H}(z) = \begin{cases} \hat{\mathcal{H}} - \mu \hat N~, &z \in \gamma^M \\ \hat{H}(z)~, &\mathrm{otherwise}~. \end{cases} \end{equation} Note that despite the ensemble average~\eqref{eq:ensemble-average} having been originally formulated for an operator with time arguments on the real-time contour, the real-time $t$ is promoted to a complex time $z$ and~\eqref{eq:operator-average} is fully general for any $z \in \gamma$. While technically not required, the path $\gamma$ can also be made independent of $z$ by inserting $\hat{\mathds{1}} \equiv \hat{U}(z, +\infty) \hat{U}(+\infty, z)$ after $\hat O(z)$ on the right-hand-side of~\eqref{eq:ensemble-average-def}, resulting in the Konstantinov-Perel' contour, shown in \cref{fig:konstantinov-contour}. \begin{figure}[!htb] \centering \includegraphics[width=\figwidth]{diagrams/sketches/konstantinov-contour} \caption{The Konstantinov-Perel' contour. There is no agreed terminology for the name of this contour, also often called the Kadanoff-Baym contour. Note that the forward $(+)$ and backwards $(-)$ branches are defined at $\im z = 0$ and only displaced for graphical purposes.} \label{fig:konstantinov-contour} \end{figure} Even though appearing trivial in its form, transforming the thermal density matrix into a part of the contour has a critical technical outcome. By incorporating the thermal averaging into the contour, perturbative expansions (via Wick's decomposition) within a thermal theory are well defined. Moreover, this general formulation embodies the Keldysh~\cite{Keldysh_1964} and Matsubara~\cite{Matsubara_1955} formalism. Due to being a positive semi-definite operator, the density matrix has a form~\cite{Stefanucci_2013} \begin{equation} \hat\rho = \frac{e^{- \hat S}}{\tr e^{-\hat S}}~, \end{equation} which suggests that further manipulation of the time contour inside the ensemble average could appropriately account for initial non-thermal interacting states. However, the treatment of such arbitrary initial conditions~\cite{Wagner_1991, Morozov_1999, Garny_2009, Berges_2004} is beyond the scope of this text. \section{Contour-ordered Green functions} \label{sec:neqft-gf} One of the central motivations behind the development of non-equilibrium quantum field theory is the calculation of $n$-point functions, also known as Green functions -- terms which will be used interchangeably. 2-point functions are fundamental objects of many-body theories, describing \textit{single-particle} excitations or statistical distributions of particles, and encode the greater part of the experimentally-accessible observables. The prototypical non-equilibrium Green function is the (2-point) contour-ordered ensemble average \begin{equation} \label{eq:contour-gf} G_{j i}(z, z') = -i \ev*{\mathds{T}_\gamma \sbr{\hat c_j(z) \hat c_i^\dagger(z')}}~, \end{equation} for $(z, z') \in \gamma$, where $\hat c^\dagger$ and $\hat c$ are creation and annihilation operators, respectively, and the indices of the operators describe some quantum number, depending on the system under consideration. Due to the algebra of the operators, \begin{equation} \mathds{T}_\gamma \sbr{\hat c_j(z) \hat c_i^\dagger(z')} = \begin{cases} \hat c_j(z) \hat c^\dagger_i(z') &\mathrm{if}\, z \succ z'~,\\ \xi \hat c^\dagger_i(z') \hat c_j(z) &\mathrm{if}\, z \prec z'~, \end{cases} \end{equation} where $\xi=1$ for bosonic ($\mathds{C}$-number algebra) and $\xi=-1$ for fermionic (Grassmann-number algebra) operators. The contour-ordered Green function can then be decomposed as \begin{equation} \label{eq:gf-time-decomposition} \begin{split} G_{j i}(z, z') &= \Theta_\gamma(z, z') \sbr{-i \ev*{\hat c_j(z) \hat c_i^\dagger(z')}} + \Theta_\gamma(z', z) \sbr{-i \xi \ev*{\hat c_i^\dagger(z') \hat c_j(z)}} \\ &= \Theta_\gamma(z, z') G_{j i}^>(z, z') + \Theta_\gamma(z', z) G_{j i}^<(z, z')~, \end{split} \end{equation} where $\Theta_\gamma(z, z')$ is a generalised Heaviside step function\footnote{The generalised Dirac delta function $\delta_\gamma(z, z')\coloneqq\od{}{z}\Theta_\gamma(z, z')$ on the contour is given by \begin{equation} \delta_\gamma(z, z') = \begin{cases} +\delta(t - t') &\textrm{if } (z \to t,z' \to t') \,\textrm{are in the forward Keldysh branch,} \\ -\delta(t - t') &\textrm{if } (z \to t,z' \to t') \,\textrm{are in the backward Keldysh branch,} \\ +i \delta(\tau - \tau') &\textrm{if } (z\to-i\tau,z'\to-i\tau') \,\textrm{are in the Matsubara branch,} \\ 0 &\textrm{otherwise}~. \end{cases} \end{equation}} on the $\gamma$ contour, \begin{equation} \Theta_\gamma(z, z') = \begin{cases} 1~, &z \succ z' \\ 0~, &\mathrm{otherwise}~. \end{cases} \end{equation} The so-called \textit{greater} ($G^>$) and \textit{lesser} ($G^<$) components owe their names to the decomposition of the contour-ordered Green function -- which equals $G^>$ when $z \succ z'$ or $G^<$ when $z \prec z'$ on the contour $\gamma$. Note, however, that the time arguments of $G^\gtrless$ do not have to fulfil such time constraints: $G^\gtrless$ are defined for \textit{all} values $(z, z') \in \gamma$. In physical terms, the characterisation of single-particle dynamics is fully encoded in $G^<$ and $G^>$, which describe amplitudes related to the propagation of a hole or particle in a fermionic many-body system, respectively. Note that the (anti-)commutation of particles holds for all times $G^>(z,z) - G^<(z,z) = -i$. \subsection{Components of contour-ordered Green functions} \label{sec:neqft-components} Depending on the contour branches the operators act on, the contour-ordered Green function can be decomposed into different \textit{components}. For $(z,z')$ lying on the Konstantinov-Perel' contour (\cref{fig:konstantinov-contour}), the components of $G(z,z')$ read \begin{equation} \begin{split} G(z, z') &= % \resizebox{0.5\linewidth}{!}{$ \begin{pmatrix} G(t_-, t'_-) &G(t_-, t'_+) &G(t_-, t_0 - i\tau')\\ G(t_+, t'_-) &G(t_+, t'_+) &G(t_+, t_0 - i \tau')\\ G(t_0 - i \tau, t'_-) &G(t_0 - i \tau, t'_+) &G(t_0 - i \tau, t_0 - i \tau') \end{pmatrix} % $} \\ &= \resizebox{0.85\linewidth}{!}{$ % \setlength\arraycolsep{-0pt} \begin{pmatrix} \Theta(t - t')G^>(t,t') + \Theta(t' - t)G^<(t,t') &G^>(t, t') &G^\rceil(t, \tau') \\ G^<(t, t') &\Theta(t - t')G^<(t,t') + \Theta(t' - t)G^>(t,t') &G^\rceil(t, \tau') \\ G^\lceil(\tau, t) &G^\lceil(\tau, t) &\mathcal{G}(\tau, \tau') \end{pmatrix} $}~, \end{split} \end{equation} with $(t_0 - i \tau, t_0 - i \tau')$ lying on the imaginary time branch and $(t_-,t_-')$ and $(t_+, t'_+)$ on horizontal time forward and backwards branches, respectively, and the relation\footnote{For most systems of interest, the Hamiltonian is the same whether the time argument is the upper or lower branch, $\hat U(t_-, t'_-) = \hat U(t_+, t'_+)$, and the relation $\hat O(t_-) = \hat O(t_+)$ similarly holds for operators in the Heisenberg picture.} $\hat O(t_-) = \hat O(t_+)$ being employed in the last equality. As a result, the distinction between horizontal branches in the time arguments of the Green functions is dropped since $G^\gtrless(t_-, z') = G^\gtrless(t_+, z')$ and $G^\gtrless(z, t'_-) = G^\gtrless(z, t'_+)$. The two time branches are simply a formal outcome of unitary time evolution -- forward and backward evolution brings the state back to its original state -- and can be dropped once operator precedence~\eqref{eq:contour-gf} is considered. This is also expected physically since there is just a single time branch in the description of reality, as can be observed in the definition of the time-evolution operator~\eqref{eq:time-evolution-op}. The contour-ordered Green function on the Konstantiv-Perel' contour then contains 5 independent components, which can be divided into three groups: $\mathcal{G}(\tau, \tau')$, $\cbr{G^<(t, t'), G^>(t, t')}$ and $\cbr{G^\rceil(t, \tau'), G^\lceil(\tau, t')}$. First is the \textit{Matsubara} Green function, which contains information about the system in its initial thermal equilibrium stage. It is hence decoupled from the other four components, which describe dynamics in real-time, and determines their initial conditions. The second group is the Keldysh Green functions, which carry information about single-particle correlations in real-time. The third group are the mixed Green functions, which carry information about vestigial single-particle thermal correlations in the system. These usually decay with time, except in specific integrable systems or meta-stable states~\cite{Wagner_1991}. Note that the components of $G(z,z')$ for $(z,z')$ on the Keldysh contour (\cref{fig:keldysh-parent}) are the same as the Keldysh components of $G(z,z')$ on the Konstantinov-Perel' contour (\cref{fig:konstantinov-contour}). \subsection{Interpretability of Keldysh Green functions} \label{sec:neqft-interpretability} \subsubsection{The Wigner basis} \label{sec:wigner-basis} There is a lack of interpretability of 2-point functions in real time. Whereas in equilibrium, by definition, they display time-translation symmetry and can be described in a frequency basis that directly relates to the system's energy spectrum, it is not immediately intuitive how to interpret a non-equilibrium Keldysh Green function. The Wigner basis \begin{equation} \label{eq:wigner-rotation} G(T, \tau)_\mathrm{W} = G(T+\tau/2, T-\tau/2)~, \end{equation} is obtained by rotating and squeezing the real-time coordinates $(t, t')$ \begin{equation} \begin{pmatrix} \sqrt{2} & 0 \\ 0 & \frac{1}{\sqrt{2}} \\ \end{pmatrix} \cdot \begin{pmatrix} \cos\frac{\pi}{4} & -\sin\frac{\pi}{4} \\ \sin\frac{\pi}{4} & \cos\frac{\pi}{4} \\ \end{pmatrix} \cdot \begin{pmatrix} t \\ t' \\ \end{pmatrix} = \begin{pmatrix} 1 & -1 \\ \frac{1}{2} & \frac{1}{2} \\ \end{pmatrix} \cdot \begin{pmatrix} t \\ t' \\ \end{pmatrix} \eqqcolon \begin{pmatrix} \tau \\ T \\ \end{pmatrix} ~, \end{equation} for which $T$ is called the \textit{centre-of-mass} time and $\tau$ the \textit{relative} time, and their derivatives are given by \begin{equation} \partial_T = \partial_t + \partial_{t'}, \quad \partial_\tau = \frac{\partial_t - \partial_{t'}}{2}~. \end{equation} Since equilibrium physics is characterised by time-translation invariance -- i.e., independence in $T$ -- equilibrium 2-point functions are just dependant on $\tau$, and are commonly expressed in the dual Fourier basis $\omega$. The Wigner-Ville transform \begin{equation} \label{eq:wigner-ville} G(T, \omega)_{\tilde{\mathrm{W}}} = \int_{-\infty}^{+\infty} \dif\tau\, e^{i \omega \tau} G(T, \tau)_\mathrm{W} = \int_{-\infty}^{+\infty} \dif\tau\, e^{i \omega \tau} \sbr{\int_{-\infty}^{+\infty} \frac{\dif\omega'}{2\pi} e^{-i \omega' \tau}G(T, \omega')_{\tilde{\mathrm{W}}}}~, \end{equation} provides a convenient linear transformation, where it is possible to establish generalisations of equilibrium properties or interpret non-equilibrium 2-point functions from the viewpoint of equilibrium 2-point functions. Note, however, that Wigner-Ville-transformed functions are \textit{non-causal} since at each time-slice $T$ the frequencies $\omega$ contain information from times \textit{before} and \textit{after} $T$. As such, when far from stationarity, interpretation of Wigner-Ville-transformed 2-point functions must be taken with a grain, or rock, of salt. \subsubsection{Other contour representations} Even though the greater/lesser components are the natural decompositions of the contour-ordered Green function~\eqref{eq:gf-time-decomposition}, other representations of non-equilibrium Green functions exist. A possible representation is given by \begin{subequations} \begin{align} \rho_{ji}(z,z') &= -i \ev*{\sbr{c_j(z), c_i^\dagger(z')}_{-\xi}} \equiv G_{ji}^>(z,z') - G_{ji}^<(z,z') \\ F_{ji}(z,z') &= \frac{1}{2} \ev*{\sbr{c_j(z), c_i^\dagger(z')}_\xi} \equiv \frac{i}{2}\del{G_{ji}^>(z,z') + G^<_{ji}(z,z')}~, \end{align} \end{subequations} where $\rho$ and $F$ are coined as the \textit{spectral} and \textit{statistical} functions, respectively. Unlike $G^\gtrless$, which describe the propagation of a hole/particle excitation, the spectral function roughly encodes the probability density of the available physical states (excitations), and the statistical function encodes the occupation of these states. This interpretation is particularly appropriate for fermionic particles in equilibrium, where Wigner-Ville-transformed $\rho(\cdot, \omega)_{\tilde{\mathrm{W}}}$ is positive definite everywhere with its integral over all $\omega$ equalling $1$, and can be thought of as the probability density of an excitation having some energy $\omega$. Note the alternative definition of the contour-ordered Green function \begin{equation} \label{eq:gf-statistical-decomp} G(z,z') \equiv \frac{1}{2}\sgn_\gamma(z-z')\rho(z,z') - i F(z,z')~. \end{equation} Another standard representation, however conceptually very similar, for non-equilibrium Green functions in the real-time coordinates $(t, t')$ can be obtained via the Keldysh rotation~\cite{Kamenev2011_bosons, Kamenev2011_fermions} \begin{subequations} \begin{align} G^\mathrm{R}(t, t') &= \Theta(t - t') \sbr{G^>(t,t') - G^<(t, t')} \\ G^\mathrm{A}(t, t') &= \Theta(t' - t) \sbr{G^<(t,t') - G^>(t, t')} \\ G^\mathrm{K}(t, t') &= G^>(t, t') + G^<(t,t')~, \end{align} \end{subequations} where the superscripts $\mathrm{R}$, $\mathrm{A}$ and $\mathrm{K}$ stand for \textit{retarded}, \textit{advanced} and \textit{Keldysh}, respectively. The retarded and advanced functions carry spectral information, and the Keldysh function statistical information about the system's single-particle excitations. There is also some historical importance due to the strong connection of the retarded and advanced with the Matsubara Green function in equilibrium (\cref{sec:analytic-continuation}). \subsubsection{The fluctuation-dissipation relation} \label{sec:fdr} After its excitation and subsequent relaxation stage\footnote{Refer to~\cite{Berges_2004, Aoki2014, Bonitz_2016} for discussions on the characteristic timescales of non-equilibrium dynamics.}, a non-equilibrium system will evolve towards a stationary state -- which is understood as displaying some time-translational invariance. However, this does not necessarily imply that the system has \textit{thermalised} and can be described by a Gibbs distribution~\eqref{eq:gibbs-ensemble}, with counterexamples found in, e.g., non-equilibrium steady-states, metastable or Floquet states. Notably, by construction~\eqref{eq:time-evolution-op}, non-equilibrium time evolution is unitary and cannot lose information about its initial state, which is at odds with the fundamental property of thermal systems being memoryless. Furthermore, the concept of temperature (and associated statistical distributions) is also completely absent from the formalism -- despite temperature possibly being encoded in the initial density matrix, it is no longer a well-defined property of the system for $t > t_0$. Nonetheless, it has often been observed that a system \textit{can} reach a state displaying identical features to thermal systems, for which an estimate of the effective temperature of the system -- or generalised statistical distribution can be determined, assuming that a fluctuation-dissipation relation holds. The fluctuation-dissipation relation establishes a deep relation between a system's statistical (fluctuation) and spectral (dissipation) information in thermal equilibrium. This is encoded in the Kubo-Martin-Schwinger (KMS) relations, for which the (anti-)periodicity $\mathcal{G}(\tau, \tau) = \xi \mathcal{G}(\tau + \beta, \tau')$ of Matsubara Green functions holds\footnote{Note that the argument $\tau$ in Matsubara Green functions denotes the imaginary time $i \tau$ in the Konstantinov-Perel' contour (\cref{fig:konstantinov-contour}) and not the relative time $\tau$ of the Wigner basis.}. Satisfying the KMS conditions for~\eqref{eq:gf-statistical-decomp} in a Wigner-Ville-transformed~\eqref{eq:wigner-ville} basis yields \begin{equation} -\frac{1}{2}\rho(\cdot, \omega)_{\tilde{\mathrm{W}}} - i F(\cdot, \omega)_{\tilde{\mathrm{W}}} = \xi e^{-\beta \omega} \sbr{\frac{1}{2}\rho(\cdot, \omega)_{\tilde{\mathrm{W}}} - i F(\cdot, \omega)_{\tilde{\mathrm{W}}}}~, \end{equation} for which the \textit{fluctuation-dissipation relation} reads \begin{equation} F(\cdot, \omega)_{\tilde{\mathrm{W}}} = i \sbr{\frac{1}{2} + \xi n(\omega)} \rho(\cdot, \omega)_{\tilde{\mathrm{W}}}~, \end{equation} where $n(\omega)$ is the Bose-Einstein distribution for bosonic or Fermi-Dirac for fermionic 2-point functions. Note that this delicate balance -- or interdependence -- between occupations and spectra results from the KMS boundary conditions, which holds only for systems in thermal equilibrium. Despite generally the absence of analogous boundary conditions in non-equilibrium, the relation can be observed for systems that have thermalised, motivating a generalisation of the fluctuation-dissipation relation \begin{equation} \label{eq:fdr} F(T, \omega)_{\tilde{\mathrm{W}}} = i \sbr{\frac{1}{2} + \xi n(T, \omega)_{\tilde{\mathrm{W}}}} \rho(T, \omega)_{\tilde{\mathrm{W}}}~, \end{equation} where $n(T, \omega)_{\tilde{\mathrm{W}}}$ describes a time-dependant distribution function and $\rho(T, \omega)_{\tilde{\mathrm{W}}}$ is a generalisation of the equilibrium spectral function that can roughly describe how the spectral density changes with the centre-of-mass time $T$. \subsubsection{Separation of timescales} \label{sec:separation-timescales} A common, albeit non-general, feature in non-equilibrium 2-point functions is the notion of \textit{fast} and \textit{slow} variables\footnote{Consider a non-interacting 2-point function of a particle with energy $\omega_0$ and average particle number $\bar{n}$ \begin{equation} G(T, \tau)_{\mathrm{W}} = -i \sbr{\Theta_\gamma(\tau) + \xi \bar{n}} e^{i \omega_0 \tau}~. \end{equation} This function has a fast dependence in $\tau$ and (infinitely) slow dependence in $T$. Adding a small disturbance $\omega_0 \to \omega_0 + \Delta\omega(T)$ will result in some dependence in $T$, however \textit{slower} than in $\tau$~\cite{Maciejko2007}.}. A separation of timescales~\cite{Picano_2021} is found when the scale $\Lambda$ of the of time-evolution -- with $\Lambda\to\infty$ for a system in a stationary state -- is much larger than the inverse of the width of the spectral features $\Omega$ \begin{equation} \envert{\frac{\partial_T G(T, \omega)_{\tilde{\mathrm{W}}}}{G(T, \omega)_{\tilde{\mathrm{W}}}}} < \Lambda^{-1}~, \qquad \envert{\frac{\partial_\omega G(T, \omega)_{\tilde{\mathrm{W}}}}{G(T, \omega)_{\tilde{\mathrm{W}}}}} < \Omega^{-1}~. \end{equation} For example, by expressing time-integrals as a Moyal product (with the exponentials emerging from the Wigner rotation followed by expressing the translations via their generator, i.e., $f(T+s) = e^{s\, \partial_T}f(T)$) \begin{equation} \int \dif \bar{t} \,A(t,\bar{t}) B(\bar{t},t') = e^{\frac{i}{2}\sbr{\partial^A_T \partial^B_\omega - \partial^B_T \partial^A_\omega}} A(T, \omega)_{\tilde{\mathrm{W}}} B(T, \omega)_{\tilde{\mathrm{W}}} \stackrel{\Lambda \gg \Omega^{-1}}{\approx}A(T, \omega)_{\tilde{\mathrm{W}}}B(T, \omega)_{\tilde{\mathrm{W}}}~, \end{equation} convolutions between 2-point functions are reduced to local frequency products in centre-of-mass time. Despite breaking causality (similarly to~\eqref{eq:fdr}), for very slow transients, the non-equilibrium problem is reduced to a quasi-equilibrium problem, which can greatly reduce the complexity of resolving non-equilibrium dynamics. \subsection{The Langreth rules} \label{sec:langreths-rules} The prescriptions for decoding simple convolutions or products of contour-ordered Green functions into their constituting components are known as the Langreth rules. The principal ingredient for their derivation is the decomposition~\eqref{eq:gf-time-decomposition}, with generalisations of these \enquote{contour calculus} rules presented in~\cite{Hyrk_s_2019}. For a simple contour-ordered product \begin{equation} C(z, z') = A(z, z') B(z', z) = \Theta_\gamma(z, z') A^>(z,z')B^<(z',z) + \Theta_\gamma(z', z) A^<(z, z')B^>(z',z)~, \end{equation} the components' products read \begin{equation} C^\gtrless(z, z') = A^\gtrless(z, z') B^\lessgtr(z', z)~. \end{equation} However, the resulting $C^\gtrless$ is not a proper greater/lesser Green function as it does not fulfil the symmetry relation (\cref{sec:symmetries}) associated with these functions. The components of the simplest contour-ordered convolution \begin{equation} \label{eq:countour-conv} \begin{split} C(z, z') &= \int_\gamma \dif \bar{z} \,A(z, \bar{z}) B(\bar{z}, z') \\ &= \Theta_\gamma(z,z')\int_{z'}^z \dif \bar{z}\,A^>(z,\bar{z})B^>(\bar{z},z') + \Theta_\gamma(z',z)\int_{z}^{z'} \dif \bar{z}\,A^<(z,\bar{z})B^<(\bar{z},z') \\ &\phantom{=} + \int_{t_0}^{\min_\gamma(z,z')} \dif \bar{z}\,A^>(z,\bar{z})B^<(\bar{z},z') +\int_{\max_\gamma(z,z')}^{t_0-i\beta} \dif \bar{z}\,A^<(z,\bar{z})B^>(\bar{z},z')~, \end{split} \end{equation} with $\min_\gamma(z,z')$ and $\max_\gamma(z,z')$ are $\min$ and $\max$ functions generalised to the contour, with the arguments being compared by their \textit{precedence} in the contour, read \begin{subequations} \begin{align} \begin{split} C^\rceil(t, \tau') &= \int_{t_0}^{t} \dif \bar{t}\, \sbr{A^>(t, \bar{t}) - A^<(t, \bar{t})}B^\rceil(\bar{t}, \tau') % \\ % &\phantom{=} - i \int_{0}^{\beta} \dif \bar{\tau}\, A^\rceil(t, \bar{\tau}) B(\bar{\tau}, \tau')~, \end{split} \\ \begin{split} C^\lceil(\tau, t') &= -\int_{t_0}^{t'} \dif \bar{t}\, A^\lceil(\tau, \bar{t}) \sbr{B^>(\bar{t}, t') - B^<(\bar{t}, t')} % \\ % &\phantom{=} - i \int_{0}^{\beta} \dif \bar{\tau}\, A(\tau, \bar{\tau}) B^\lceil(\bar{\tau}, t')~, \end{split} \\ \begin{split} C^\gtrless(t, t') &= \int_{t_0}^t \dif\bar{t}\, \sbr{A^>(t, \bar{t}) - A^<(t, \bar{t})}B^\gtrless(\bar{t}, t') - \int_{t_0}^{t'} \dif\bar{t}\, A^\gtrless(t, \bar{t}) \sbr{B^>(\bar{t}, t') - B^<(\bar{t}, t')} \\ &\phantom{=} - i \int_{0}^{\beta}\dif\tau A^\rceil(t, \tau) B^\lceil(\tau, t')~. \end{split} \end{align} \end{subequations} \subsection{Calculation of contour-ordered Green-functions} \subsubsection{Initial conditions} \label{sec:gf-initial-conditions} For problems on the Konstantinov-Perel' contour, the initial conditions of the Keldysh and mixed Green functions are implicitly determined by the Matsubara Green function $\mathcal{G}(\tau, \tau')$ \begin{equation} \label{eq:gf-initial-conditions} \begin{split} G^<(t_0, t_0) = \mathcal{G}(0, 0^+)~, &\qquad\qquad G^>(t_0, t_0) = \mathcal{G}(0^+, 0)~, \\ G^\rceil(t_0, \tau') = \mathcal{G}(0, \tau')~, &\qquad\qquad G^\lceil(\tau, t_0) = \mathcal{G}(\tau, 0)~. \end{split} \end{equation} Problems with contours such as \cref{fig:keldysh-contour} have their initial conditions explicitly determined by $\hat \rho(t_0)$. For problems with contours such as \cref{fig:velicky-contour}, their initial conditions are Green functions defined for all times smaller than $t_0$. For example, for thermal initial conditions, the preparatory stage of the system can be obtained via an inverse Wigner-Ville transform (cf.~\eqref{eq:wigner-ville}) of the equilibrium $2$-point function $G^\gtrless(\omega)$, calculated in real frequency (cf.~\cref{sec:project-dyson}) \begin{equation} \label{eq:wigner-ville-inv} G^\gtrless(T < t_0, \tau)_\mathrm{W} = \int \frac{\dif\omega}{2\pi} e^{-i \omega \tau} G^\gtrless(\omega)~, \end{equation} followed by an inverse Wigner rotation (cf.~\eqref{eq:wigner-rotation}). Encoding the system's initial condition in such a manner precludes the use of mixed Green functions and can greatly simplify the complexity of the problem. \begin{figure}[!htb] \centering \includegraphics[width=\figwidth]{diagrams/sketches/wigner-rot} \caption{An inverse Wigner-Ville transformation followed by an inverse Wigner rotation of an equilibrium spectral function $\rho(\omega)$.} \label{fig:wigner-rot} \end{figure} \subsubsection{Equations of motion} Contour-ordered Green functions are nothing more than ensemble averages of operators in the Heisenberg picture, which obey the equations of motion \begin{equation} \label{eq:martin-schwinger} \od{}{t} \ev*{\hat O(t)}= \od{}{t} \ev*{\hat{U}(t_0, t) \hat{O}(t) \hat{U}(t, t_0)} = i \ev*{\sbr{\hat H, \hat O(t)}} + \ldots \end{equation} Nonetheless, calculating Green functions for an interacting $\hat H$ in this manner results in an infinite hierarchy of differential equations, with an ever-increasing order of operator averages. This hierarchy is known as Martin-Schwinger's, and while theoretically describing any many-body system exactly, it constitutes an intractable problem. Truncating the hierarchy via bare perturbative schemes in non-equilibrium settings can violate conservation laws due to spurious \textit{secular} terms, which grow with time~\cite{Berges_2004}. It also may not preserve non-linear features of the theory, which are necessary for coherent effects~\cite{Cornwall_1974}. This non-linearity -- e.g., self-consistency in the equations describing the system's dynamics -- is ultimately related to the emergence of \textit{universality}, where for $t \gg t_0$ the dynamics should be insensitive to the initial conditions, and the system may thermalise. A more sophisticated approach is required to address these issues appropriately. \section{Non-equilibrium two-particle-irreducible effective action} \label{sec:2pi} The $n$-particle-irreducible ($n$PI) effective actions are a class of field theories that provide practical and systematic approximations for classical~\cite{Bode_2022} and quantum physics~\cite{Berges_2004}. Their main advantage is that, unlike in bare perturbative approaches where all connected diagrams (graphs) contribute to the perturbative expansion, in $n$PI approaches, only $n$PI diagrams do -- these are graphs that do not become disconnected once $n$ lines are cut. Despite still being of perturbative nature, $n$PI approaches resum an infinite number of diagrams belonging to a particular class. The 2PI effective action is the simplest action~\cite{Carrington_2004} that both generates the non-linearity needed for the \textit{universality} requirement of non-equilibrium theories and eliminates the spurious \textit{secular} terms present in bare perturbative schemes. These ensure a \textit{conserving} -- also known as $\Phi$-derivable~\cite{Baym_1961, Baym_1962} -- approximation where a series of conservation laws, such as conservation of particle number are fulfilled, despite truncations of the perturbative series. \subsection{Path-integral construction} \label{sec:path-integral} The starting point is the construction of a functional path-integral representation of the partition function (cf.~\eqref{eq:ensemble-average}) \begin{equation} \label{eq:partition-function} Z = \tr \cbr{\hat\rho(t_0) \mathds{T}_\gamma \sbr{e^{-i \int_\gamma \dif \bar{t}\,\hat{H}(\bar{t})}}}~. \end{equation} For the sake of brevity, $\gamma$ is taken as the Keldysh contour (\cref{fig:keldysh-contour}) -- other contours, such as Konstantinov-Perel's (\cref{fig:konstantinov-contour}) follows a similar derivation, however with more contour degrees of freedom (\cref{sec:neqft-components}). Considering a \textit{normal-ordered}\footnote{A normal-ordered product of operators has all creation operators to the left of the annihilation operators.} Hamiltonian $H(t) = H\sbr{\hat b^\dagger(t), \hat b(t),t}$, the path integral is formulated via a Trotter-type limiting decomposition~\eqref{eq:trotterization} of~\eqref{eq:partition-function}. A complete set \begin{equation} \mathds{1} = \int \dif\sbr{\phi^*,\phi}\, e^{-\envert{\phi}^2} \ket{\phi}\bra{\phi}~, \end{equation} of appropriately time-labelled coherent states\footnote{A coherent state $\ket{\phi} \coloneqq e^{\xi \phi \hat b^\dagger}\ket{0}$ is the eigenstate of the annihilation operator $\hat b \ket{\phi} = \phi \ket{\phi}$, where $\ket{0}$ denotes the vacuum, with $\xi=1$ for a complex scalar and $\xi=-1$ for a complex Grassmann number $\phi$.} are inserted between the time-evolution operators of the Trotterized form of the partition function, with $\dif\sbr{\phi^*,\phi} = \frac{\dif \phi^*\,\dif \phi}{2\pi i}$ for complex scalar and $\dif\sbr{\phi^*,\phi} = \dif \phi^*\,\dif \phi$ for Grassmann fields. Expressing the trace as \begin{equation} \tr \hat O = \int \dif\sbr{\phi^*,\phi}\, e^{-\envert{\phi}^2} \braket{\phi | \hat O | \phi}~, \end{equation} this results in the discrete path-integral formulation \begin{equation} \begin{split} Z &= \int \lim_{N\to\infty} \prod_{j=0}^{N-1} \dif\sbr{\phi_{j_+}^*,\phi_{j_+}} \dif\sbr{\phi_{j_-}^*,\phi_{j_-}} \braket{\phi_{1_-} | \hat \rho(t_{-\infty}) | \phi_{1_+}} e^{i S\sbr{\phi}}~, \end{split} \end{equation} with, \begin{equation} \label{eq:discrete-action} \begin{split} \resizebox{0.95\linewidth}{!}{$S[\phi] = \sum_{j=1}^{N-1} \delta t_j \sbr{\del{+i \phi^*_{j_-}\, \frac{\phi_{j_-} - \phi_{j_- -1}}{\delta t_j} - H\sbr{\phi^*_{j_-}, \phi_{j_- -1}}} % \\ - \del{-i \phi^*_{j_+}\, \frac{\phi_{j_+} - \phi_{j_+ -1}}{\delta t_j} - H\sbr{\phi^*_{j_+}, \phi_{j_+ -1}}} } + i \phi^*_{(N-1)_+} \phi_{(N-1)_-}$}~, \end{split} \end{equation} where $\delta t_j = t_j - t_{j-1}$ and $j_\mp$ denotes $j$-th Trotter slice, and the $-$ and $+$ subscripts denote the time forward and backward branches of the Keldysh contour, respectively. In the limit ${N \to \infty}$, $\phi$ is promoted to a \textit{field} $\phi(z)$. The appearance of the last term in the discrete action $S$ arises at the inflexion of the contour. Here, the field at $\phi(t_{N_-})$ and $\phi(t_{N_+})$ is indistinguishable, and hence there is no time evolution operator between the associated coherent states. The initial distribution of the system is encoded by $\hat \rho(t_{0})$, for which can be shown~\cite{Kamenev2011_bosons, Kamenev2011_fermions} that $\braket{\phi_{0_-} | \hat \rho(t_{0}) | \phi_{0_+}} = \exp\del{\xi \phi^*_{0_-}\phi_{0_+}\rho}$. Even though the Keldysh action in continuum form \begin{equation} \label{eq:keldysh-action} \begin{split} S[\phi] % &= \int_{-\infty}^{+\infty} \dif t\, \cbr{\phi^*(t_-) \sbr{i \partial_{t_-} - h(t_-)} \phi(t_-) - \phi^*(t_+) \sbr{i \partial_{t_+} - h(t_+)} \phi(t_+)} % \\ &= \int_{\gamma} \dif z\, \phi^*(z) G^{-1}(z,z) \phi(z)~, % \\ % &\equiv % \int_{-\infty}^{+\infty} \dif t\, % \begin{pmatrix} % \phi^*(t_-) & \phi^*(t_+) % \end{pmatrix} % \begin{pmatrix} % i \partial_{t_-} - h(t_-) &0 % \\ % 0 &i \partial_{t_+} - h(t_+) % \end{pmatrix} % \begin{pmatrix} % \phi(t_-) % \\ % \phi(t_+) % \end{pmatrix} \end{split} \end{equation} where \begin{equation} G^{-1}(z,z') = \delta_\gamma(z, z') \sbr{i \partial_z - H(z)}~, \end{equation} appears to have both contour branches decoupled, that they are connected through the boundary terms of the discrete action. On a similar vein, $\delta_\gamma(z, z')$ is not a proper Dirac $\delta$-distribution. The \textit{causal} construction of the path-integral connects the fields between adjacent Trotterized slices of the path-integrals~\eqref{eq:discrete-action} and is technically defined as $\delta_\gamma(z,z'+0^+)$. \subsection{Two-particle-irreducible effective action construction} Omitting the explicit dependence on the initial average and sourcing $\phi$, considering it to be a multi-component field, the cumulant generating function $W$ is defined as \begin{equation} \label{eq:path-integral-Z} \begin{split} Z\sbr{\bm{j},\bm{K}} &= \int D\bm{\phi} \exp\cbr{i \sbr{S\sbr{\bm{\phi}} + \del{ \bm{\phi}^*_z \bm{j}_z + \bm{j}^*_z \bm{\phi}_z} + \bm{\phi}^*_z \bm{K}_{z z'} \bm{\phi}_{z'}}} \equiv e^{i\,W\sbr{\bm{j}, \bm{K}}}~, \end{split} \end{equation} The $n$-point functions can be generated through functional derivatives with respect to the 1- and 2-point external source fields $\bm{j}$ and $\bm{K}$, obeying the same algebra as $\bm{\phi}$. Note the DeWitt notation where repeated, and continuous indices are integrated and summed over \begin{equation} \bm{j}_z\bm{\phi}_z = \sum_i \int_\gamma \dif z\, j_i(z) \phi_i(z)~. \end{equation} The cumulant-generating functional $W\sbr{\bm{j},\bm{K}}$ generates the $1$- and $2$-point functions via \begin{subequations} \begin{align} \xi \frac{\delta}{i \delta j_i(z)} i W\sbr{\bm{j},\bm{K}} &= \overline{\phi^*}_i(z)~, \\ \frac{\delta}{i \delta j^*_i(z)} i W\sbr{\bm{j},\bm{K}} &= {\overline\phi}_i(z)~, \\ \xi \frac{\delta}{i \delta K_{ij}(z,z')} i W\sbr{\bm{j},\bm{K}} &\equiv i G_{ij}(z,z')~, \end{align} \end{subequations} where the $\xi$ factor arises due to the anti-commutative algebra of complex Grassmann fields (cf.~\eqref{eq:contour-gf}). Note that $G_{ij}(z,z')$ is not a second cumulant but a second \textit{moment} of the quantum distribution. Even though these distinctions are often awkwardly overlooked\footnote{This is rooted in the fact that in the absence of mean-fields -- i.e., vanishing 1-point functions -- the second cumulant and moment are identical. This equality always holds for 2-point functions of fermionic or bosonic fields with no condensate part. However, it does not generally hold for higher-order-point functions.} in many-body physics, there is an important distinction when characterising \textit{correlations}, which becomes especially relevant for higher-order-point functions. Cumulants are the objects that describe \textit{de facto} $n$-point interactions as they cannot be factorised into products of lower-order interactions. This has the formal consequence that when expressing an $n$-point interaction as a sum of graphs (Feynman diagrams), only \textit{connected} graphs are considered. These cumulants are termed connected\footnote{In the jargon of field-theory, this terminology is intimately related to the linked-cluster theorem, where perturbative corrections to the $n$th cumulant are determined by connected graphs with $n$ lines. Note that vacuum bubble diagrams are factored out through derivatives of $W\sbr{\bm{j}, \bm{K}}$. Notwithstanding, $n$-point functions can still have disconnected contributions, i.e., products of lower-order point diagrams, which ought to be removed, resulting in \textit{connected} $n$-point functions.} $n$-point functions, $G^\mathrm{(connected)}(z_1, \ldots, z_n)$, and are calculated via \begin{equation} \xi^n \frac{\delta}{i \delta j_{i_1}(z)}\ldots\frac{\delta}{i \delta j^*_{i_{2n}}(z_{2n})} i W\sbr{\bm{j},\bm{K}}= i G^\mathrm{(connected)}_{i_1, \ldots, i_{2n}}(z_1, \ldots, z_{2n})~, \end{equation} for which the alternative expression \begin{equation} G^\mathrm{(connected)}_{i j}(z,z') = -i \xi \sbr{\frac{\delta}{i \delta K_{ij}(z,z')} i W\sbr{\bm{j},\bm{K}} - \overline\phi_i(z) \overline{\phi^*}_j(z')} \end{equation} is deduced. The superscript $\mathrm{(connected)}$ will be dropped, and all $n$-point functions from this point onwards are considered connected. In analogy with transformations in thermodynamics, a new generating functional describing the same physics can be formulated~\cite{Cornwall_1974} depending on the conjugate variables $\bm{\overline\phi}$ and $\bm{G}$. The so-called two-particle-irreducible (2PI) effective action $\Gamma\sbr{\bm{\overline\phi}, \bm{G}}$ is obtained through the triple Legendre transform \begin{equation} \begin{split} \Gamma\sbr{\bm{\overline\phi}, \bm{G}} &= \resizebox{0.85\linewidth}{!}{$W\sbr{\bm{j}, \bm{K}} - \int_\gamma \dif z\, \fd{W\sbr{\bm{j},\bm{K}}}{j_i(z)} j_i(z) - \int_\gamma \dif z\, \fd{W\sbr{\bm{j},\bm{K}}}{j^*_i(z)} j^*_i(z) - \int_\gamma \dif z \dif z'\, \fd{W\sbr{\bm{j},\bm{K}}}{K_{ij}(z,z')}K_{ji}(z',z)$} \\ &= W\sbr{\bm{j}, \bm{K}} - \bm{\overline{\phi^*}}_z \bm{j}_z - \bm{j^*}_z \bm{\overline\phi}_z - \bm{\overline{\phi^*}}_{z} \bm{K}_{z z'} \bm{\overline\phi}_{z'} - i \bm{G}_{z z'} \bm{K}_{z' z}~. \end{split} \end{equation} Apart from the advantageous functional dependence on the observables of interest, the 2PI effective action obeys a variational principle for vanishing sources $\bm{j}$ and $\bm{K}$. In this limit, the sourceless physical system is recovered, and the equation of motion for the $2$-point functions is encoded as a stationary condition of the effective action, \begin{subequations} \begin{align} \fd{\Gamma\sbr{\bm{\overline\phi}, \bm{G}}}{\overline\phi_k(z)} &= - \xi j^*_k(z) - \int_\gamma \dif{z'}\,\overline{\phi^*}_i(z') K_{ik}(z',z) \\ \fd{\Gamma\sbr{\bm{\overline\phi}, \bm{G}}}{\overline{\phi^*}_k(z)} &= - j_k(z) - \int_\gamma \dif{z'}\, K_{kj}(z,z') \overline\phi_j(z') \\ \label{eq:2pi-fd} \frac{\delta\Gamma\sbr{\bm{\overline\phi}, \bm{G}}}{\delta G_{kk'}(z,z')} &= -i K_{k'k}(z',z)~. \end{align} \end{subequations} The \textit{effective action}, consisting of a classical action plus quantum corrections, is obtained by expanding~\cite{Calzetta_2008} the original quantum fields $\bm{\phi}$ as a classical plus a fluctuation term $\bm{\phi} = \bm{\overline{\phi}} + \bm{\varphi}$, \begin{equation} \begin{split} e^{i\Gamma\sbr{\bm{\overline\phi}, \bm{G}}} &= \exp\cbr{i W\sbr{\bm{j}, \bm{K}} - i \del{\bm{\overline{\phi^*}}_z \bm{j}_z + \bm{j^*}_z \bm{\overline\phi}_z + \bm{\overline{\phi^*}}_{z} \bm{K}_{zz'} \bm{\overline\phi}_{z'} + i \bm{G}_{zz'} \bm{K}_{z'z}}} \\ &=\int D\bm{\varphi} \exp\resizebox{0.75\linewidth}{!}{$\cbr{i \sbr{S\sbr{\bm{\overline{\phi}}+\bm{\varphi}} + \bm{\varphi}^*_z (\bm{j}_z + \bm{K}_{zz'} \bm{\overline\phi}_{z'}) + (\bm{j}^*_z + \bm{\overline{\phi^*}}_{z'} \bm{K}_{z'z}) \bm{\varphi}_z + \bm{\varphi}_z^* \bm{K}_{zz'} \bm{\varphi}_{z'}}} e^{\bm{K}_{zz'} \bm{G}_{z'z}}$}~. \end{split} \end{equation} At their core, $n$PI methods, much like most perturbative theoretical techniques, are based on splitting the action into a non-interacting part (bilinear in the fields) and an interacting part, \begin{equation} S\sbr{\bm{\phi}} = S_0\sbr{\bm{\phi}} + S_\mathrm{int}\sbr{\bm{\phi}} = \bm{\phi}^*_z \bm{G}_{0_{zz'}}^{-1} \bm{\phi}_{z'} + S_\mathrm{int}\sbr{\bm{\phi}}~, \end{equation} where \begin{equation} \label{eq:inv-gf} \bm{G}_{0_{zz'}}^{-1} = \delta_\gamma(z,z') \sbr{i \partial_z \mathds{1} - \bm{h}_0(z)}~, \end{equation} and $\bm{h}_0(z)$ is the non-interacting part of the Hamiltonian $\bm{H}(z)$. For the sake of brevity, consider vanishing 1-point functions $\bm{\overline\phi} = \bm{\overline{\phi^*}} = 0$. The one-loop order 2PI effective action is obtained by discarding the interacting part and evaluating the complex Gaussian integral~\cite{Berges_2004, Stoof_2009, fraboulet:tel-03466730} \begin{equation} \begin{split} \Gamma^{\mathrm{(1 loop)}}\sbr{\bm{G}} &= -i \xi \tr \log \sbr{i \del{\bm{G}^{-1}_0 - \bm{K}}} - i \xi \tr \bm{K} \bm{G} % \\ % &= -i \xi \tr \log i \bm{G}^{-1} - i \xi \tr \sbr{\del{\bm{G}^{-1}_0 - \bm{G}^{-1}} \bm{G}} % \\ % & = i \xi \tr \log i \bm{G} - i \xi \tr \sbr{\bm{G}^{-1}_0 \bm{G}} + \const~, \end{split} \end{equation} where~\eqref{eq:2pi-fd} was used to remove the explicit dependency on $\bm{K}$. Going beyond one loop, the 2PI effective action reads \begin{equation} \Gamma\sbr{\bm{G}} = i \xi \tr \log i \bm{G} - i \xi \tr \sbr{\bm{G}^{-1}_0 \bm{G}} + \Gamma_2\sbr{\bm{G}} + \const~. \end{equation} The functional $\Gamma_2\sbr{\bm{G}}$ contains all contributions beyond one-loop, i.e., the sum of all scattering effects due to interactions (\cref{sec:neqft-2pi-loop}). For a theory with both bosonic (complex scalar) and fermionic (complex Grassmann) fields, the 2PI effective action can be written as \begin{equation} \label{eq:gamma2-functional} \Gamma\sbr{\bm{G}, \bm{D}} = +i \tr \ln i\bm{G} - i \tr \bm{G}_0^{-1}\bm{G} - i \tr \ln i\bm{D} + i \tr \bm{D}_0^{-1}\bm{D} + \Gamma_2\sbr{\bm{G}, \bm{D}} + \const~, \end{equation} where $\bm{G}$ denotes 2-point bosonic functions and $\bm{D}$ the fermionic 2-point functions, assuming that 1-point bosonic functions vanish. \subsection{Contour Dyson and Kadanoff-Baym equations} \label{sec:dyson-kbe} Deriving the equations of motion for the 1-point and 2-point functions through a variational procedure ensures that the global symmetries and conservation laws of the effective action are preserved. For vanishing source fields, these are given by the stationarity conditions of the 2PI effective action: \begin{equation} % \fd{\Gamma\sbr{\bm{\overline\phi}, \bm{G}, \bm{D}}}{\bm{\overline\phi}} = 0~, % \quad \fd{\Gamma\sbr{\bm{G}, \bm{D}}}{\bm{G}} = 0~, \quad \fd{\Gamma\sbr{\bm{G}, \bm{D}}}{\bm{D}} = 0~. \end{equation} Notably, one finds the Dyson equations for the 2-point functions \begin{subequations} \begin{align} \label{eq:inv-dyson} G^{-1}_{ij}(z,z') &= G^{-1}_{0, ij}(z,z';\bm{\overline\phi}) - \Pi_{ij}(z,z'; \bm{\overline\phi}, \bm{G}, \bm{D}) \\ D^{-1}_{ij}(z,z') &= D^{-1}_{0, ij}(z,z';\bm{\overline\phi}) - \Sigma_{ij}(z,z'; \bm{\overline\phi}, \bm{G}, \bm{D}) ~, \end{align} \end{subequations} where the self-energies $\bm{\Pi}$ and $\bm{\Sigma}$ are defined as \begin{subequations} \begin{align} \label{eq:2pi-self-energy} \Pi_{ij}(z,z') &\equiv +i \fd{\Gamma_2\sbr{\bm{G}, \bm{D}}}{G_{ji}(z',z)} \\ \Sigma_{ij}(z,z') &\equiv -i \fd{\Gamma_2\sbr{\bm{G}, \bm{D}}}{G_{ji}(z',z)}~. \end{align} \end{subequations} The Dyson equations are of little use for time-dependent non-equilibrium problems due to requiring the inversion of dense operators. While appropriate in equilibrium problems, where owing to time-translation invariance (\cref{sec:neqft-interpretability}), the equations become diagonal in the Fourier basis and can easily be inverted, this is far from the case in non-equilibrium where there is no symmetry in time. However, noting that \begin{equation} \bm{G}^{-1}_{z \bar{z}} \bm{G}_{\bar{z} z'} = \mathds{1}_{zz'} \Leftrightarrow \int_{\bar{z}} G^{-1}_{ik}(z, \bar{z}) G_{kj}(\bar{z}, z') = \delta_\gamma(z, z')\delta_{ij}~, \end{equation} the equations can be transformed into an initial value problem: \begin{subequations} \begin{align} \bm{G}_{0_{z\bar{z}}}^{-1}\bm{G}_{\bar{z}z'} &= \mathds{1}_{zz'} + \bm{\Pi}_{z\bar{z}} \bm{G}_{\bar{z}z'}~, \\ \bm{G}_{z\bar{z}}\bm{G}_{0_{\bar{z}z'}}^{-1} &= \mathds{1}_{zz'} + \bm{G}_{z\bar{z}}\bm{\Pi}_{\bar{z}z'}~. \end{align} \end{subequations} Given the form of the inverse non-interacting 2-point functions~\eqref{eq:inv-gf}, the resulting equations of motion are two-time integrodifferential equations, also known as the Kadanoff-Baym equations: \begin{subequations} \label{eq:kb} \begin{align} \sbr{i \vec{\partial}_z \mathds{1}- \bm{h}_0(z)} \bm{G}(z, z') &= \delta_\gamma(z,z') \mathds{1} + \int_\gamma \dif\bar{z}\, \bm{\Pi}(z,\bar{z})\bm{G}(\bar{z}, z') \\ \bm{G}(z, z') \sbr{-i \cev{\partial}_{z'} \mathds{1} - \bm{h}_0(z')} &= \delta_\gamma(z,z')\mathds{1} + \int_\gamma \dif\bar{z}\, \bm{G}(z, \bar{z})\bm{\Pi}(\bar{z}, z') ~. \end{align} \end{subequations} The equations of motion for $\bm{D}(z,z')$ follow a similar derivation and structure. A distinct feature of Kadanoff-Baym equations is their non-Markovian structure, evident from the integrals on the right-hand side of the integrodifferential equations~\eqref{eq:kb}. The appearance of such \textit{memory} effects is a formal consequence that dynamics of (2-point) correlations can be fully described with knowledge from all other (2-point) correlations. This is a natural trade-off from the reduction of the state space from, e.g., the differential equations~\eqref{eq:martin-schwinger} generating the Martin-Schwinger hierarchy -- with Markovian structure but dependence on \textit{all} $n$-point functions. Similar to the equations generated by the hierarchy, the Kadanoff-Baym equations are formally exact. However, in practical terms, the exact calculation of $\Gamma_2$, which generates the self-energies, is an intractable problem. \subsection{Two-particle-irreducible loop expansion} \label{sec:neqft-2pi-loop} Being a Legendre transform of the cumulant generating functional $W$, the 2PI effective action is an exact method and, in principle, contains the full information about the system. However, due to the great complexity of evaluating the $\Gamma_2$ functional exactly, approximations to the 2PI effective action arise from the truncation of $\Gamma_2$. By definition, $\Gamma_2$ contains all closed, topologically distinct, 2PI diagrams constructed from the bare vertices of the theory. This can succinctly~\cite{Rammer_2007_ea} be expressed through \begin{equation} \label{eq:gamma2} \begin{split} \Gamma_2\sbr{\bm{G}, \bm{D}} &= -i \ev*{\sum_{n=1}^\infty \frac{(i S_{\mathrm{int}})^n}{n!}}_{\bm{G}, \bm{D} \,\mathrm{\& 2PI}}~, \end{split} \end{equation} where $\ev*{}_{\bm{G}, \bm{D} \,\mathrm{\& 2PI}}$ denotes that Wick's decomposition\footnote{The Wick theorem or decomposition is neatly understood in the path-integral formalism: for a \textit{Gaussian} distribution -- i.e., a non-interacting theory, any high-order moment can be decomposed in products of second moments (or cumulants). That is, $\ev*{X_1 X_2 \ldots X_{n-1} X_{n}} = \sum_{p_n} \prod_{\cbr{i,j} \in p} \ev*{X_i X_j}$ for even $n$ and zero for odd $n$, where $p_n$ denotes all possible ways of arranging $\cbr{1,\ldots,n}$ in pairs $\cbr{i,j}$.} is used to express the products of field operators into sums of products of pairs of field operators ($2$-point functions), that these $2$-point functions are the connected Green functions $\bm{G}$ and $\bm{D}$, and that only 2PI vacuum bubble (closed) diagrams are considered. The simplest $\Gamma_2$ truncation scheme is the loop expansion. This essentially mimics a coupling expansion, i.e., a perturbation expansion in powers of the interaction vertex, where a small coupling parameter $|V_0| \lesssim 1$ is required for convergence of the series. Since the $\Gamma_2$ diagrams are composed of the full, interacting $\bm{G}$ and $\bm{D}$, the self-consistency arising from the back-coupling to the self-energies~\eqref{eq:kb} sums arbitrarily high powers of the coupling term. For illustrative purposes, consider $S_{\mathrm{int}}$, containing the vertices of some theory with 3 quantum fields -- $c$, $f$ and $b$ -- in diagrammatical form: \begin{equation} i S_{\mathrm{int}} = \int_\gamma \dif \bar{z} \sbr{V_0 c^*(\bar{z}) b^*(\bar{z}) f(\bar{z}) + \hc} = \input{diagrams/feynman-diagrams/phi3-vertex} \end{equation} Given a specific power of $i S_\mathrm{int}$, the vertices are connected in all possible ways, and all non-2PI diagrams are discarded, resulting in \begin{equation} \label{eq:gamma2-diagrams} \Gamma_2 = -i \sbr{\input{diagrams/feynman-diagrams/gamma2_nophoton} + \frac{1}{3}\input{diagrams/feynman-diagrams/gamma2_higher-order} + \ldots} \end{equation} It can also be argued that due to the small coupling, higher-order diagrams will have a smaller contribution to the physics of the problem; hence, only the lowest order terms of $\Gamma_2$ have to be taken into account. Self-energy diagrams must be 1PI, which is fulfilled as the functional derivative required to calculate the self-energies~\eqref{eq:2pi-self-energy} amounts to cutting one line of the $\Gamma_2$ function graphs, which are 2PI by construction.