Adaptive Numerical Solution of Kadanoff-Baym Equations

Reply to comments and reported weaknesses

Report 1

Seems to be limited to special classes of models, where the interaction can be treated perturbatively (...)

We have added a brief theory discussion on how to solve the integral equations arising in non-perturbative treatments of the self-energies and added an explicit solution of the $T$-matrix approximation for the Fermi-Hubbard model.

Worked examples not clearly demonstrating the advantage of the new method (...) A more precise / extensive discussion of the advantage obtained from the adaptive method would be... a service to practitioners... The only statement I could find about that comparison is in the first example but is very short and qualitative. Are there not scaling arguments available?

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 time-dependent interaction (quench) in the Fermi-Hubbard 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.

the paper would benefit very much from more impressive examples.

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

Using adaptive step sizes, while technically challenging to implement, is by no means a new development in the integration of differential equations.

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 T-matrix of the Fermi-Hubbard model seem to confirm this conclusion.

However, due to the unavailability of the advertised GitHub repository and the lack of a deeper documentation of the alleged repository

We understand this point of criticism. The code, alongside a detailed documentation, is now made available.

However, if the paper was amended with a (short) overview of the julia code

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.

No new physics results, advance is purely numerical (...) I believe that there is too little physics advancement to meet the requirements of SciPost Physics.

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.

Only perturbative approximations, no memory integrals in self-energy evaluation considered

We have added a brief theory discussion on how to solve the integral equations arising in non-perturbative treatments of the self-energies and added an explicit solution of the $T$-matrix approximation for the Fermi-Hubbard model.

Only small systems studied

Not addressed in the re-submission.

Report 3

Figure 5 and Table I... Furthermore, one wonders why a similar analysis is not performed in the other cases

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

Why should it be enough to study the two-time correlator?

We have added the following sentence to the introduction: "Such two-point functions are fundamental objects of many-body physics as they describe, for instance, single-particle excitations and statistical particle distributions, which represent the essential part of the experimentally accessible observables." and "The distinct non-Markovian structure of the KB (integro-differential) equations arises from the reduction of the state space from the (differential) equations generating the Martin-Schwinger hierarchy -- with Markovian structure but dependence on \textit{all} $n$-point functions. Similar to the computation of Martin-Schwinger hierarchy, an exact computation of the KB equations is an intractable problem, for which truncations of the interaction diagrams will be required."

Reduction Schwinger-Keldysh - MSR: a word of caution should be placed in what concerns this reduction in the case of interacting problems

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 phase-space distribution with derivatives beyond second order, which is thus not of Fokker-Planck type []. In Schwinger-Keldysh field theory, $H_{int}$ in turn leads to non-classical 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 self-energy 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. []."

Memory truncation: should be re-emphasized that the Green's function are connected, therefore one has clustering properties

We have added the following sentence to section 3.3.2: "The clustering decomposition principle~[] ensures that at a large-enough time separation of the physical operators, any $n$-point function factorizes. In terms of connected 2-point functions this has the signature of an exponential or power-law decay in the $\tau$ direction, for massive and massless fields, respectively~[]. This principle should hold for any stable, long-lived state, an example being thermalised systems [] described by a Gibbs ensemble."

Interesting would be a discussion of the difference between algebraic and exponential decay in that coordinate, how does this affect the scaling of the method?

We have added a short discussion addressing this point: "The sensitivity of the algorithm to values off the two-time 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."

What is the reason to choose the excitation transfer model? the logic of choice of models should be anticipated.

We have moved this part_ "... intra-species 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.

The applicability of the method to generic situations (i.e. beyond perturbation theory in the interactions) is not clear, and not even an effort is made to comment on that. This is key in non-trivial out-of-equilibrium dynamics, such as the description of non-thermal fixed points. Some comparison to the so far successfully developed techniques to tackle such problems would be useful. An answer could be that the latter rely on an effectively classical description due to high mode occupations. Could the present technique also benefit from that?

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 Kadanoff-Baym (KB) equations. There is per se no obstacle to evaluating non-perturbative KB equations, such as following from $1/N$ or $T$-matrix approximations, for instance, on the non-equidistant time grids introduced by our method. This is exemplified by our solution of the T-matrix for the Fermi-Hubbard 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

Unify the basis for the Green's functions (GF)

We have made the contour-GF 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-) time-ordered 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.

When discussing the fact that GF decay as a function of relative time, the GKBA is often referenced as an interpretation. To me, that's just another numerical method, not a physical interpretation. The reason that GFs decay as a function of rel. time is the fact that closed, interacting quantum systems thermalize, which implies a loss of memory of the initial state. See e.g. Berges, Cox, Physics Letters B 2001. This could be emphasized.

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 many-body systems [65], and is also related to the generalised Kadanoff-Baym ansatz [49]."

I'm not sure I understand why in section 3.3.2. it is said that the numerical effort with a memory cut scales as $n^2$.

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 time-stepping the two-time plane. The referee is totally correct in pointing out that, once this truncation is applied, it also becomes unnecessary to continue time-stepping 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)$."

_ In the intro, 1/N expansions are mentioned (...) it would be great to at least comment if not explicitly calculate how the adaptive algorithm can be generalized to this case._

We have added a brief theory discussion (3.3) and an explicit solution (4.1.3) of the t-matrix approximation for the Fermi-Hubbard model.

The excellent review of Schlünzen et al...

We thank the referee for pointing this out and have added this reference to section 2 and 3.

The efficiency of the implementation also depends on the scaling with system size L. It would be great to show that it is also optimal, i.e. the various matrix/tensor contractions were implemented well.

As mentioned in the introduction, KB equations do not scale exponentially with the system size, as other numerical methods do. Rather, the 2-point 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 out-of-the-box 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 repo should be made public before publication and a short documentation

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.

A central point of the paper is the connection between quantum time evolution and classical stochastic evolution. There is however another limit, namely the classical limit...

We have added a respective discussion to Section 2 (see reply to the first report).

I'm not sure Eq.(6,10) is consistent with Eq.(1). Shouldn't there be a factor of 2 if the Grassmann Fermionic fields are decomposed into two real Majoranas?

We thank the referee for pointing this out, and have made the necessary adjustments.

Report 3

I suggest the authors to expand and clarify the discussion around Eq.39.

We have added some clarifications to the roles of the tolerances in the error control of the timestepper.

Figure 5 and Table I...: the table itself is rather hard to read without an appropriate caption. Various symbols used there have been introduced in previous sections and it is hard to keep up for the reader.

We have adjusted the caption accordingly.

Overall the question of whether this approach should be preferred to more conventional ones need to be addressed more directly.

We have added Fig.[], which shows how the algorithm adapts to a time-dependent interaction in the Fermi-Hubbard 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 time-dependent interaction parameter $U(t)$ is switched on and off periodically, thus inducing, in particular, different decay rates in the relative-time 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 short-time evolution with rapid changes is typically followed by a very slow long-time evolution. "