Studying dynamics in two-dimensional quantum lattices using tree tensor network states

1- The manuscript is overall well written and and the results are critically discussed,

2- Finding efficient methods to simulate the real time evolution of D>1 dimensional systems is a timely and important challenge.

1- Lack of references (many of them have already been mentioned by the other two referees). Also the isometric form had been used early on D. Nagaj, E. Farhi, J. Goldstone, P. Shor, and I. Sylvester, Phys. Rev. B 77, 214431 (2008).

2- From the introduction it is not getting completely clear why we should expect tree-tensor networks to perform better than matrix-product states. In particular, they have the same 1D area law restrictions.

3- Given that the manuscript promises to provide a benchmark, a more quantitative comparison would be useful.

In my opinion, the work represents a useful exploration of tree-tensor networks as a tool to study the real time dynamics. While I found the results not too surprising, it is still useful to have a benchmark of the method.

The existing reports already include all my criticisms, and thus I have no additional requests.

The manuscripts propose to use tree tensor networks to perform the time evolution of a 2D quantum system. There could be several good reasons to follow such a procedure but very few of them are carefully reported.

The main weakness is
1) the lack of a critical analysis of their results,
2) the lack of an accurate benchmark of their algorithm (against i.e. ground-state physics)
3) the consideration of only single-site operators
4) the lack of description of the physics of the models considered (e.g. location of the localization transition etc)

I have read the paper with interest and cannot recommend it for publication in its present form.
The authors implement the TDVP for tree tensor network (something that was already done in the literature) and use it to time evolve 2D states.

They constantly mention exact results, but their benchmarks are incomplete, show a strong dependence on the bond dimension of the TTN and there is no criterium the reader can use to establish when the results are reliable.

As a side remark they miss several key references on the topic (including those where the TDVP has been used to optimize 2D TTN), and seem to be completely unaware of the current trends in the TN community where people try to understand how to establish when the results of approximate time evolution are reliable. On this specific aspect, I suggest they read the introduction and references in our recent paper J. Surace, M. Piani, and L. Tagliacozzo Phys. Rev. B 99, 235115 – Published 7 June 2019 where we have tried to give an overview of the effort in the community.

In the following I list the concerns that lead me to formulate the above opinion in detail starting from the analysis of the results presented.

**Regarding the results,

The authors start studying an exactly solvable case of free fermions with disorder.
I am quite puzzled that there is no plot of the exact evolution but only of the error with respect to the exact result.

This does not provide a good insight on how large or small the error is (since it is not even the relative error they plot).
1) They should add the plot of the exact results (even if boring) and superimpose the results of their simulations on the top of it, at least in a inset of their current plots.

2) What convergence within two means in the following sentence?
" Without disorder, convergence of the local density within 2 is obtained for both quaternary and binary TTNS only for t ≤ 1.
They decide to only characterize local expectation values of the fermionic occupation (something pretty restrictive given that they have access to the full many-body wave function). Anyway there is a complete lack of discussion about what are the expected results for such
local observable. I would assume that at least for weak disorder they should equilibrate to a value given by the Gaussian Diagonal Ensemble constructed from their initial state.

3) What is the equilibration time?

4) Is it shorter than the recurrences?

5) When the disorder is ramped up, they enter an Anderson localized phase, what is the critical strength for the transition?

6) What is the effect of such a transition in their evolution? Could they observe it?

7) The values of the used bond dimensions are far too small (in the binary tree they can be pushed to several hundred/ a thousand), they should repeat the simulations with larger bond dimension until lines for the errors with different bond dimension superimpose and then they can attribute all the remaining deviation from the exact evolution to the Trotterization errors.

8) Given the availability of exact results in this context, why haven't the authors tried to extrapolate the finite bond dimension results and see if they are able to obtain the exact ones?

9) I am also concerned with the lack of fundamental checks on their evolution. For example the evolution is unitary and should conserve the energy of the initial state. Actually the deviation of the conservation of the energy has been widely used as a measure on how good an approximate evolution is.
In the present version there are no plots about the energy conservation, and the lack of discussion about the use or not of a symplectic integrator and if and for how long the energy is conserved.
 
10) As a result the plots presented are not really conclusive, they show that if you run the tdvp on tree tensor network you get some results out of it. But what these results mean is far from obvious, at least for me. How does this approach compare with say a simple spin-wave analysis or semiclassical approaches like the truncated Wigner methods?

11) Using the Jordan Wigner transformation for simulating fermions is definitely a possibility but I would suggest at least to mention all the theoretical development about fermionic tensor networks.

12) The sentence  ."..The application of TDVP formally requires the TTNS corresponding to the initial condition to possess a full tree rank of r." is formally wrong. The fact that TDVP works on subspaces with fixed bond dimension has been addressed already in the literature by passing to a two sites algorithm (and some initialization stage). I believe that adding noise to the initial state in order to fill the rank is a pretty dangerous strategy, that could lead to unexpected results
My belive seems to be confirmed by the sentence
"For times beyond the convergence time we also observe a strong influence on the initialization procedure
of the redundant parameters of the TTNS. Increasing the bond dimension systematically reduces this effect."

 13) What is the convergence time they talk about?
All the above observations are even more relevant once they study the non-solvable system of hard-core bosons. Here we are completely lost and do not know how to address the results.

14) The quantity rho discussed in the text is defined in a label of fig. 5 with a site index (already pretty strange, I would have expected a definition in the main text). Furthermore in the definition of the anisotropy A \rho is defined with two indices that are summed over (so they supposedly represent the two coordinates of a point). However it is not clear to me where the zero of coordinates is and thus which asymmetry A is measuring.
I guess it is the crossing of the two lines in fig 5 but the reader should not guess.

15) Furthermore there is no explicit indication about how physical sites in the lattice are mapped to the TTN, since one could do it in infinitely many ways. I guess that the obvious choice has been made, but it would be important to understand if the 0,0 point is in between the tensors or in the center of a 4 leg tensor, in order to understand better the role of the structure of the network on this geometry.

16) I am very puzzled by the sentence "However, given the systematic improvement in the convergence with respect to the bond dimension combined with reduction of the anisotropy of
the results, we consider the latter to achieve numerically exact results for t < 2." 
Their plots show huge variations between TTN with different bond dimensions and between TTN with different structures, what is exact from their point of view?

***General mistakes and lack of references
First of all there are several errors in the text.

1) In the introduction they claim "the logarithm of the latter quantity gives an upper bound on the entanglement entropy for every bipartition of the lattice" this is simply wrong and should be reformulated correctly. (Think of a bipartition in an MPS when one separates even from odd sites)

2) Saying that PEPS algorithms involve uncontrolled approximations seems to me a bit too strong. But this is actually an opinion.

3) There is a typo "however describing to higher spatial dimension"

4) In the list of reference about MPS with long-range interactions key references are missing, in particular, the foundational papers are as far as I am aware

Gregory M. Crosswhite, A. C. Doherty, and Guifré Vidal
Phys. Rev. B 78, 035116 – Published 14 July 2008

Fröwis, F.;Nebendahl, V.;Dür, W.
Physical Review A, vol. 81, Issue 6, id. 062337

5) I have co-authored two papers in this context that deal with 1) the TDVP for systems with long-range interaction, in a version where the MPS is interpreted as a TTN (as you can appreciate in the drawings in the supplementary material). The algorithm is characterized both at equilibrium in  Koffel, Lewenstein, Tagliacozzo, Physical Review Letters, vol. 109, Issue 26, id. 267203
and out of equilibrium in  Hauke Tagliacozzo  Physical Review Letters, vol. 111, Issue 20, id. 207202.
In the supplementary material of both papers we have introduced a version of the single-site TDVP these two papers should be cited in conjunction with References 16 and 36 since they both were published before the two cited papers.

6) The TDVP for studying ground states of 2D TTN was implemented and benchmarked in Andrew J. Ferris Phys. Rev. B 87, 125139 – Published 25 March 2013  the paper should also be referenced.

7) As a matter of fact, I strongly recommend using the TDVP they implement to first find the ground state of 2D systems described with TTNs (by just moving to imaginary time) and compare the results they obtain with the one available in the literature about bond-dimensions and precision (e.g. the above paper is a good starting point but they could also compare with the results of their reference 23).

Once they are sure that their algorithm does not contain any error, they can use it to time evolve the systems and perform accurate comparisons with free systems (as indicated in the first part of this report) and then present the results for the interacting cases.

Alternatively, if they do not want to compare their algorithm with known results at equilibrium, they should at least compare it carefully with the exact results, where available, and run extensive tests against exact diagonalization results on small lattices where the analytical results are not available.

They are already listed in detail in the previous

1- Tackling time-evolution of two-dimensional systems is an open problem, where tentatives have been done, but no standard recipe is available yet -- therefore definitely with topic to be further investigated

2- The manuscript is self-contained, with a large didactic introduction of the whole notation (see however weak point 3) and even a pseudo-code snippet

3- The chosen benchmarks are fully reasonable: an exactly solvable problem (non-trivial for tensor-network (TN) methods, though), and another system recently tackled in similar studies

4- The overall structure of the manuscript and the aims of the work are clear, and there is no apparent trace of over-selling results

1- The analysis of the benchmark results is not clear in some aspects, and seems performed in a bit of a rush -- see Report for Details

2- Results are presented for two different tree-tensor-network (TTN) structures, but the "quaternary" is not mentioned or sketched in the whole text (though it is easy to imagine what is that...)!

3- A couple more illustrations in Section II, where quite some definitions and symbols are introduced, would be desirable: without knowing Ref. 36 almost by heart, it would be quite though for the reader to reconstruct the derivation of the method.

As highlighted in the list of "Strengths", the manuscript deals with an interesting problem, the description of time-evolution in two-dimensional quantum systems via tensor-network methods, for which a solution is far from being fully established.
Moreover, it illustrates the method and the results of sensible benchmarks in a self-contained manuscript.
Therefore it definitely has the potential for being published.
However, as summarised under "Weaknesses", there are a number of reasons that prevent its publication in its present form and demand a further iteration by the Authors.
After having addressed them, I would be in favour of publication.

-----
In general, the manuscript seems to give proper credit to other contributions and to state similarities and differences with works dealing with a very similar problem (11-21 & 35-41). Only I was a bit surprised not to find some more works about trying to extend the limitations of time-evolution, e.g.,
C. Krumnow, L. Veis, Ö. Legeza, and J. Eisert, Phys. Rev. Lett. 117, 210402 (2016);
C. Krumnow, J. Eisert, and Ö. Legeza, arXiv:1904.11999;
M.M. Rams and M. Zwolak, Phys. Rev. Lett. 124, 137701 (2020)
and under the ones listed about "TTNS are rarely used in the context of interacting lattice systems", e.g.,
V. Murg, F. Verstraete, Ö. Legeza, and R. M. Noack, Phys. Rev. B 82, 205105 (2010)
W. Li, J. von Delft, and T. Xiang, Phys. Rev. B 86, 195137 (2012)
M. Gerster, et al., Phys. Rev. B 96, 195123 (2017)
E. Macaluso, et al., Phys. Rev. Research 2, 013145 (2020)
Yes, the latter are own references, sorry about that, but they deal with 2D strongly correlated systems and thus fit particularly well in such a list. A more extensive search in the literature is definitely welcome.

-----
Some statements could probably be sharpened: e.g., "Since the entanglement of a generic system grows linearly with time, the accessible timescales are limited" should be phrased in terms of global quenches, and the distinction with respect to local ones made, etc.
By the way, Ref. 1 was among the first ones providing numerical evidence, but there are theoretical predictions of the linear growth in 1D, too:
e.g., P. Calabrese and J. Cardy, J. Stat. Mech.: Theory and Experiment 2005, P04010 (2005) and Journal of Statistical Mechanics: Theory and Experiment 2007, P10004 (2007) (the second one focussing on the mentioned distinction between local and global quenches).

-----
A graphical representation of Eqs. 5-8 would be desirable for the readers not immediately recalling the whole Ref.36 by heart
Incidentally, a graphically more polished version of Figs. 1-3 could be nice to see (though only decorative, I admit).

-----
About the time evolution procedure itself:
i) Should Eq. (12) not contain a minus prefactor with respect to Eq. (11), according to Eq.(5) and to what happens for the MPS case in Ref.36?
ii) Below Eq. 12 it is stated that, in general, "the result of the contraction of the Hamiltonian with the environment tensors is not a set of matrices but a compact tensor network". Do the Authors mean, as it should be the case, that a MPO/TPO Hamiltonian would result in a MPO/TPO connection of the three terms in Fig.3b? Why then not illustrate it in Fig. 3 directly, and simplify considerably the discussion?
iii) I must admit that in Sec.II.C-D I do not see the difference between the adopted "splitting integrator" and "propagation carried out while descending and ascending on the tree": Could the Authors clarify their point, possibly via an illustration?
iv) What supports the statement "it is unlikely that the quadratic error bound in the time-step for the total time-evolution of the latter carries over"? It would be highly desirable to substantiate this with a quantitative analysis in the following, which I cannot find.
v) Same applies for "This intuition, however, is largely based on systems with local, or at least smoothly decaying interactions, and can break down for interactions which result from mapping a 2D lattice to a 1D chain, for example."
vi) I find the remarks about the "redundant parameters" scheme interesting, but how does this procedure compare with a two-tensor optimization, or a local inflation of the bond dimension in a one-tensor scheme, as it is done for optimization purposes already? [see White, Phys. Rev. B 72, 180403 (2005), Hubig, et al., Phys. Rev. B 91, 155115 (2015); Silvi, et al., SciPost Phys. Lect. Notes 8 (2019)]

-----
Detailed remarks about the analysis of the benchmark:
a) The role of the time-step is not discussed and clearly separated from the approximation due to bond-dimension truncation.
b) The statements about better reliability of the binary structure with respect to the quaternary one do not seem to be supported from the plots in Fig. 6, unless I am overlooking something very naive. In particular, I really do not see why the Authors say that data from the binary structure display a "systematic improvement wrt bond dimension", implying that the ones from quaternary do not...
c) An hint about the relative error, not only its absolute value, would be desirable; even better would be the typical plot of the exact (or best) solution vs all other approximate curves as a function of time: I mean in the style that is typically used to illustrate the explosion of errors in a 1D setup -- this would probably help about point a), too. At this level, and in view of the point iv) above, it is not clear that the choice of integration method and time-step does not play an important role, too...
d) The results for non-interacting fermions are all for an initial random state (by the way, at which average filling?): What would be different for a more structured initial state? Say a cluster away from or around the sample boundaries, a checkerboard pattern, etc.? At least one would naively expect to see a comparison between fermions and hard-core bosons in the same style of Fig.5, not?
e) Maybe a too ambitious question, but would other kind of transformations, like Bravyi-Kitaev instead of Jordan-Wigner, impact on the performances of the algorithm?
f) Incidentally, are symmetries, namely particle-number conservation, encoded in the TTNs structure?
g) Why do the Authors not ask the raw data of Ref.16 for a direct comparison? This would be of great interest, as stated somewhere in the text, too.

-----
Minor points:
- At line 5 of Algorithm 2, it should read t_{1/2} instead of t_0, right?
- What does "convergence of the local density within 2" mean at page 6, column 1?
- A part of the caption of Fig.4 got probably lost: it should read similar to Fig.6, in which the reference to quaternary TTNs is clear (this should be added in the main text, too).
- When comparing the best binary and quaternary data, the Authors seem to infer a lot from their difference being of order 1%: where exactly is the surprise, given that both have an error estimate of roughly that order? Can the Authors find a better illustration to convey their message?
- Alongside with refs. 36 and 41 on the TDVP method, I would have expected to see also J. Haegeman, et al., Phys. Rev. Lett. 107, 70601 (2011): is there a particular reason behind the choice of omitting it?
- I am honestly puzzled by "For times beyond the convergence time we also observe a strong influence ...": are data of any meaning beyond that point? Maybe it is a mismatch in what I understood, but the term "convergence time" is not defined anywhere...
- What do the Authors exactly intend with dynamics of "critical one-dim. systems"? What is going to vary and what is quenched or tuned?
- By the way, there are good reasons for a stable time-evolution algorithm be missing for MERAs (though notice Phys. Rev. A 77, 052328 (2008)), given that their structure is full of loops and already extremely hard to manipulate for static optimisations...