SciPost Submission Page
Adaptive Numerical Solution of KadanoffBaym Equations
by Francisco Meirinhos, Michael Kajan, Johann Kroha, Tim Bode
This Submission thread is now published as
Submission summary
Authors (as registered SciPost users):  Tim Bode · Johann Kroha · Francisco Meirinhos 
Submission information  

Preprint Link:  https://arxiv.org/abs/2110.04793v2 (pdf) 
Code repository:  https://github.com/NonequilibriumDynamics/KadanoffBaym.jl 
Data repository:  https://github.com/NonequilibriumDynamics/KadanoffBaym.jl 
Date accepted:  20220510 
Date submitted:  20220405 08:32 
Submitted by:  Meirinhos, Francisco 
Submitted to:  SciPost Physics 
Ontological classification  

Academic field:  Physics 
Specialties: 

Approaches:  Theoretical, Computational 
Abstract
A timestepping scheme with adaptivity in both the step size and the integration order is presented in the context of nonequilibrium dynamics described via KadanoffBaym equations. The accuracy and effectiveness of the algorithm are analysed by obtaining numerical solutions of exactly solvable models. We find a significant reduction in the number of timesteps compared to fixedstep methods. Due to the at least quadratic scaling of KadanoffBaym equations, reducing the amount of steps can dramatically increase the accessible integration time, opening the door for the study of longtime dynamics in interacting systems. A selection of illustrative examples is provided, among them interacting and open quantum systems as well as classical stochastic processes. An opensource implementation of our algorithm in the scientificcomputing language Julia is made available.
Author comments upon resubmission
List of changes
See "Author comments".
Published as SciPost Phys. Core 5, 030 (2022)
Reports on this Submission
Report
Dear Editor,
I have read the Authors reply to my previous report and the new manuscript. I think the Authors have successfully addressed the points I previously raised and improved their manuscript accordingly. I think this work is now suitable for publication. Given the acceptance criteria of Scipost I think Scipost Physics Core would be more appropriate in this case.
Anonymous Report 2 on 202246 (Invited Report)
 Cite as: Anonymous, Report on arXiv:2110.04793v2, delivered 20220406, doi: 10.21468/SciPost.Report.4880
Report
The authors have substantially improved the paper in the aspects I had commented on. The code has now been published. However, in its current form it will be of little use to the community because it is not enough documented, in particular the examples are not usable in its current form. It would probably increase the reach of this code, the number of people using it and therefore citing this paper, if the authors improved on the documentation.
The central novelty of this paper is the application of adaptive time stepping to KBE, which is a purely algorithmic/numerical progress and no new results on physical phenomena is presented. Therefore, I would still recommend rather publishing this paper in either Codebase or Core.
Francisco Meirinhos on 20220405 [id 2355]
Reply to comments and reported weaknesses
Report 1
We have added a brief theory discussion on how to solve the integral equations arising in nonperturbative treatments of the selfenergies and added an explicit solution of the $T$matrix approximation for the FermiHubbard model.
To underline the advantage of the adaptive method more strongly than conveyed by Tab. 1, we have also added Fig. 11, which shows how the algorithm adapts to a timedependent interaction (quench) in the FermiHubbard model.
We would like to point out that section 4.1.1 does contain a quantitative scaling argument in the sense that we show how the error scales as a function of the tolerances and the order of the integration. For interacting systems, this explicit scaling relative to the exact solution usually cannot be obtained, yet the general argument holds in the same way.
We certainly agree with the referee that solving more involved models is interesting and is the main target of our numerical approach. However, the discussion of such examples would involve substantially more physical discussion of the problems and break the scope of the present paper. Several projects on such systems using the adaptive KB method are presently pursued in our group and will be published in forthcoming papers.
Report 2
We agree that it is not a new development in the integration of ODEs per se, where adaptivity is rather the standard. Yet the additional technical difficulty to implement this for KB equations had to the day precluded reaching this standard here. Hence, we think that solving this problem is an important stepping stone for the physics community working with these equations. Our latest performant results for the Tmatrix of the FermiHubbard model seem to confirm this conclusion.
We understand this point of criticism. The code, alongside a detailed documentation, is now made available.
We prefer not to amend the paper with code documentation because the code may be subject to updates by us or the community, which cannot be included in the manuscript after publication. We are, of course, going to include such an overview in the repository.
Regarding the first point, we refer to our answer to report 1. Regarding the second, we strongly believe that the paper should be published in SciPost Physics, which possesses a suitable audience because it presents a numerical method that will be used by physicists. SciPost Physics Codebases would not reach the wide community of physicists working on the various topics presented in our papers, namely fermions, bosons, classical, open and quantum problems. Regarding SciPost Physics LectureNotes, the method presented is directed to active researchers and, although our paper includes some pedagogical aspects, it is less suitable for beginning researchers in the field. We also mention that the other 2 referees recommend SciPost Physics as well. However, we would be open to the editor's decision for SciPost Physics or SciPost LectureNotes.
We have added a brief theory discussion on how to solve the integral equations arising in nonperturbative treatments of the selfenergies and added an explicit solution of the $T$matrix approximation for the FermiHubbard model.
Not addressed in the resubmission.
Report 3
An error analysis as in Figure 5 and Table I is not performed for all examples because there would not be much more information in treating another exactly solvable quadratic system, while for interacting systems there are usually no analytical solutions available to compare with.
List of changes
Report 1
We have added the following sentence to the introduction: "Such twopoint functions are fundamental objects of manybody physics as they describe, for instance, singleparticle excitations and statistical particle distributions, which represent the essential part of the experimentally accessible observables." and "The distinct nonMarkovian structure of the KB (integrodifferential) equations arises from the reduction of the state space from the (differential) equations generating the MartinSchwinger hierarchy  with Markovian structure but dependence on \textit{all} $n$point functions. Similar to the computation of MartinSchwinger hierarchy, an exact computation of the KB equations is an intractable problem, for which truncations of the interaction diagrams will be required."
We have added the following sentence to section 2.2.2: "Note, however, that the analogy between Eq. [] and a classical action breaks down for interacting quantum systems. For instance, an interaction term such as $H_{int} \sim a^\dagger a^\dagger a a$ results in an equation for the Wigner phasespace distribution with derivatives beyond second order, which is thus not of FokkerPlanck type []. In SchwingerKeldysh field theory, $H_{int}$ in turn leads to nonclassical vertices, e.g. $\phi^{{*}^2} \phi^2$, which are no longer Gaussian in the response fields. Detailed discussions of the differences between a classical approximation neglecting these vertices, which on the level of the selfenergy essentially amounts to dropping contributions from the spectral function, and the full quantum case can be found in Refs. [] for closed systems, i.e. for $\lambda = 0$ in Eq. []."
We have added the following sentence to section 3.3.2: "The clustering decomposition principle~[] ensures that at a largeenough time separation of the physical operators, any $n$point function factorizes. In terms of connected 2point functions this has the signature of an exponential or powerlaw decay in the $\tau$ direction, for massive and massless fields, respectively~[]. This principle should hold for any stable, longlived state, an example being thermalised systems [] described by a Gibbs ensemble."
We have added a short discussion addressing this point: "The sensitivity of the algorithm to values off the twotime diagonal can be explicitly adjusted via the parameter
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. For slow (algebraic) decay, the same tolerances will result in a larger number of grid points."We have moved this part_ "... intraspecies particle transport even in the absence of a direct hopping term between the sites" _of the sentence from the end toward the beginning of the section, and also adjusted the later part.
Regarding the general applicability of our method, we would like to emphasize that it can be applied to any system described by a set of KadanoffBaym (KB) equations. There is per se no obstacle to evaluating nonperturbative KB equations, such as following from $1/N$ or $T$matrix approximations, for instance, on the nonequidistant time grids introduced by our method. This is exemplified by our solution of the Tmatrix for the FermiHubbard model.
Regarding the high mode occupation: in case there is a mapping to an effective classical field theory, resorting to lattice simulation methods is, to our knowledge, the preferred solution. It that sense, our method does not benefit from this since it becomes a different problem. Concerning the classical limit of the quantum field theory, we have added a paragraph to section 2.2.2.
Report 2
We have made the contourGF decomposition (Eq 13) to only depend on $G^>$ and $G^<$. We would like to point out that we stick to the $G^<, G^>$ basis throughout the quantum part, except for the section on open systems where the introduction of (anti) timeordered Green functions is a natural step. Moreover, their relation to the greater and lesser functions is made transparent. For classical stochastic processes, in turn, there is no obvious definition of greater and lesser Green functions, for which reason we fall back on classical analogues of the statistical, retarded and advanced Green functions.
Please see list of changes from Report I. We have included the mentioned reference in 3.3.2 and 4.1.3 and changed: "As discussed in 3.3.2, this observation holds quite generally for quantum manybody systems [65], and is also related to the generalised KadanoffBaym ansatz [49]."
The fact that we refer to a scaling with $n^2$ is because we considered the case where one only truncates the integrals yet keeps on timestepping the twotime plane. The referee is totally correct in pointing out that, once this truncation is applied, it also becomes unnecessary to continue timestepping those points, which then results in the mentioned linear scaling with $n$. We have adjusted our discussion of this in 3.2.2 as follows: "Moreover, these grid points can then also be excluded from future time evolution, which further reduces the overall complexity to $\mathcal{O}\left(n N_\tau k(N_\tau)\right)$."
We have added a brief theory discussion (3.3) and an explicit solution (4.1.3) of the tmatrix approximation for the FermiHubbard model.
We thank the referee for pointing this out and have added this reference to section 2 and 3.
As mentioned in the introduction, KB equations do not scale exponentially with the system size, as other numerical methods do. Rather, the 2point Green functions can always be mapped to $L \times L$ matrices that need to be multiplied at every point in time. The timestepper is agnostic to the operator structure (apart from being by construction dense in the time domain). For the case of dense operators in, e.g., position space, these are implemented efficiently outofthebox by default in Julia (resorting to BLAS calls). The user is free to use other (such as sparse or custom) data structures that appropriately capture the structure of the quantum numbers of the problem in question.
The Github repository has been made public (https://github.com/NonequilibriumDynamics/KadanoffBaym.jl), containing the source code, all examples in notebook form as well as documentation.
We have added a respective discussion to Section 2 (see reply to the first report).
We thank the referee for pointing this out, and have made the necessary adjustments.
Report 3
We have added some clarifications to the roles of the tolerances in the error control of the timestepper.
We have adjusted the caption accordingly.
We have added Fig.[], which shows how the algorithm adapts to a timedependent interaction in the FermiHubbard model, along with the following explanation: "The main strength of our adaptive algorithm can be understood from Fig.[], which shows a prototypical example of a system with varying time scales. The now timedependent interaction parameter $U(t)$ is switched on and off periodically, thus inducing, in particular, different decay rates in the relativetime direction. The resulting effective step size $h$ closely follows the development of $U(t)$, relaxing back to larger values when feasible, while quickly contracting again when the interaction is rapidly ramped up. Thus, when dealing with varying time scales, adaptivity appears to be the natural solution. Note that even though the approximation defined by Eq. (58) is not expected to be accurate in the regimes where $U(t)$ is large, this does not affect the validity of our argument regarding adaptivity. Physical examples where adaptivity could be particularly beneficial are systems displaying prethermalisation [] or the condensation thresholds in photonic condensates out of equilibrium [], where a shorttime evolution with rapid changes is typically followed by a very slow longtime evolution. "