Generating function for projected entangled-pair states

Diagrammatic summation is a common bottleneck in modern applications of projected entangled-pair states, especially in computing low-energy excitations of a two-dimensional quantum many-body system. To solve this problem, here we extend the generating function approach for tensor network diagrammatic summation, a scheme previously proposed in the context of matrix product states. Taking the form of a one-particle excitation, we show that the excited state can be computed efficiently in the generating function formalism, which can further be used in evaluating the dynamical structure factor of the system. Our benchmark results for the spin-$1/2$ transverse-field Ising model and Heisenberg model on the square lattice provide a desirable accuracy, showing good agreement with known results. We then study the spin-$1/2$ $J_1$-$J_2$ model on the same lattice and investigate the dynamical properties of the putative gapless spin liquid phase. We conclude with a discussion on generalizations to multi-particle excitations.


I. INTRODUCTION
Strongly correlated quantum systems occupy a central position in condensed matter physics, often triggering exotic behavior at low temperatures.For low-dimensional systems, where quantum effects are pronounced, the entanglementbased tensor network method is now widely recognized as an ideal tool for both analytical and numerical studies of these systems [1,2].Beyond its immense success in exploring ground state properties due to its structure that adheres to the area law, tensor network methods also enable the study of low-energy properties above the ground state, including dynamical correlations and entanglement dynamics.The former are often directly measurable in spectroscopic experiments within condensed matter physics, such as inelastic neutron scattering experiments in quantum magnets [3], while outof-equilibrium properties are accessible in quantum simulators [4].However, in systems extending beyond one spatial dimension, efforts to fully understand low-energy excitations in quantum many-body systems through tensor networks are still in their early stages.Consequently, there is a high demand for an efficient and accurate method to compute lowlying excitations and related dynamical correlation functions in two-dimensional (2D) systems.
One physically motivated way of studying excitations is to use the tensor network generalization [46,47] of the Feynman-Bijl ansatz [48][49][50][51] or the single mode approximation [52,53], which is a variational ansatz for an excited state in the tangent space of the ground state tensor manifold [54,55].In this construction, the ground state is perturbed locally by introducing an "impurity" tensor to represent a local quasiparticle.By making a momentum superposition of such a local perturbation, one obtains a natural representation of a low-energy excited state.This approach has been successfully applied to build excitations on top of ground states of matrix product state (MPS) form in 1D [56][57][58][59][60][61] and quasi-1D [60,62,63] systems, where contractions arXiv:2307.08083v2[cond-mat.str-el]4 Mar 2024 can be performed exactly and efficiently.For a 2D quantum many-body system in the thermodynamic limit, this construction leads to the sum of infinitely many copies of PEPS with infinite size, so that computing expectation values or optimizing the variational energy is not straightforward.This can still be achieved through summing a geometric series of channel operators built from the boundary MPS approach [28,29] or alternatively through summing diagrams within the corner transfer matrix renormalization group method [30,31].Nevertheless, the number of tensor diagrams one needs to take care of in these approaches is substantial.Additionally, it is not clear how to apply these approaches to finite-size systems, since both approaches rely on a fixed-point iteration which is designed for an infinite system.
Instead of directly conducting the tensor diagram summation, in Ref. 64 some of the authors and collaborators have introduced a set of generating functions for tensor network diagrammatic summations, and applied this scheme to study the excitation spectrum in 1D systems with MPS.The key idea of the generating function is that the relevant tensor diagram summations can be expressed as low-order derivatives of a single diagram.The origin of this scheme can be traced to the fact that interactions in quantum many-body systems are local and the low-energy excited states only contain one or few quasi-particle excitations.Thus, inspired by the generating functional method in quantum field theory, one can introduce a source term in the tensor network and express the excited state as a first-order derivative of a new tensor diagram.With straightforward extensions, the effective Hamiltonian and norm matrices in the variational parameter space of excitations can be obtained, making it simple to evaluate expectation values or optimize the variational parameters.
In this work, we will extend this idea to PEPS in the thermodynamic limit.We begin with a short explanation of the infinite PEPS (iPEPS) ground state ansatz and how to construct low-energy excited states on top of a PEPS ground state in Sec.II.Then we introduce the generating function for PEPS and illustrate how to solve several practical issues in Sec.III.With impurity tensors obtained through solving an eigenvalue problem, it is straightforward to study the low-energy properties of 2D systems, e.g., the dynamical spin structure factor.In Sec.IV we will first benchmark our approach using the wellstudied transverse-field Ising model and Heisenberg model on the square lattice, and then apply the method to the more challenging spin-1/2 J 1 − J 2 model on the same lattice.We conclude with some perspectives on generalizing this approach to more general excitations and other settings in Sec.V.

A. PEPS as a variational ground state
For a 2D quantum many-body system, despite an exponentially large Hilbert space, the ground state of a gapped system obeys the so-called entanglement area law [65,66].This guarantees the validity of the direct extension of MPS for 1D quantum systems to PEPS as a suitable numerical tool for studying ground states of gapped 2D systems.For 2D critical systems, which could violate the entanglement area law, PEPS have also been shown to be a powerful tool which can capture the critical properties via a so-called finite entanglement scaling approach [67][68][69].
On the square lattice, a PEPS is defined by a set of rank-5 tensors, with one tensor for each site.Using translation symmetry with a suitable unit cell of tensors, PEPS can be defined directly in the thermodynamic limit, giving rise to the so-called iPEPS method.Taking an iPEPS with a one-site unit cell as an example, the wave function is given by the following tensor network form: where the black (red) bonds represent the virtual (physical) indices with bond dimension D (d).The accuracy of approximating the ground state of a 2D quantum system with the PEPS ansatz can be controlled by the virtual bond dimension D.
Although in general exact contraction of PEPS in the thermodynamic limit is not possible, various controlled approximate contraction schemes for PEPS have been developed, including the boundary MPS approach [70][71][72], the tensor renormalization group and its various generalizations [73][74][75], and the corner transfer matrix renormalization group (CTMRG) method [23,[76][77][78].In this work we will use the CTMRG method, where the basic idea is to approximate the surroundings of every site using a set of tensors with finite bond dimension, typically denoted as the environment bond dimension χ.As an illustration, we first trace out the physical index of A and Ā, where Ā represents the complex conjugate of tensor A, creating the norm of the PEPS wave function: where The thick black bonds in Eq. ( 3) thus have dimension D 2 .Through the CTMRG procedure, one can obtain a series of environment tensors that serve as the approximation of all surrounding tensors: Here, the C i tensors are the corner tensors and the T i 's represent the edge tensors (i = 1, 2, 3, 4).The green bonds (bonds that connect C and T tensors) have dimension χ and typically one would choose χ ≥ D 2 to ensure the accuracy.Note that, to obtain the environment tensors, one would need to introduce projectors to truncate the growing environment bond dimension in each CTMRG iteration.Apart from the ground state calculation, we will show that the environment tensors and the projectors will also play important roles in the calculation of generating functions for excited states.Using environment tensors, the reduced density matrix of a local region can be constructed and observables can be evaluated.E.g., given a local Hamiltonian of the form H = i h i , the energy density can be calculated as: Note that, when computing Eq. ( 5), the approximations due to projectors in CTMRG are the same for the numerator and denominator.This point will be illuminating when considering the contraction scheme for generating functions of PEPS.The remaining task is to optimize the ground state tensor A to minimize the variational energy.Conventionally, this can be achieved through imaginary time evolution, where the local Hamiltonian is used to construct Suzuki-Trotter gates, and after a sufficiently long imaginary time evolution a wellapproximated ground state ansatz can be obtained [70,79,80].However, the initial ansatz for imaginary time evolution tends to significantly affect the final outcome, hindering the simulation without any tentative knowledge of the target model and sometimes causing a serious issue of getting trapped in the local minima of Hilbert space.Additionally, due to the approximations involved in truncating the PEPS after applying Suzuki-Trotter gates, the imaginary time evolution does not always converge to the variationally optimal PEPS tensor.
Instead of doing imaginary time evolution, minimizing energy through the energy gradient, i.e., variational optimization, has been proven to be a more accurate optimization scheme [81,82].When the number of variational parameters is small (for example, when the PEPS tensors are strongly constrained by symmetries), variational optimization can be done using a simple finite difference approach to compute the energy gradient [19,20,83,84].For generic tensors, a direct evaluation of the gradient requires a summation of a large number of tensor diagrams [72,81,82], but the numerical derivatives can be obtained through automatic differentiation (AD) [85].By storing the computational graph of PEPS [14,15] in the memory during the forward-mode calculation, the numerically exact gradients can be evaluated through the back propagation [86].We then make use of those gradients to further optimize our tensors in order to lower the variational energy.
In the following simulations, we have checked that for each bond dimension D, the chosen value of χ is either large enough and the results barely change when further increasing χ, or the largest χ within computational reach.In search of the ground state, we have imposed C 4v symmetry of the square lattice on the local tensor, since the models we studied do not break this symmetry.In the actual numerical computation, the ground state optimization is carried out using the peps-torch package [87], where the energy gradient is obtained using the back propagation of AD provided in PyTorch.For the excited state computation (see later sections), we have written our codes in Python and used PyTorch for back propagation.To lower the bar for PEPS practitioners, we have made a sample code publicly available [88].

B. Quasiparticle excitation with PEPS
The excited state for a given many-body Hamiltonian can be thought of as a quasiparticle excitation on top of the ground state [89].As mentioned in the introduction, in real space such an excitation can be approximated through the tensor network version of the single mode approximation, where we replace one ground state tensor with some "impurity" tensor at a certain location [90].We then sum up all different tensor network configurations generated by the translation operator, taking into account the corresponding phase factors to construct an eigenstate of the translation operator with a given momentum.
With the ground state expressed as a one-site translation invariant iPEPS, the one-particle excitation with PEPS can be written as: where Tx ( Ty ) represents the translation operator in the x (y)direction with its eigenvalue being e ikx (e iky ), and B is the impurity tensor to be determined.
Due to the linear dependence of excited states on the impurity tensor, the variational optimization boils down to a generalized eigenvalue problem: where H and N are the effective Hamiltonian and norm matrix in the variational space.With translation invariance, N can be obtained by hollowing out the Ā tensor in the center of the bra layer, and the summation goes over every tensor graph with an empty site in the ket layer.Using channel operators with boundary MPS [28], this sum can be carried out at the same computational cost as ground state computation.A slightly different summation scheme using CTMRG has also been proposed in Ref. [30].Similarly, the effective Hamiltonian H can also be computed, where an additional summation for the Hamiltonian operator appears.These summations are the main bottleneck of studying excitations using the quasiparticle ansatz.
Because of the gauge freedom in the PEPS representation (Eq.( 1)) in terms of the local tensor A, the linear subspace of excitations contains a few zero modes [28].These are reflected as zero eigenvalues in the norm matrix and should be projected out to make the above eigenvalue equation wellconditioned.We will come back to the conditioning of N below.

A. Algorithm
The issue of summing up many diagrams also appears in the MPS study of the one-dimensional quantum system, where some of the authors and collaborators have proposed generating functions to tackle this problem [64].Now we adopt this method in the current setting of PEPS in the thermodynamic limit.The basic idea of generating function is that the sum of extensive tensor diagrams can be expressed as a low-order derivative of a new tensor diagram.This is possible for excited states due to the fact that tensor diagrams only differ by the location of the impurity tensor.
For iPEPS, following Ref.64, we define the generating function for the one-particle excited state as: where G B stands for the generating function of the impurity tensor and takes the following position-dependent form: The one-particle excited state in Eq. ( 6) can now be evaluated by taking the derivative: One question immediately follows: namely, what would be the right environment tensors for ⟨G Φ (λ)|G Φ (λ)⟩ when computing physical observables?Note that, for PEPS in the thermodynamic limit, the contraction scheme relies on the translation invariance to interpret the environment tensors as fixed points of the one-dimensional transfer operator.For the generating function in Eq. ( 8), this translation invariance is lost for any nonzero λ.The solution to this issue comes from the following fact: after taking derivatives at λ = 0, all contributions of the sum can be viewed as correlation functions of local operators.Recall that when computing observables or correlation functions in the infinite MPS, one uses the fixed point of the transfer matrix as the boundary, suggesting that fixedpoint tensors of the transfer matrix of the infinite MPS ground state are the right environment tensors for generating function of infinite MPS.The same is true for the PEPS ground state in the thermodynamic limit.Indeed, we have tested that using a fixed point of the MPS transfer matrix in the generating functions gives the same result as directly summing all tensor diagrams.Thus, similarly for PEPS, we can use the environment tensors from CTMRG of the ground state PEPS (i.e., the state |G Φ (λ)⟩ at λ = 0) as the environment tensors for the generating function Eq. ( 8).
With the generating function for the excited state, the norm matrix and effective Hamiltonian can be similarly expressed as derivatives of a single network.Using translation symmetry, we can lower the order of derivatives by introducing two slightly new networks [64].For that, we first construct two (idealized) generating functions, G N (λ, B) and G H (λ, µ, B) for norm matrix and effective Hamiltonian, respectively.G N takes the following form: The green tensors in Eq. ( 11) are given by: where the tensor with a cross mark lies in the center of the lattice.The norm matrix can then be evaluated by taking the derivative: For the generating function of the effective Hamiltonian, we need to insert the local Hamiltonian between the bra and the ket layers.Note that, although the local Hamiltonian can be represented as a projected entangled-pair operator, it is rarely used in practice, due to the increased computational cost.Here we use the generating function for a local Hamiltonian, as introduced in Ref. [64], to represent the full Hamiltonian as a derivative of a new network.Then we sandwich the generating function of the local Hamiltonian term between the bra and ket layers.To illustrate, assuming that we only have the nearest neighbor terms, a new tensor can be defined as: where I is the identity matrix and î = x, ŷ.Then the (idealized) generating function for the effective Hamiltonian is given by the following: From Eq. ( 15), the effective Hamiltonian H can be obtained by taking a second order derivative Denoting the linear size of the bulk, i.e., the region encircled by the dashed box in Eqs.(11) and (15), as L x and L y , it appears that L x and L y have to be chosen carefully.In the case of infinite MPS, we have checked that the linear size in the bulk has to be compatible with the momenta, i.e., k = 2πm/L x (m = 0, 1, 2, • • • , L x − 1), to ensure exact zero modes in the norm matrix.We expect that the same requirement is also true in the two-dimensional case.Moreover, L x and L y should be large enough (compared to the maximal correlation length of the ground state) to mimic an infinite system, so that the error due to finite linear size becomes negligible.At the same time, with larger linear sizes, more k points can be considered.Note that, here we restrict to a finite bulk size, instead of an infinite one, due to the fact that for a given ground state tensor A and impurity tensor B, the generating function and its derivatives may not converge within the same number of iterations.This is because the forward calculation is computing the zeroth order of the expansion of generating functions at λ = 0, which typically resembles the tensor diagram for a local observable (see Eq. ( 11) and ( 13) for example), while the derivatives contain multi-point correlators.Therefore, it is conservative but safe to work with a finite bulk size.
With the norm matrix and effective Hamiltonian obtained from derivatives of generating functions, we can now solve the generalized eigenvalue equation in Eq. ( 7), and obtain the impurity tensor B for excited states.As mentioned before, due to gauge degree of freedom, exact zero modes appear in the norm matrix and one needs to project out the corresponding subspace, leading to the following generalized eigenvalue equation: where the subspace projector P projects out the exact zero modes of the norm matrix.Note that, up to this point, there is no approximation involved and the scheme can be directly applied to the infinite MPS context where computation can be made essentially exact.However, there is one more issue for PEPS to address: for large tensor network graphs in 2D (Eqs.( 11) and ( 15)), its computational complexity grows exponentially with system size.Therefore, the exact contraction is not feasible.By utilizing the projectors during the CTMRG process, on the other hand, we can approximate tensors on each column as the following: where the tensor projectors (dark green triangles) can be evaluated through the singular value decomposition [23] and a similar process can be conducted for the row and corner tensors.With the projectors inserted, now we show the complete tensor network for the norm matrix in Fig. 1.Note that here we show explicitly how to compress our tensor graphs vertically but in fact the similar action is also taken along the horizontal direction.We can then take the average of these two compressed tensor graphs as G N (λ, B).
The motivation for the above mentioned splitting the contraction of G N (λ, B) into two is due to the insertion of local Hamiltonian terms for G H (λ, µ, B).Note that the tensor projectors intuitively would reduce the correlation between nearby sites, which is undesirable when inserting the local Hamiltonian operators.Therefore, in practice, we conduct the insertion of Ĥr,r+x and Ĥr,r+ŷ separately, resulting in the tensor contractions along the horizontal and vertical directions, < l a t e x i t s h a 1 _ b a s e 6 4 = " K s Z r i n u r S 2 S q S b M k 9 q e T i O 9 R q j 9 y f z U K T m z S p + E s b I l D Z m r v y c m N N J 6 H A W 2 M 6 J m q J e 9 m f i f 1 0 l N e O N P u E x S g 5 I t F o W p I C Y m s 7 9 J n y t k R o w t o U x x e y t h Q 6 o o M z a d g g 3 B W 3 5 5 l T Q r Z e + q f P l w U a r e Z n H k 4 Q R O 4 R w 8 u I Y q 3 E M N G s B g A M / w C m + O c F 6 c d + d j 0 Z p z s p l j + A P n 8 w f Z j Y 2 H < / l a t e x i t > T 2 < l a t e x i t s h a 1 _ b a s e 6 4 = " 8 c e y t C t A n i r U L Y M g q k W L 9 W 4 S y N k = " > A A A B 6 n i c b V D L S g N B E O y N r x h f U Y 9 e B o P g K e x K f B y D X j x G z A u S J c x O e p M h s 7 P L z K w Q Q j 7 B i w d F v P p F 3 v w b J 8 k e N L G g o a j q p r s r S A T X x n W / n d z a + s b m V n 6 7 s L O 7 t 3 9 Q P D x q 6 j h V D B s s F r F q B 1 S j 4 < l a t e x i t s h a 1 _ b a s e 6 4 = " 8 c e y t C t A n i r U L Y M g q k W L 9 W 4 S y N k = " > A A A B 6 n i c b V D L S g N B E O y N r x h f U Y 9 e B o P g K e x K f B y D X j x G z A u S J c x O e p M h s 7 P L z K w Q Q j 7 B i w d F v P p F 3 v w b J 8 k e N L G g o a j q p r s r S A T X x n W / n d z a + s b m V n 6 7 s L O 7 t 3 9 Q P D x q 6 j h V D B s s F r F q B 1 S j 4 R q j 9 y f z U K T m z S p + E s b I l D Z m r v y c m N N J 6 H A W 2 M 6 J m q J e 9 m f i f 1 0 l N e O N P u E x S g 5 I t F o W p I C Y m s 7 9 J n y t k R o w t o U x x e y t h Q 6 o o M z a d g g 3 B W 3 5 5 l T Q r Z e + q f P l w U a r e Z n H k 4 Q R O 4 R w 8 u I Y q 3 E M N G s B g A M / w C m + O c F 6 c d + d j 0 Z p z s p l j + A P n 8 w f Z j Y 2 H < / l a t e x i t >  < l a t e x i t s h a 1 _ b a s e 6 4 = " F / 0 5 x q / 0 P M S m C X I f W q K N s j J R q m 4 = " > A A A B 6 n i c b V D L S g N B E O y N r x h f U Y 9 e B o P g K e z 6 P g a 9 e I y Y F y R L m < l a t e x i t s h a 1 _ b a s e 6 4 = " F / 0 5 x q / 0 P M S m C X I f W q K N s j J R q m 4 = " > A A A B 6 n i c b V D L S g N B E O y N r x h f U Y 9 e B o P g K e z 6 P g a 9 e I y Y F y R L m < l a t e x i t s h a 1 _ b a s e 6 4 = " F / 0 5 x q / 0 P M S m C X I f W q K N s j J R q m 4 = " > A A A B 6 n i c b V D L S g N B E O y N r x h f U Y 9 e B o P g K e z 6 P g a 9 e I y Y F y R L m < l a t e x i t s h a 1 _ b a s e 6 4 = " F / 0 5 x q / 0 P M S m C X I f W q K N s j J R q m 4 = " > A A A B 6 n i c b V D L S g N B E O y N r x h f U Y 9 e B o P g K e z 6 P g a 9 e I y Y F y R L m respectively.More specifically, we include G Ĥ in all the vertical (horizontal) bonds and compress the tensor graph for each column (row).
The scheme is now clear and all we need to do is to first contract the finite tensor network in Eqs. ( 11) and ( 15) (with tensor projectors inserted), and then compute the derivatives using either AD or any other ways of numerical finite difference approach.Then, by solving the eigenvalue equation Eq. ( 17) after getting N and H, we obtain the impurity tensors with corresponding energies.With zero modes in the norm matrix projected out, the excited states constructed with the impurity tensors are then orthogonal to each other, and also orthogonal to the ground state.
Note that, due to the projectors inserted in Fig. 1, the zero modes in the norm matrix are no longer exact, and instead become fuzzy.Therefore, we use the subspace projector P to truncate the basis by discarding the norm matrix eigenvector with eigenvalue close to zero [30].Denoting the eigenvalue decomposition of N as N = vΛv † , we keep the desired subspace of v and construct ṽ, which serves as P in Eq. (17).One subtle point is that the selection of subspace cannot be simply determined a priori.Practically, we start from the leading eigenvalue of norm matrix and gradually enlarge the size of the subspace by including next leading ones.If adding certain norm matrix eigenvector changes the energy eigenvalues drastically, we would discard this vector when constructing the subspace.So far, we cannot conclude a strict protocol in choosing the subspace projectors, and occasionally need to examine different threshold for the subspace to obtain stable energy eigenvalues.Note also that, a suitable gauge in the ground state tensor may be helpful for improving the conditioning of the norm matrix [91].With these measures, we solve Eq. ( 17), and the excited state can then be constructed with B ν for further investigation of its properties.

B. Discussion of the algorithm
Comparing to the previous generating functions for the finite-size MPS [64], two new ingredients are included in the infinite PEPS case.The first one is concerning the infinite size in typical PEPS calculations.Although our scheme can be easily generalized to the finite PEPS, due to the large computational cost, iPEPS is typically used in the literature.For that purpose, we have introduced ground state environment tensors in computing the norm matrix and effective Hamiltonian of the variational space.
The second ingredient is about the contraction of the generating functions, where we have used ground state projectors to achieve an efficient contraction.The idea of using projectors from the ground state CTMRG procedure can further be put on firm grounds as follows.That is, all the correlations in our ansatz are essentially mediated by the ground state tensors, and all contributions to the effective Hamiltonian and norm matrix can be viewed as correlation functions of local operators in the ground state.Unraveling the CTMRG procedure, one would find that the CTMRG procedure is essentially inserting tensor projectors (which was computed selfconsistently) into the original network.Moreover, when conducting the forward calculation of generating functions with AD, we assign the impurity tensor being a tensor with all elements equal to zero (see Eqs. ( 13) and ( 16) where the derivatives are computed at B = 0).Therefore, the contraction is akin to the normal ground state CTMRG process.Thus it is natural to use the same projectors for the generating functions, which would not lower the accuracy.However, if one computes the generating function for excited state at nonzero λ or nonzero B tensor, which would indeed be the case when part of the derivatives is obtained with finite difference approaches, the projectors may need to be recomputed, taking into account the effect of B tensor.In that case, the computation time and memory consumption would also be higher.Nevertheless, here we refrain from making a detailed comparison between various possibilities, but leave them to future developments.We also note that this scheme is not restricted to CTMRG and in fact can be generalized to other PEPS contraction methods, e.g., tensor renormalization group [75] and boundary MPS approach [70].
In comparison with previous approaches using PEPS for constructing one-particle excited states, our approach largely reduces the number of tensor diagrams to be considered.In Ref. 28 the tensor contraction is done by evaluating the cornershaped transfer matrix, with norm matrix and effective Hamiltonian obtained by summing up the tensor diagrams whose impurity tensors in the bra and ket layer are located in different relative positions.While the corner transfer matrices are also adopted in Ref. 30 for the environment tensors, again due to the position change many tensor graphs need to be taken into account.
The approach proposed in Ref. 31 comes closer to the approach presented here, which relies on AD to do the tensor diagram summations.Let us now compare with that one here.The most prominent feature of our approach is that we only need to deal with one tensor diagram for N or H, where all the tensors are fixed before contraction.In Ref. 31, one needs to take care of the additional environment tensors (coming from impurity tensors), which were computed iteratively.In this sense, the approach presented here gets rid of all the intricate summations and updating steps in PEPS excitations.This simplification will play an essential role when considering multiparticle excitations with PEPS.
Before showing the applications, let us also mention that the generating functions are in fact independent of AD.In some cases, the derivative is taken only with respect to one parameter λ and therefore can be easily carried out using the finite difference method.Indeed, using finite difference may become crucial when the bond dimension is large, since in that case AD might be limited due to memory issues.Nevertheless, in this work we use AD for two reasons.Firstly, in obtaining the norm matrix and effective Hamiltonian, multiple parameters are involved and thus the back propagation mode of AD is more efficient.Secondly, being numerically exact, AD can avoid finite difference error.While the full expression of norm matrix and effective Hamiltonian may not be necessary for Lanczos type diagonalization algorithms (and thus only one parameter is involved in the derivative), we will not explore this possibility here and leave it to future considerations.

IV. APPLICATIONS
Now we benchmark our algorithm with two well-studied models, namely the transverse-field Ising model and spin-1/2 Heisenberg model on the square lattice, and then move to the more challenging J 1 − J 2 model.

A. Transverse-field Ising model
The spin-1/2 quantum Ising model with a transverse field in 2D is described by: where σ x,z are the Pauli matrices and the summation of ⟨i, j⟩ runs over all nearest-neighbor pairs.As the most widely known model for benchmarking, its phase transition can be accurately captured by PEPS [70,92].A high precision estimate of the transition point was made by the cluster Monte Carlo method with the value being h c = 3.04438 [93].When h → ∞, the ground state preserves the global Z 2 symmetry Û = i σ x i with a non-zero energy gap, while at h = 0 the ground state is two-fold degenerate with all spins aligning up or down, breaking the Z 2 symmetry.Close to the transition point, the energy gap decreases and becomes zero at h c , which is a Lorentz-invariant critical point.
In order to compare with the results in Refs.29 and 30, we first compute the lowest-energy excited state at each k and plot its energy versus momentum in Fig. 2(a).Here and below we have chosen the high symmetry path M → X → S → Γ → M → S in the first Brilliouin zone of the square lattice, shown in the right inset of Fig. 2(a).
Let us first examine the choice of bulk size N s = L x × L y in computing excited states.As shown in the left inset of Fig. 2(a), with h closer to the critical point and larger bond dimension D, the maximal correlation length (ξ) of the ground state (measured from transfer matrix spectrum of ground state PEPS) becomes larger, and one would need larger bulk size N s to converge the energy gap ∆ (at momenta Γ(0, 0)), in agreement with previous discussion on bulk size.Nevertheless, due to the relatively small bond dimension (and therefore small ξ) used in this work, the excited states shown in this work are relatively easy to converge.
From Fig. 2(a), one can further find that away from the critical point (e.g., at h = 2.5), the finite D effect is rather small, and the excitation energy does not show significant changes from D = 2 to D = 3. Close to the critical point (e.g., at h = 3.0), the finite D effect is most prominent in the lowest excited state, where the energy gap ∆ decreases with larger bond dimension.Our results show both qualitative and quantitative agreement with Ref. 29 and 30, and also agree with the series expansion result [94].
Besides the quasiparticle dispersion relation, we can also fit the numerically computed energy gap ∆ versus magnetic field h, taking a function form ∆ ∝ |h − h c | ν .In Fig. 2(b), we compare the numerical data with the analytically known critical exponent ν = 0.629971, and indeed find a good agreement, further confirming the accuracy of the generating function method.

B. Heisenberg model
Next, we examine another model widely used for the purpose of benchmarking, the antiferromagnetic (AFM) Heisenberg model: where we take J = 1 as the energy unit, and the sum runs over all nearest-neighbor pairs.While the AFM nature demands at least a two-site unit cell for the ground state ansatz, by rotating the local spin on one of the two sublattices, we can preserve the one-site translation symmetry for PEPS.
Before showing the excited states, as a good benchmark we compute the static structure factor, which can be measured by elastic neutron scattering experiments.The static structure factor is defined by the following equation [95]: where ⟨•⟩ c denotes the connected part of the correlator and α = x, y, z.Because of the AFM nature, the static structure factor possesses a peak at k = (π, π), and becomes zero at k = (0, 0).Given a variational PEPS ground state, the static structure factor can be evaluated efficiently using generating functions.Here we show S(k) with D = 4 in Fig. 3. Our result is in good agreement with the previous one using iPEPS [81], as well as those by other methods [95,96].More importantly, the results of different N s all look similar, fortifying the fact that for our method the choice of bulk size casts little effect on the final outcomes, as long as the bulk linear size is large comparing to the ground state maximal correlation length ξ.
In Fig. 4 we show our results for the lowest-energy dispersion of the Heisenberg model.One clear feature that we can see is the vanishing of excitation energy at Γ(0, 0) and X(π, π) with increasing D, which is desirable since the Heisenberg model is known to be gapless at the these two points.To further confirm the gapless nature, in the inset we provide the scaling of energy gap ∆ (at Γ(0, 0)) with 1/ξ.We can see that the energy gap indeed extrapolates to a value close to zero when ξ → ∞.Except for the gapped feature due to the finite bond dimension, our dispersion is close to the ones by earlier iPEPS results [29,30] and Gutzwiller-projected trial wavefunctions [97].Also, the well-known feature that the gap at (π, 0) is smaller than the one at (π/2, π/2) [97] is recovered in our simulation.In fact, the nature around M(π, 0) point is of interest and has been under investigation [98][99][100].In Ref. 30 it has been demonstrated using PEPS that the repulsion from the multi-magnon branches pushes the magnon at M(π, 0) to a lower energy [99].This feature is also captured in our simulation by showing a dip in the single magnon branch near the M point, which cannot be captured by the spin wave theory (dashed line).In the next section we will show that af- In the previous subsections, we have studied models with only nearest neighbor couplings using our method.Both the transverse-field Ising model and Heisenberg model have served as good benchmarks and have been studied earlier using conventional iPEPS construction for excited states [28][29][30][31].Here, we aim to push further and study the J 1 −J 2 model where the NNN coupling introduces frustration and thus exotic ground states [101,102].The Hamiltonian is given by: where the J 1 (J 2 ) term runs over all (next-)nearest-neighbor pairs on the square lattice.The J 1 − J 2 model has been under intense investigation.Since the frustration induces a sign problem for quantum Monte Carlo method, alternative numerical methods have been applied, suggesting that between J 2 ≈ 0.5J 1 and J 2 ≈ 0.6J 1 a ground state without magnetic order is realized (hereafter we fix J 1 = 1).Using PEPS methods, this vanishing of the staggered magnetization has been confirmed [14] using advanced scaling techniques -for any finite bond dimension, it appears that a variational PEPS ground state yields a nonzero magnetization.Despite an abundant amount of related studies for the ground state, very few results have been reported for the excited state and there is no result so far using the two-dimensional tensor network ansatz [103].We will now investigate the excitation spectrum for this frustrated model.
In order to include the NNN coupling, î in Eq. ( 14) now represents x, ŷ, x + ŷ, and x − ŷ.As a result, the tensor contraction that one needs to apply when preparing G H becomes slightly cumbersome.For example, the tensor graphs that we need to consider now involve two columns at the same time: which consumes more memory when storing the computational graph for the back propagation.Therefore, our calculation in this work is restricted to D = 4 for the J 1 − J 2 model and we will leave further consideration on increasing the available bond dimension, e.g., using symmetries in tensor networks to reduce computational cost, to future work.We note in passing that, the NNN coupling may also be handled using diagonal channel environments in Ref. 81, which however differs from the CTMRG environment used in this work.
In Fig. 5 we show the lowest-energy quasiparticle dispersion for the J 1 −J 2 model with J 2 = 0.3 and J 2 = 0.5 along a high symmetry path in the Brillouin zone.Similar to the case for the Heisenberg model, with increasing bond dimension, the energy gap becomes smaller, in agreement with the fact that the phase is gapless [104].Moreover, these dispersions are also close to those by the variational Monte Carlo (VMC) method [105].Besides the shape of the dispersion, we can clearly see that the energy gap becomes smaller at M(π, 0) when J 2 increases.More importantly, at J 2 = 0.5 our results using D = 2 indicate that even the M(π, 0) point may become gapless.This softening of the dispersion at the M point is suggestive of the formation of a spin liquid phase at intermediate values, where the mode at the M point would correspond to a two-spinon state [105].
Based on the variational wavefunctions for the excited states, it is straightforward to compute their contributions to the dynamical structure factor (DSF).The DSF is defined as ⟩|. E 0 and E k n are the ground-state and excited-state energies for |Ψ(A)⟩ and |Φ k (B n )⟩, respectively.And S α k is the Fourier transform of local spin operator with α = (x, y, z).Note that, in this case the forward mode AD could be computationally cheaper than reverse mode AD, since only a few parameters are involved in the derivative.In Fig. 6, we show the DSF for J 2 = 0.1, 0.3, 0.5, respectively, where all data are computed with PEPS bond dimension D = 3.The delta function in DSF has been replaced by a normalized Gaussian with broadening width σ = 0.1.Two clear features can be seen from our data: (1) the energy gap gets softened at the M(π, 0) point, and (2) the spectral weights get closer to the magnon branch (lowest excited energies) in the frustrated region.We have noted that in the previous results by VMC the authors also observed a similar trend in their DSF results, where a clear energy continuum is formed after J 2 /J 1 ≈ 0.45.In Fig. 6 we can clearly see that as J 2 goes to 0.5, various dispersions gather toward the lowest magnon branch, featuring the potential formation of the energy continuum.

V. CONCLUSION
In this work, we have extended the generating functions introduced in Ref. 64 to PEPS.After obtaining the well approximated ground state through variational optimization, we solve the eigenvalue problem in Eq. ( 17) where effective Hamiltonian and norm matrices are evaluated by taking the derivatives of their corresponding generating functions.We then utilize the eigenstates as the impurity tensors to construct the oneparticle excited state and employ the transverse-field Ising model and the Heisenberg model as benchmarks.Our results are consistent with previous ones [29,30] using PEPS that requires the serial summation on different parts of tensor graphs.We extend our consideration to the model with the NNN coupling and study the J 1 − J 2 model.Despite the increasing computational cost which hinders the computation with larger D, our results already show good agreement with previous ones by VMC [105].
Before closing, let us discuss how to generalize this idea to multi-particle excitations, especially those associated with spinons.It is known that spin liquids are characterized by fractionalized excitations, where the quasiparticles cannot appear alone but come in pairs, which are intrinsically twoparticle excitations.In the previous study, spinon excitation (among other fractional excitation) has been treated in MPS calculation using good quantum numbers, which appears as a local excitation with a non-local string attached [106,107].In the case of PEPS, fractionalized excitations can also be defined and their static properties, e.g., correlation functions and correlation lengths have been studied in Ref. 20.In that case, the local tensor of PEPS needs a certain gauge symmetry.Once constructed, the excited states containing spinon excitations can be obtained.To further study dispersion relation of spinon excitation, one may put one spinon infinitely far away or use a suitable boundary condition to reduce the problem into an effective one-particle problem.Generating functions can then be used for efficiently optimizing the oneparticle problem as we have shown above.
In summary, our proposal of combining the generating function with tensor network ansatz helps reduce the complexity in obtaining the physical properties or observables in-volving the summation of tensor networks, especially in computing quasiparticle excitations.The present method provides not only the dispersion relation of the quasiparticles but also a systematic way of knowing what the relevant quasiparticle for a given system is even if it is not clear a priori.It may then be possible to discover completely new quasiparticles in 2D quantum many-body systems.From another perspective, with the generating function method, the PEPS representation for excited state now takes almost the same form as ground state PEPS.Thus it may also be helpful for preparing excited state in quantum computers, which we leave to future study.
Note added -When finalizing this work, we became aware of a related work using AD methods for diagram summations in PEPS [108].

T 4 <
l a t e x i t s h a 1 _ b a s e 6 4 = " S v p L O j h 4 d R z 8 c H q l + R Z 2 6 r C C x e w = " > A A A B 6 n i c b V D L T g J B E O z F F + I L 9 e h l I j H x R H Z 9 H 4 l c P G K U R w I b M j s 0 M G F 2 d j M z a 0 I 2 f I I X D x r j 1 S / y 5 t 8 4 w B 4 U r K S T S l V 3 u r u C W H B t X P f b y a 2 s r q 1 v 5 D c L W 9 s 7 u 3

3 <
6 s B g A M / w C m + O c F 6 c d + d j 3 p p z s p l D + A P n 8 w f B K 4 1 3 < / l a t e x i t > C l a t e x i t s h a 1 _ b a s e 6 4 = " A 3 J 1 7 1 3 2 3 P 9 9 3 1 z Q A 1 H 6 C N o I d A E = " > A A A B 6 n i c b V D L T g J B E O z F F + I L 9 e h l I j H x R H Y N P o 5 E L h 4 x y i O B D Z k d e m H C 7 O x m Z t a E E D 7 B i w e N 8 e o X e f N v H G A P C l b S S a W q O 9 1 d Q S K 4 N q 7 7 7 e T W 1 j c 2 t / L b h Z 3 d v f 2 D 4 u F R U 8 e p Y t h g s Y h V O 6 A a B Z f Y M N w I b C c K a R Q I b A W j 2 s x v P a H S P J a P Z p y g H 9 G B 5 C F n 1 F j p o d a r 9 I o l

T 2 <
l a t e x i t s h a 1 _ b a s e 6 4 = " R k k l w f q l U s H P 3 F T w U n / D w E R U / B 0 = " > A A A B 6 n i c b V D L S g N B E O y N r x h f U Y 9 e B o P g K e w G X 8 e g F 4 8 R 8 4 J k C b O T 3 m T I 7 O w y M y u E k E / w 4 k E R r 3 6 R N / / G S b I H T S x o K K q 6 6 e 4 K E s G 1 c d 1 v J 7 e 2 v r G 5 l d 8 u 7 O z u 7 R 8 U D 4 + a q j 9 y f z U K T m z S p + E s b I l D Z m r v y c m N N J 6 H A W 2 M 6 J m q J e 9 m f i f 1 0 l N e O N P u E x S g 5 I t F o W p I C Y m s 7 9 J n y t k R o w t o U x x e y t h Q 6 o o M z a d g g 3 B W 3 5 5 l T Q r Z e + q f P l w U a r e Z n H k 4 Q R O 4 R w 8 u I Y q 3 E M N G s B g A M / w C m + O c F 6 c d + d j 0 Z p z s p l j + A P n 8 w f Z j Y 2 H < / l a t e x i t > T 2

O s 7 F 2 V 3 FIG. 1 .
FIG. 1. Complete tensor network for generating function of norm matrix with bulk size Ns = 4 × 4. The projectors inserted are taken from the CTMRG of the ground state.Similarly for the generating function of the effective Hamiltonian.

FIG. 2 .
FIG. 2. (a) The energy dispersion of the lowest-lying excitation for transverse-field Ising model at distinct k points for h = 2.5 and h = 3.0.We adopt Ns = 24 × 24 for D = 2 and Ns = 16 × 16 for D = 3 and confirm that further enlarging Ns does not change the results.The red triangles show the results from series expansion (SE).The left inset shows the energy gap with different Ns for (h, D, ξ) = (2.5, 2, 0.603), (h, D, ξ) = (3.0,2, 1.194), and (h, D, ξ) = (3.0,3, 2.531), where ξ is the maximal correlation length.The right inset shows the high symmetry path in the Brillouin zone.(b) The gap function along with magnetic field h.We adopt D = 4, Ns = 12 × 12, and the solid curve takes the scaling relation of 3D Ising universality class with the critical exponent ν = 0.629971.