Automatic differentiation applied to excitations with Projected Entangled Pair States

-Detailed description of the methodology.

-Just a simple benchmark is reported.
-Some issues on the thermodynamic extrapolations to be fixed.

In this work, the authors develop an efficient method to optimize tensor networks for excited states, with an application to the single-band Hubbard model.

The paper is interesting and well written, with a detailed despription of the optimization part.

I only have a few comments:

1) I do not like (and fully understand) Fig.1. First of all a single inset for two extrapolations (U/t=4 and 12) is a bit confusing, especially
because the thermodynamic value for U/t=12 is on the left, while the labels are on the right.
In addition, $\omega_s$ is not defined in the main text, but only on the appendix. Most importantly, I suspect that the shift in the main panel
is -U/2 and not -U. Otherwise, the values in the inset are not compatible with the ones in the main panel.

2) At page 13, it is written that "the gap can be extracted in a controlled way" for iPEPS simulations. However, in Fig.1 the value of the gap strongly bends down for D=6 and the extrapolated value could be as small as 3.5 or so. I suspect that D=7 or 8 are beyond the present possibilities, but it would be nice to have a fair discussion on the extrapolation issues.

3) Related to the previous points: are the iPEPS gaps reported for the thermodynamic extrapolation or for the largest D values? The thermodynamic extrapolation should have much larger errorbars.

4) I should not understand why in Eq.(23) the upper index "1p" is used. When inserting the complete basis of single-particle states, no approximations are done for the spectral function of Eq.(22).

5) Reference [11] has been published on SciPost Physics in 2021.

In summary, the paper can be published after the authors have addressed these points.

See report.

1) Very clear and (almost) self-contained presentation of the developed method and how it can be implemented.

The authors explain how one can apply automatic differentiation (AD) techniques to optimize 2D tensor networks (PEPS) for the recently introduced excitation ansatz for quasi particles on top of ground state. They benchmark and compare their results for the 2D Hubbard model against unbiased quantum monte carlo (QMC) results and demonstrate that they can comfortably go into the regime $U\gg t$ inaccessible with QMC.

With its detailed explanations, the manuscript enables a broader use of the PEPS-based excitation ansatz by the community. The numerical exploration of spectral properties of strongly-correlated 2D systems is generically very challenging, such that the described technique has the potential to become a state-of-the-art method; further follow-up work based on this technique is to be expected. Hence, the presented paper satisfies the acceptence criteria for Scipost Physics, and I recommend a publication once the following points have been addressed.

1) The explanation of the CTM method is almost self-contained, but it lacks the somewhat important detail of how the projectors to smaller bond dimension around equ. (12) are obtained. Looking into the provided example code reveals that it's basically just an SVD, but I think it would be beneficial to briefly comment on that in the main text as well.
Of course, it's the authors decision to which extent they want explain it.

2) I don't understand the pseudo code of Algorithm 2.
Line 25 uses apply_H_eff, but only apply_H_eff_fixed_point is defined. On the other hand, the converge_boundaries function is unused.
Do you call converge_boundaries at the beginning of apply_H_eff_fixed_point instead of taking C as an extra argument, and use apply_H_eff_fixed_point for the eigensolver?

3) My understanding is that you use an eigensolver that obtains multiple (orthogonal) exciations for each k point at once, correct? Could you maybe stress this in the text?
How many excitations did you actually calculate to plot Fig. 3?

4) Upon the first read, it wasn't immeditately clear to me that the $\omega_S$ in Fig. 3 refers to the energy of the excation ansatz at momentum S=(pi/2,pi/2), because the labeled path through the BZ is only introduced later in Fig. 3.
Maybe it could be clarified that you know/expect that the excitation energy is minimized at the S point due to the AFM order.

5) I think the term "complete basis of single-particle states" is a bit misleading/confusing/not strictly correct.
Really, the basis you use are the excitation ansatz states $|\Phi(B_A)_k>$ on top of the ground state, as becomes evident when you plug them in in eq. 25.
These states are still strongly-corrleated many-body wave functions, and they don't form a complete basis - as you also mention in the discussion/comparison with QMC, the PEPS calculations do miss out the higher-energy continuum part of the spectral function, and this is the origin of it.

6) Typos:
- ylabel of Fig. 2: comparing with Fig. 3, I think this should also be $(\omega_S - U/2)/t$.
- caption of Fig. 3: $S(\pi/2,\pi,2) \rightarrow S(\pi/2,\pi/2)$

1) Development of a generic and efficient PEPS based method for computing excitation spectra of strongly correlated systems in 2D
2) Method is clearly presented

1) Somewhat limited discussion of the spectral functions of the 2D Hubbard model, that are evalu-ated with the proposed technique

In this manuscript, the authors introduce an excitation ansatz for PEPS states that are optimized using automated differentiation (AD) techniques. AD-based techniques have already been successfully used to optimize iPEPS ground states. Here, AD is applied to optimize an excitation ansatz to obtain spectral information of strongly correlated 2D systems; a notoriously challenging task. The proposed method is used to compute the single-particle spectra of the 2D Hubbard model at half filling. The spectra are benchmarked against Quantum Monte Carlo results obtained using Wick rotation techniques (stochastic Max Ent methods).

Overall, the paper is clearly written and the developed computational method is well presented. The discussion of the results is a bit short and could profit from a deeper connection to the existing literature. Here are some suggestion and questions that seem interesting:
1) Could one use rotational symmetry of the state to avoid computing horizontal and vertical optimizations of Eq. 11? In the later example of the Hubbard model a larger unit cell is used, but similar considerations could be transferred to this case too.
2) On page 13 it is briefly discussed why it is hard for QMC to go to large U/t values. As it is written now, it is not immediately clear what is meant without consulting the appendix. It could be beneficial to improve this discussion in the main text.
3) Going from Eq 22 to Eq 23 is exact. It is suggested otherwise in the paragraph in between these two equations. This part should be clarified. Maybe, the convergence generating factor could already be introduced here too.
4) How accurate is the fit to the spin-polaron dispersion in Fig 3? It seems that at the Gamma point, where faint spectral weight is visible in the ground state branch, the spin-polaron dispersion fit strongly overestimates the energy (spin-polaron: w-U/2 ~ -4.5t vs numerics: -3.8t). The authors should discuss implications. Is the proper bandwidth of the two lowest lying branches similar (it looks on the first sight when ignoring the spin-polaron fit).
5) It is surprising that the lowest energy excitation branch is carrying less weight than higher branches. This is somewhat counterintuitive. Do the authors have an intuition why this is so? Is there more spectral weight pushed toward the low energy excitation with increasing D?
6) With increasing U/t, the polaron quasi-particle becomes increasingly extended in real space. What are the implications of using a 2x2 unit-cell for the excitation ansatz?
7) Even though the excitation ansatz is tailored toward extracting low energy properties, often high energy properties are reproduced very nicely too. I would be interested in seeing the spectra at an extended frequency regime (maybe this could be shown in an Appendix?).

Once these questions have been carefully addressed, I think this work is suitable for publication in SciPost, as it develops a novel computational method for computing spectra of 2D strongly correlated systems at zero temperatures.