\chapter{Numerical solution of Kadanoff-Baym equations}
\label{sec:vide}

\paragraph{Contribution statement} This \namecref{sec:vide} is based on the section \enquote{\nameref*{sec:vide}} of the following article
\fcite[true]{Meirinhos_2022}
\noindent which is presented here with minor editorial changes and with the formatting is adapted to match the thesis style. Such reproduction is permitted as the journal articles published by SciPost are licensed under the \href{https://creativecommons.org/licenses/by/4.0/}{\textit{Creative Commons Attribution 4.0 International (CC BY 4.0)}} license, which allows adaptation, reproduction and redistribution of the original source for any purpose, as long as the licence terms are followed. In~\cite{Meirinhos_2022}, the concept of adaptivity in the numerical solution of ordinary differential equations is extended to Kadanoff-Baym equations, which are two-time (integro)differential equations that arise as equations of motion of non-equilibrium 2-point functions. Ensuring the time-step size is large can dramatically increase the numerically-accessible integration times since the complexity of the numerical solution of Kadanoff-Baym equations scales at least cubically and quadratically in the number of time steps, regarding numerical operations and memory requirements, respectively. Following rather well-established literature on the numerical solution of non-stiff ordinary differential equations, adaptivity is implemented in the control of the time-step size and the order of the interpolation polynomials. The former ensures the time-step size is optimal -- i.e., as large as possible -- given the tolerance criteria, and the latter ensures the polynomial order that incurs the least amount of local error is chosen. This is achieved by mapping a discretised two-time to a one-time stepping scheme and  standard adaptivity algorithms and heuristics. It is shown how an adaptive scheme in the numerical solution of 2-point functions in time has a significant edge over fixed step-size/order schemes at a minimal computational and implementation overhead. First, it is observed that an adaptive time-stepping scheme can successfully integrate an exactly solvable model in considerably fewer time steps than a fixed scheme while incurring a smaller global error. Second, it is shown that an adaptive time-stepper effectively changes the time-step size in time regions where, for example, a stronger drive is acting on the system, and the integrator should step slower to keep the local error small. In this example, an overall factor of $5$ between the smallest and largest time step critically accelerates the numerical integration. This adaptiveness could be particularly beneficial in studying transients such as non-equilibrium quenches, where there is a violent change in the system at early times, which induces rapid changes and hence requires small time steps. Afterwards, this is typically followed by a pre-thermalisation regime, where the system changes slowly while evolving towards an equilibrium state. The evolution to an equilibrium state should permit larger time steps since the dynamics are almost static in the centre-of-mass direction. Note, however, that the adaptivity is not formulated in Wigner coordinates (\cref{sec:wigner-basis}) and hence slow time evolution in the centre-of-mass direction cannot be fully leveraged. The versatility of the adaptive scheme is tested in several physical models, ranging from quantum to classical stochastic systems. Moreover, within the quantum systems, the adaptive algorithm is tested in fermionic tight-binding and Fermi-Hubbard models as well as in bosonic mixture and open boson dimer models. Within the publication~\cite{Meirinhos_2022}, I was responsible for developing the adaptive scheme's idea, theory and numerical implementation. Furthermore, I wrote several parts of the text and took part in developing and analysing the examples provided. Specifically, within the section \enquote{\nameref*{sec:vide}}, transcribed in this \namecref{sec:vide}, I was responsible for the entirety of its contents.

The computation of solutions to the Kadanoff-Baym equations consists formally in finding numerical solutions to an integrodifferential equation of the form
\begin{equation}
\label{eq:volterra_general}
    i \partial_t g(t,t') = h_0(t) g(t,t') + \int_{\gamma'} \dif\bar{t}\, K(t,\bar{t}) g(\bar{t},t')~,
\end{equation}
which, together with its adjoint, spans the entire $(t, t')$ plane. Note that $g(t, t')$ and the kernel $K(t, t')$ are assumed to be either skew-Hermitian or symmetric with respect to their arguments. While at first glance~\eqref{eq:volterra_general} may look like a Fredholm integral equation~\cite{Wazwaz_2011}, in physical systems the integrals $\int_{\gamma'} \dif\bar{t}$ are always reduced to Volterra form, i.e.\ $\gamma'=[t_0, t]$  or $\gamma'=[t_0, t']$ (cf.~\eqref{eq:countour-conv}), the deeper reason for this being causality. Since $K$ is usually a functional of $g$, \eqref{eq:volterra_general} belongs to the class of generic non-linear Volterra integrodifferential equations (VIDE). For the rest of the analysis,  assume that the integral kernel is smooth and non-singular, as the converse is rarely encountered in the class of physical problems considered here and would require problem-dependent modifications of the quadrature rules to be properly accounted for~\cite{Brunner_1982}.

The fact that the VIDE~\eqref{eq:volterra_general} is defined on a two-dimensional domain has not only obfuscated its analysis it has also impeded a direct application of most existing numerical algorithms, which have primarily been focused on univariate VIDEs. An appropriate discretisation scheme is presented, which allows the application of general linear methods for solving the Kadanoff-Baym equations and an exposition of the variable Adams method, the preferred multi-step method for solving these equations.

\section{Stepping Scheme for Kadanoff-Baym Equations}
\label{sec:stepping-scheme}

Due to the causal structure of the Volterra initial-value problem, the Kadanoff-Baym equation at the point $(t,t')$ is only dependent on time arguments smaller or equal to $(t,t')$. By taking the Cartesian product of a (non-equidistant) one-dimensional grid 
\begin{equation}
    \mathcal{T} \coloneqq \cbr{t_0 < t_1 < \ldots < t_i < \ldots < t_N}
\end{equation}
with itself, a symmetric mesh $\mathcal{T}\times\mathcal{T} = \cbr{(t, t') \mid t\in\mathcal{T}, t' \in\mathcal{T}}$ for the two-time domain is obtained. Within such a discretisation, the time-stepping procedure can be regarded as a fan-like stepping in the symmetric two-time mesh, as depicted in \cref{fig:time-stepping}. Accordingly, this can be understood as a system of \textit{univariate}, vector-valued differential equations
\begin{equation}
\label{eq:vert_horiz_diag}
\begin{split}
    \phantom{-}i \partial_{t_i} \mathbf{g}^v(t_i) &= h_0(t_i) \mathbf{g}^v(t_i) + \del{\mathbf{K} \circ \mathbf{g}}^v(t_i) \quad \ \, \textrm{(vertical step)}
    \\
    -i \partial_{t^{}_i} \mathbf{g}^h(t^{}_i) &= \mathbf{g}^h(t^{}_i) h_0(t^{}_i)^\dagger + \del{\mathbf{g} \circ \mathbf{K}}^h(t^{}_i) \quad \textrm{(horizontal step)}
    \\
    \phantom{-}i \partial_{t_i} \mathbf{g}^d(t_i) &= h_0(t_i) \mathbf{g}^d(t_i) - \mathbf{g}^d(t_i) h_0(t_i)^\dagger + \del{\mathbf{K} \circ \mathbf{g}- \mathbf{g} \circ \mathbf{K}}^d(t_i)  \quad \textrm{(diagonal step)}~,
\end{split}
\end{equation}
where 
\begin{equation}
\begin{split}
    \mathbf{g}^v(t_i) &= \sbr{g(t_i, t_0), g(t_i, t_1), \ldots, g(t_i, t_{i})}~,
    \\
    \mathbf{g}^h(t_i) &= \sbr{g(t_0, t_i), g(t_1, t_i), \ldots, g(t_{i}, t_i)}~,
    \\
    \mathbf{g}^d(t_i) &= \sbr{g(t_i, t_i)}~,
\end{split}
\end{equation}
and $\circ$ denotes the element-wise Volterra integration
\begin{equation}
    (\mathbf{A} \circ \mathbf{B})^v(t_i) = \sbr{\int_{\gamma'}{\dif\bar{t}\,} A(t_i, \bar{t}) B(\bar{t}, t_0), \int_{\gamma'}{\dif\bar{t}\,} A(t_i, \bar{t}) B(\bar{t}, t_1), \ldots, \int_{\gamma'}{}{\dif\bar{t}\,} A(t_i, \bar{t}) B(\bar{t}, t_i)}~,
\end{equation}
with analogous definitions for the $h$ and $d$ components.

Kadanoff-Baym equations are set apart from univariate ordinary differential equations (ODEs) or VIDEs by the fact that their dimension grows with each time-step — the size of $\mathbf{g}^v(t)$ and $\mathbf{g}^d(t)$ grows by one when stepping from $t_{i}$ to  $t_{i+1}$. This requires a continued resizing of the equations and is one reason why such equations are not straightforwardly compatible with the extensive amount of available ODE solvers. Moreover, unlike population-growth problems, for example, where the size of the equations may also grow with time, the new equations that are added when solving Kadanoff-Baym equations have a \textit{past}. This can be visualised via \cref{fig:time-stepping} by noting that, for example, when stepping vertically or horizontally from $g(t_4, t_4)$, the right-hand side of the differential equations for the new elements in $\mathbf{g}^v(t)\vert_{t=t_4}$ and $\mathbf{g}^h(t)\vert_{t=t_4}$ involve in general non-zero terms at times $t < t_4$. For multi-step methods, in particular, this may necessitate additional care (cf.~\cref{sec:adams}).

\begin{figure}[t]
  \centering
  \includegraphics[]{diagrams/sketches/time-stepping}
  \caption{Time-stepping procedure for~\eqref{eq:kb}.}
  \label{fig:time-stepping}
\end{figure}

Viewing the Kadanoff-Baym integration procedure effectively as a one-time ODE problem has two main benefits: First, it opens up the possibility of applying virtually any general linear method to solve Kadanoff-Baym equations. And second, additional one-time functions such as mean fields (first cumulants) can be solved simultaneously and in a unified manner, allowing direct method implementations with well-defined local error estimations.

\section{Univariate Volterra Integrodifferential Equations}

Following the structure presented in~\eqref{eq:vert_horiz_diag}, consider a univariate non-linear VIDE in standard form, i.e.
\begin{equation}
    y'(t) = F[t, y(t)] + \int_{t_0}^{t}{\dif s\,} K[t, s, y(s)]~,
\end{equation}
which can also be seen as a system of two equations, of which one is an ordinary differential equation and the other a Volterra integral equation,
\begin{equation}
\label{eq:vide-analytic}
\begin{split}
    y'(t) &= F[t, y(t)] + z(t)~, \\
    z(t) &= \int_{t_0}^{t}{\dif s\,} K[t, s, y(s)]~,
\end{split}
\end{equation}
subject to the initial condition
\begin{equation}
    y(t_0) = y_0~.
\end{equation}
In some cases, it is possible to solve such equations with analytic methods~\cite{Wazwaz_2011}, yet this usually requires the integral kernel to have specific properties such as linearity $K\sbr{t, s, y(s)}=K(t,s)y(s)$, which is not the case for most physical systems of interest. Hence, one must resort to discrete methods.

While there are many methods one can employ to solve ODEs, \textit{a priori}, there is no best method. Its choice strongly depends on factors such as stiffness, desired accuracy and function evaluation cost. A variable order and variable step size Adams (predictor-corrector) method provides a good trade-off between cost (two function evaluations per step) and overall accuracy, even when the number of equations is very large, as is indeed the case with Kadanoff-Baym equations, where the number of equations roughly equals the dimension of $G(t_0, t_0)$ times the number of time-steps.

In methods based on integration, \eqref{eq:vide-analytic} is integrated from $t_n$ to $t_{n+1}$
\begin{equation}
\label{eq:vide-discrete}
\begin{split}
    y(t_{n+1}) &= y(t_n) + \int_{t_n}^{t_{n+1}}{\dif s\,}  \cbr{F[s, y(s)] + z(s)}~,
\end{split}
\end{equation}
and the integrals are then evaluated with interpolating quadrature formulas. Here it becomes clear that the main computational bottleneck in solving these equations is in the computation of $z(t)$, which can be evaluated with a so-called direct quadrature method
\begin{equation}
\label{eq:vide-quadrature}
    z(t_n) = \int_{t_0}^{t}{\dif s\,} K\sbr{t_n, s, y(s)} = \sum_{\ell=0}^{n-1} \int_{t_\ell}^{t_{\ell+1}}{\dif s\,}  K\sbr{t_n, s, y(s)}~.
    \end{equation}
Nonetheless, it is possible to differentiate $z(t)$ further and treat $\cbr{y'(t), z'(t)}$ as a system of coupled differential equations~\cite{Houwen_1983}, which would be more suitable in cases where the integral equation is stiff ($-{\partial K}/{\partial y} \gg 1 $)~\cite{van_der_Houwen_1981}. It was opted for the former due to its more straightforward implementation and because most physical systems of interest do not satisfy such stiffness criterion. 

\section{Variable Adams method}
\label{sec:adams}

The variable Adams method~\cite{Hairer_1993} is a predictor-corrector scheme where the integrand of~\eqref{eq:vide-discrete} is approximated by a Newton polynomial, that is, an interpolation polynomial for previously computed points. A prediction $y^*_{n+1}$ for the solution of $y(t_{n+1})$ (note that here $^*$ denotes the prediction, not complex conjugation) is obtained via an explicit method with a $(k-1)$-th order polynomial
\begin{equation}
\label{eq:adams-explicit}
y^*_{n+1} = y_n + \int_{t_n}^{t_{n+1}} \dif s\, \sum_{j=0}^{k-1} \sbr{\prod_{i=0}^{j-1}(s - t_{n-i})} \delta^j \cbr{F\sbr{t_n, y(t_n)} + z(t_n)}~,
\end{equation}
and the divided differences are defined recursively as
\begin{equation}
\label{eq:recursive-derivative}
\begin{split}
    \delta^0 F\sbr{t_\ell, y(t_\ell)} &= F\sbr{t_\ell, y(t_\ell)}~,
    \\
    \delta^{j} F\sbr{t_\ell, y(t_\ell)} &= \frac{\delta^{j-1} F\sbr{t_\ell, y(t_\ell)}- \delta^{j-1} F\sbr{t_{\ell-1}, y(t_{\ell-1})}}{t_\ell - t_{\ell-j}}~.
\end{split}
\end{equation}
The prediction for $y(t_{n+1})$ is now corrected via an implicit method, where the $k$-th order interpolation polynomial of the integrand depends on the predicted value $y^*_{n+1}$:
\begin{equation}
\label{eq:adams-implicit}
    y_{n+1} = y^*_{n+1} + \int_{t_n}^{t_{n+1}} \dif s\, \sbr{\prod_{i=0}^{k-1}(s - t_{n-i})} \delta^k \cbr{F[t_{n+1}, y(t_{n+1})] + z(t_{n+1})}~.
\end{equation}

The integrals in~\eqref{eq:vide-quadrature} can be evaluated in the same predictor-corrector manner:
\begin{equation}
\begin{split}
    z^*_n &= 
    \sum_{\ell=0}^{n-1}\int_{t_\ell}^{t_{\ell+1}} \dif s\, \sum_{j=0}^{k-1} \sbr{\prod_{i=0}^{j-1}(s - t_{\ell-i})} \delta^j K_n\sbr{t_\ell, y(t_\ell)}~,    \\
    z_n &= z^*_n + \sum_{\ell=0}^{n-1}\int_{t_\ell}^{t_{\ell+1}} \dif s\, \sbr{\prod_{i=0}^{k-1}(s - t_{\ell-i})} \delta^k K_n\sbr{t_{\ell+1}, y(t_{\ell+1})}~,
\end{split}
\end{equation}
with divided differences defined as
\begin{equation}
\begin{split}
    \delta^0 K_n\sbr{t_\ell, y(t_\ell)} &= K\sbr{t_n, t_\ell, y(t_\ell)}~,
    \\
    \delta^{j} K_n\sbr{t_\ell, y(t_\ell)} &= \frac{\delta^{j-1} K_n\sbr{t_\ell, y(t_\ell)}- \delta^{j-1} K_n\sbr{t_{\ell-1}, y(t_{\ell-1})}}{t_\ell - t_{\ell-j}}~.
\end{split}
\end{equation}
The main difficulties when evaluating the predictor-corrector equations~\eqref{eq:adams-explicit} and~\eqref{eq:adams-implicit} are that it is challenging to obtain a closed formula for the integrals and that it is algorithmically expensive to calculate the divided differences via recursive formulas~\eqref{eq:recursive-derivative}. While for equidistant time grids, the equations find a simple and compact form~\cite{Hairer_1993}, in the non-equidistant case, the expressions rapidly become convoluted and complicated to implement. These problems can be circumvented by recurrence formulas~\cite{Hairer_1993}, making the evaluation of the integrals and $j$-th derivatives more efficient. 

In between time steps, an estimate of the local truncation error can be obtained by computing $\tilde{y}_{n+1} - y_{n+1}$, where $\tilde{y}_{n+1}$ is the result of the implicit step using a $(k+1)$-th order formula. It is assumed that as $k\to\infty$, the error approaches zero (in which case the integral quadrature formula is said to be \textit{convergent}). A measure of this error satisfying specific tolerances is obtained via 
\begin{equation}
\label{eq:error-def}
    le_{k}(n+1) \coloneqq \frac{\tilde{y}_{{n+1}} - y_{{n+1}}}{\texttt{atol} + \texttt{rtol} \cdot \max\del{\envert{y_{n}}, \envert{y_{{n+1}}}}}~,
\end{equation}
for which the integration step is accepted if
\begin{equation}
    \enVert{le_{k}(n+1)}_p \le 1~,
\end{equation}
and the norm is defined as
\begin{equation}
\label{eq:ode_norm}
    \enVert{x}_p = \del{\frac{1}{n}\sum_i^n \envert{x^i}^p}^\frac{1}{p}~,
\end{equation}
where typically $p=2$. Given this acceptance criterion, the roles of the tolerances $\texttt{rtol}$ and $\texttt{atol}$ in~\eqref{eq:error-def} can be better understood considering them separately under the infinity-norm. In this scenario, $-\log_{10}\texttt{rtol}$ controls the minimum number of correct digits between time steps, while $\texttt{atol}$ is a threshold for the magnitude of the elements of $y$ for which the minimum number of correct digits is guaranteed.
This local error is then used to adjust both the step size $h_{n} \coloneqq (t_{n+1} - t_{n})$ and the order $k$. The next time step is chosen as the largest possible step that still satisfies the local error being $\lesssim 1$. Given the current local error $\enVert{le_{k}(n+1)} \simeq C h_{n}^{k+1}$ for some constant $C$, and assuming that the subsequent error is maximal, i.e.\ $ \enVert{le_{k}(n+2)} \simeq C h_{n+1}^{k+1} \approx 1$ , the next time-step can be chosen optimally as~\cite{Hairer_1993}
\begin{equation}
    h_{n+1}= h_n \, \enVert{le_{k}(n+1)}_p^{-\frac{1}{k+1}}~.
\end{equation}
Obtaining the optimal order $k$ is slightly more involved, and the reader is referred to~\cite{Hairer_1993} for an excellent and self-contained explanation of heuristic mechanisms for order selection. However, regardless of the order $k$, the number of required function evaluations per time step is constant. Hence $k$ is ideally set to a large value ($\gtrsim 5$) such that the integrator can take larger steps and the overall computational cost is reduced.

\section{Volterra Integral Equations of the Second Kind}
\label{subsec:vie2}

More elaborate self-energy approximations ($GW$, $T$-matrix \cite{vonFriesen_2010}, $1/\mathcal{N}$ \cite{Berges_2004}), which comprise resummations of particular classes of diagrams, require the solution of Volterra integral equations of the second kind~\cite{Jones_1985}:
\begin{equation}
    \label{eq:vie2}
    I(t, t') = \Phi(t, t') - \int_{\mathcal{C}}\dif\bar{t}\, \Phi(t, \Bar{t}) I(\Bar{t}, t')~.
\end{equation}
In the mentioned self-energy approximations, the kernel $K(t,t')$ of~\eqref{eq:volterra_general} then typically depends linearly on $I(t,t')$ and $\Phi(t,t')$ is a function of $g(t,t')$.

There are several ways of solving~\eqref{eq:vie2}: by inversion of the triangular system of equations obtained when discretising in the same manner as in~\eqref{eq:vert_horiz_diag}, by reduction to a VIDE through differentiation, or by iteration of the equation~\cite{Schluenzen_2019}. Since~\eqref{eq:vie2} has to be solved simultaneously with~\eqref{eq:volterra_general}, reducing it to a VIDE would be ideal, yet this generally results in stiff equations~\cite{Houwen_1983} for which the variable Adams method (\cref{sec:adams}) is not appropriate. Therefore, to achieve congruity with the method previously presented, \eqref{eq:vie2} is solved iteratively at every predictor and corrector step, i.e.\ following the same evolution procedure as depicted in \cref{fig:time-stepping}.

\section{Leveraging Symmetries and Physical Properties}

Leveraging symmetries and other physical properties of a system can significantly reduce the computational effort on top of what can be achieved by adaptive time-stepping.

\subsection{Symmetries in the Two-Time Domain}
\label{sec:symmetries}

Apart from the symmetries of the Hamiltonian, the two-time Green functions encountered in quantum and classical systems possess symmetries in the two-time domain $(t, t')$. For example, in the quantum case, the numerical implementations are based on the greater and lesser Green functions that are skew-Hermitian in time,
\begin{equation}
    \label{eq:skew_hermitian_sym}
    \sbr{\bm{G}^\lessgtr(t, t')}^\dagger = -\bm{G}^\lessgtr(t', t)~.
\end{equation}
Hence, the solutions are entirely determined by either the upper- or lower-triangular elements, which essentially cuts in half the number of equations by requiring only the integration of $\mathbf{G}^d$ and either $\mathbf{G}^v$ or $\mathbf{G}^h$. However, similar relations hold for classical systems with different symmetry relations~\cite{Meirinhos_2022}.

\subsection{Memory truncation}
\label{sec:memory-truncation}

The clustering decomposition principle~\cite{Weinberg_1995_cluster} ensures that at a large-enough time separation of the physical operators, any $n$-point function factorises. In terms of \textit{connected} 2-point functions, this has the signature of exponential or power-law \textit{decay} in the relative-time direction (\cref{sec:wigner-basis}) for massive and massless fields, respectively~\cite{Banks_2008}. This principle should hold for any stable, long-lived state, an example being thermalised systems \cite{Berges_2001} described by a Gibbs ensemble. This effect is likewise present in physical systems connected to some kind of reservoir (e.g. as in open quantum systems \cite{lappe2021non}, or quantum impurity systems described by dynamical mean-field theory \cite{schuler2018truncating}). The VIDE can hence often be approximated by a Volterra \textit{delay}-integrodifferential equation with 
\begin{equation}
    z(t) = \int_{t-\tau_\mathrm{cutoff}}^{t}{\dif s\,} K\sbr{t, s, y(s)}~,
\end{equation}
where $\tau_\mathrm{cutoff}$ is some cut-off time. This is rooted in the fact that the physical Green functions in such systems display long-time decay and, thus
\begin{equation}
    \enVert{\int_{t-\tau_\mathrm{cutoff}}^{t}{\dif s\,} K\sbr{t, s, y(s)}}
    \gg
    \enVert{\int_{t_0}^{t - \tau_\mathrm{cutoff}}{\dif s\,} K\sbr{t, s, y(s)}}~.
\end{equation}
Since one bottleneck when solving Kadanoff-Baym equations is in the evaluation of the integrals, introducing a cutoff time can dramatically reduce the computational complexity from $\mathcal{O}\del{n^2 k(n)}$ to $\mathcal{O}\del{n^2 k(N_{\tau_\mathrm{cutoff}})}$ where $n$ denotes the number of time-steps, $N_{\tau_\mathrm{cutoff}}$ the number of time points in the interval $[t-\tau_\mathrm{cutoff}, t]$ and $k$ is the complexity of integrating the kernel as a function of the number of required time points. Moreover, these grid points can then also be excluded from future time evolution, which further reduces the overall complexity to $\mathcal{O}\del{n N_{\tau_\mathrm{cutoff}} k(N_{\tau_\mathrm{cutoff}})}$. This point and its relation to the generalised Kadanoff-Baym ansatz \cite{Schluenzen2019} are taken up again in our discussion of the Fermi-Hubbard model in \cite{Meirinhos_2022}.

The algorithm's sensitivity to values off the two-time diagonal can be explicitly adjusted via the parameter \texttt{atol}, irrespective of the nature of the decay of the Green functions away from the diagonal. For rapid (exponential) decay, e.g. in a driven system, a given value of this parameter will lead to a small number of grid points. The same tolerances for slow (algebraic) decay will result in more grid points.