Reply to Referee

The referee writes:

(1) The paper aims to study the real-time dynamics of quantum impurity problems. The bath is cast into an effective description, allowing the authors to treat the impurity in an effective environment. Quantum impurity models represent a long-standing important class of problems in condensed matter. Hence it is interesting to see new developments. The authors, however, already have three other papers along this line already within the last half a year, all with very similar titles: [41] R. Chen, X. Xu, and C. Guo, Grassmann time-evolving matrix product operators for quantum impurity models, PRB (2024). [42] R. Chen, X. Xu, and C. Guo, Grassmann time-evolving matrix product operators for equilibrium quantum impurity problems, New J. Phys. (2024). [43] R. Chen, X. Xu, and C. Guo, Real-time impurity solver using grassmann time-evolving matrix product operators, arXiv:2401.04880 (2024). This raises the question of whether the present paper represents a rather incremental or even duplicate step. Having several similar near-simultaneous papers and referencing these, makes the present paper rather unreadable by itself. Figures like Fig. 1 are not understandable the way they are currently presented.

Our response: We thank the Referee for refereeing this work and for his/her comments. The previous works focus on impurity solvers using our newly propose GTEMPO framework. Concretely, Ref.[41] solves non-equilibrium impurity problem on the real-time axis. Ref.[42] solves the equilibrium impurity problem on the imaginary-time axis. Ref.[43] extends the method in Ref.[41] to solve equilibrium impurity problem on the real-time axis with a quench. In those works, a central numerical calculation is to build the influence functional as a GMPS, which is done using the partial-IF method. The number of GMPS multiplications required in the partial-IF algorithm scales as O(N) for N discrete time steps. The current work does not aim for any new impurity solver. Instead, it aims to replace the partial-IF algorithm with a more efficient algorithm, referred to as the TTI-IF method, to build the influence function as a GMPS. Therefore, one could use the TTI-IF method in all the previous works to significantly improve the numerical efficiency. Compared to the partial-IF algorithm, the TTI-IF method respects the time-translational invariance of the influence functional, and the number of GMPS multiplications required in the new method does not scale with N. As we understand from the referee’s comment, the difficulty of understanding Fig.1 is mainly due to that the GMPS multiplication is not well explained in this work. The multiplication of two GMPSs result in a new GMPS, which is different from usual MPO or MPS arithmetic. This operation is introduce in Ref.[41], to make the current work more self-contained, we have added a short explanation of this operation in Appendix. A. Roughly speaking, the GMPS multiplication originates from the Grassmann tensor multiplication, which is the only operation required to build the MPS-IF. We also mention that the bosonic counterpart of the Grassmann tensor multiplication is the element-wise product of two normal tensors. In principle, one could also convert the MPSs into MPOs by “copying” its physical indices, and the GMPS multiplication will be converted into MPO multiplications which is a standard operation. But we prefer to directly implement the GMPS multiplication for efficiency without the intermediate conversion step (the element-wise operation is more natural for manipulating the path integral).

The referee writes:

(2) The paper is aware of other recent works on MPS representation of the Feynman-Vernon influence functional (Refs. 47-51). However, the precise overlaps of and differences with these other approaches are not discussed much at all. A few lines in the introduction would be helpful. It appears there is a distinction based on a many-body treatment vs. Gaussian-like states which effectively can work in a single-particle picture. This needs to be explained in more detail. Hence while the present paper appears to have merit, I strongly recommend that the authors address in detail the points raised below in the paper itself prior to me deciding on a recommendation for its publication.

Our response: We thank the referee for this comment. We have added in the introductory section a brief discussion about the difference between our GTEMPO and the tensor network IF method, which is attached as follows:

Formalism-wise, these methods differ from GTEMPO in that in their numerical calculations the PI is converted into a fermionic operator expression in the Fock state basis, thus avoiding directly dealing with Grassmann variables (GVs). In the meantime, the algorithm design in GTEMPO could be more straightforward as it directly translates the Grassmann expression of the fermionic PI into MPS calculations.

Briefly speaking, in the tensor network IF method, the Grassmann path integral is converted into a fermionic operation expression, which thus avoids to directly deal with Grassmann variable numerically. However, the addition transformation would make it less convenient for numerical implementation. Moreover, as shown in Ref.[41], for transport problem the computational cost of the tensor network IF method scales exponentially against the number of baths, while the cost of GTEMPO is essentially independent of the number of baths.

The referee writes:

(3) More detailed remarks: Overall, I find it difficult to follow the arguments. In Eq. (5): why do the partial IFs commute with each other? After all, k is summed over, and the a's anticommute. Where does Eq. (7) come from?

Our response: We thank the referee for this comment. This is because that any Grassmann expression with an even number of Grassmann variables commute with each other, which is added in the end of Page 4 in the revised manuscript. Thinking of Eq.(7) requires some intelligence and experience. However, it is easy to verify that it is correct: One can simply do the multiplication of the right hand side of Eq.(7) and will find that the result is exactly Eq.(6).

The referee writes:

(4) Fig. 1 is not understandable without having to look up other references of the author, like Refs 41-43. And even then I am not sure what is shown. Naively, MPS cannot be muliplied unless one calculates an overlap. But multiplying them to get another MPS? That usually requires MPOs. Do the authors have tensor product in mind? Also, how does this relate to a Keldysh contour?

Our response: We thank the referee for this comment. As in our answer to comment (1), we think that the gap in understanding Fig.1 is to understand the definition of GMPS multiplication, which results in a new GMPS. We have thus added in Appendix.A the definition of GMPS multiplication, which is the only operation needed to build the MPS-IF.

The referee writes:

(5) The bath of an impurity is described by a continuum of states. Hence from a numerical perspective, I'd expect this needs to be discretized, only then resulting in discrete Grassmann variables. But this also comes with a discretization error that appears discussed nowhere. This seems hidden behind the reference to the QuaPI'method(bess,shod––thisratherbe𝑄𝑢𝑎𝑃𝐼′𝑚𝑒𝑡ℎ𝑜𝑑(𝑏𝑒𝑠𝑠,𝑠ℎ𝑜𝑑̲𝑡ℎ𝑖𝑠𝑟𝑎𝑡ℎ𝑒𝑟𝑏𝑒QuAPI'?) For ED, the paper mentions uniform discretized of the bath into 8000 intervals. How does this compare to GTEMPO? If there is a coarse-graining of the bath in GTEMPO, does this coarse-graining have to be uniform over the spectral range? Or may one choose a logarthmic grid as common for the numerical renormalization group (NRG) for the benefit of reaching much lower energies / longer time scales? Note that there are also mixed log-linear discretization schemes used in the literature for non-equilibrium transport through correlated impurities.

Our response: We thank the referee for this comment. We would like to clarify some misunderstanding of the referee about the GTEMPO method here. First of all, the major difference of GTEMPO, compared to conventional MPS methods, is that the bath degrees of freedom is analytically integrated out in GTEMPO using the Feynmann-Vernon influence functional, and one is only left with the impurity degrees of freedom at different times, which is represented as a GMPS. Thus there is no bath discretization error in GTEMPO, the major sources of error are the time discretization and the MPS bond truncation error. Second, the QuAPI method does not handle the bath discretization neither. It merely discretize the continuous hybridization function \Delta into a discrete hybridization matrix used in the exponent of Eq.(4). We also thanks the referee for pointing out the mistake of abbreviation QuaPI to QuAPI, we have corrected it in the article.

The referee writes:

(6) The paper writes: `With the more efficient TTI IF method to construct the MPS-IF, we can easily reach longer evolution time' [up to t=8.4/Gamma]. For impurity models which can have an exponential dynamically generated low-energy scale (the Kondo temperature) this is still a very modest time range. What is the limiting case in time here?

Our response: We thank the referee for this comment. The computational cost of the TTI-IF method scales as O(N\chi^4) as pointed in the manuscript. For t=8.4/\Gamma with \delta t=0.014/\Gamma, we have N=600. And since there are 8 Grassmann variables per time step, the size of the GMPS is about 4800 already in our numerical simulations. It is easily to do better than the partial-IF but the cost will still be very significant if we further enlarge t. Therefore we stop at t=8.4/\Gamma. In the revised manuscript we have changed the sentence “we can easily reach longer evolution time” into “we can reach longer evolution time” to weaken the statement.

The referee writes:

(7) The paper talks about steady-state currents, but Fig. 4 still mostly shows transient behavior. If t=8.4/Gamma `can be easily reached', why not go to much longer times to reach a steady state? It would be helpful to quantitatively compare to steady-state currents for interacting quantum impurity models as found already in the literature, e.g., PRL 101, 140601 (2008), or PRL 121, 137702 (2018).

Our response: We thank the referee for this comment. The reason that we do not go to even larger times is answered in our reply to comment (6): there is a linear scaling of the cost with t, and even with current choice of t the size of GMPS is close to 5000, thus it is hard to consider a much larger t and we decide to stop at this value. We also note the numerical simulation in Fig.4 does not aim particularly for the steady state. In fact, one would still need algorithmic development if one is interested in the steady state only, since the time to reach steady state may be very long with the current method, which can be seen from Fig.4(a,b,c), especially for small chemical potential bias. For example, one could use infinite MPS techniques which directly targets at the steady state and completely neglects the transient dynamics (our later work, arXiv: 2403.16700 extends the technique of the current work to directly aim for the steady state). Similarly, benchmarking the steady state currents with “PRL 101, 140601 (2008), or PRL 121, 137702 (2018)” is likely to be very challenging for the current method, which is left to investigate in future developments. To stress the scope of our numerical simulation of the transport problem, we have added the sentence “Our goal in this study is mainly to verify whether the existing finite-time calculations have reached steady state or not, by extending the evolution time with our new method” in the revised manuscript. From Fig.4, we can clearly see that the answer to this question is no for small V < 0.71 and is yes otherwise.

The referee writes:

(8) Also why not show the full transient in Fig. 4 from t=[0,8]/Gamma rather than only the second half [4,8]?

Our response: We thank the referee for this comment. We start from t=4.2/\Gamma, as the existing finite-time calculations stop here, and our goal is to verify whether these existing finite-time calculations has reached steady state or not. We have added this sentence in the revised manuscript, which is attached as follows:

(which is the longest time that has been reached in previous studies [41, 48])

The referee writes:

(9) Minor details: The paper writes: `we note that our method can be directly applied to general QIMs as long as the Feynman-Vernon IF applies'. Please be more specific: When does it not apply? Does this require a non-interacting bath?

Our response: We thank the referee for this comment. We have added an explanation after this sentence, which is attached in the following:

we note that our method can be directly applied to general QIMs as long as the Feynman-Vernon IF applies, e.g., the bath is noninteracting and is linearly coupled to the impurity.

The referee writes:

(10) After Eq. (4), the `bath spectrum density', by its name, appears a bath property, and therefore should be independent of impurity parameters such as Vk which already refers to hybridization.

Our response: We thank the referee for this comment. The name of this function differs from literature to literature. We agree that the name “bath spectrum density” is misleading as it not only contrains the bath information. We have replaced it by “coupling strength function”, which has been used in [Erpenbeck et. al., Phys. Rev. Lett. 130, 186301 (2023)].

The referee writes:

(11) It appears to me that the Prony algorithm is closely related to `linear prediction'. The authors may briefly comment on the relation between the two in their Appendix A.

Our response: We thank the referee for this comment. These two algorithms are indeed intimately connected: both of them learn the function in Eq.(B4). The difference is that in the Prony algorithm Eq.(B4) is explicitly used (the parameters in this equation is directly used to construct the MPO afterwards), while in linear prediction, one often does not need the explicit form of Eq.(B4), but only the recursion relation between the history and future, which is implicitly determined by Eq.(B4). Of course, in linear prediction, one can also obtain the explicit form of Eq.(B4) first and then use it to predict the future, instead of predicting the future based on a recursion from the history data. We have added this discussion in the end of Appendix.B in the revised manuscript.