Generating Function for Tensor Network Diagrammatic Summation

The understanding of complex quantum many-body systems has been vastly boosted by tensor network (TN) methods. Among others, excitation spectrum and long-range interacting systems can be studied using TNs, where one however confronts the intricate summation over an extensive number of tensor diagrams. Here, we introduce a set of generating functions, which encode the diagrammatic summations as leading order series expansion coefficients. Combined with automatic differentiation, the generating function allows us to solve the problem of TN diagrammatic summation. We illustrate this scheme by computing variational excited states and dynamical structure factor of a quantum spin chain, and further investigating entanglement properties of excited states. Extensions to infinite size systems and higher dimension are outlined.


I. INTRODUCTION
The study of quantum many-body systems using tensor networks (TNs) has witnessed great success in the last three decades [1][2][3]. Originally, TN methods were developed to efficiently capture ground state properties of many-body lattice models with short-range interaction [4][5][6][7]. Later on, numerous progress has been made in various directions, including determining low-energy excited states [8], exploring dynamical and finite temperature properties [9], and finding valuable applications in long-range interacting systems [10,11]. These developments not only deepen our theoretical understanding of many-body systems [12], but also bridge TN methods to real experiments [13].
New developments also bring challenges, however. Both quasiparticle excited states [5,[14][15][16] and global observables contain contributions with a sizable number of tensor diagrams, due to the fact that quasiparticle or local operators of global observables can be on arbitrary patch of the lattice. Except in a few cases where efficient summation techniques have been proposed [17][18][19][20], most notably the matrix product operator (MPO) representation of global observables [19], extensive and costly tensor diagram manipulation seems to be unavoidable and becomes the bottleneck in modern TN applications. Thus, an efficient and universal approach for TN diagrammatic summation is highly called for.
Another domain where diagrammatic summations frequently appear is the perturbation theory of interacting quantum fields [21,22]. There, for correlation functions containing summations of Feynman diagrams, one can introduce a source field and formally represent correlation functions as derivatives of the perturbed partition function, known as the generating functional method [21]. Given the close relation be- * Ji-Yao.Chen@mpq.mpg.de tween the trio of TN methods, many-body systems, and quantum field theory (QFT), and the pictorial similarity of tensor diagram and Feynman diagram, it is tempting to look for a generating function formalism in TNs, where certain derivatives can compactly represent the summations of TNs. This is plausible, as partition functions of classical statistical models are known to be representable as TNs [23].
In this work, inspired by the generating functional method in QFT, we propose a set of generating functions for TNs, which encode TN diagrammatic summations as leading order expansion coefficients. It then requires taking derivatives of the generating functions, which can be accomplished with automatic differentiation (AD) [24][25][26]. To illustrate the scheme, we investigate the low-lying spectrum of a quantum spin chain with periodic uniform matrix product state (MPS) and the excitation ansatz [14], and subsequently study entanglement properties of excited states which, to our knowledge, were rarely studied due to the overwhelmingly large number of tensor diagrams involved.

II. PERIODIC UNIFORM MPS AND EXCITATIONS
Let us consider a translationally invariant quantum spin chain with N sites. For simplicity, we assume the ground state is unique, and can be approximated by a periodic uniform MPS with the form: where the same rank-3 tensor A with dimension d × D × D is repeated on every site, and the left boundary is contracted with the right boundary. Here s i = 1, . . . , d represents basis of d-dimensional local Hilbert space, while D is the virtual bond dimension, controlling the accuracy of the MPS ansatz. We further denote the one-site translation operator as T withT |s 1 , s 2 , . . . , s N = |s N , s 1 , . . . , s N −1 , which sat-isfiesT N = 1. By construction, |Ψ(A) is translationally invariant with momentum k = 0, i.e.,T |Ψ(A) = |Ψ(A) . For a given model, the ground state tensor A can be optimized using the conjugate gradient method, with gradient obtained from AD of the computation graph for energy [26] [27].
With ground state tensor at hand, excited states can be constructed using the single mode approximation [28], which correspond to one-particle excitation and work well for a broad range of models [29][30][31]. In full generality, one can perturb the ground state by replacing one site tensor A with a new tensor B which is yet to be determined, and then build up a Bloch state using translation operator [5,14], taking the form: where tensor B contains variational parameters for excited state. |Φ k (B) then is an eigenstate of translation operator with eigenvalue e ik , where momentum k = 2πm/N, m = 0, 1, . . . , N − 1. Due to momentum superposition, a summation of N different tensor diagrams appears in Eq. (2), which will be our main focus. Since |Φ k (B) depends on tensor B linearly, variationally optimizing B boils down to a generalized eigenvalue problem: where E is the generalized eigenvalue, and H (N) is the effective Hamiltonian (norm) matrix in the variational space, with H µν = Here B is complex conjugate of B, whose component after vectorization is denoted as B ν . Since momentum is a good quantum number, we have suppressed the dependence of H, N, E, and B on momentum k. Solving the generalized eigenvalue equation in each momentum sector, one recovers the low-energy spectrum [32].
To construct H and N, one needs to sum over N different tensor diagrams for each, with MPO representation of the HamiltonianĤ. These extensive tensor diagram summations are the main obstacles of computing the excitation ansatz, rendering manipulating excited states unfavorable. Below we introduce our formalism based on simple yet powerful generating functions with the following strategy: to compute H or N, we will first construct a suitable generating function, and then use AD to compute the derivative [26], which will reproduce H or N and is much simpler than directly summing all diagrams. In this way, we will get rid of all the tensor diagram summations, making it possible to investigate detailed properties of excited states. Note that, unlike generating functionals in QFT, whose closed-form expressions are rare, the TN generating functions and their derivatives can be computed numerically exact with AD. We find that, depending on the origins of diagrammatic summation, the generating functions can be divided into two classes, one for TN state and the other for TN operators, which we introduce separately.

III. GENERATING FUNCTION FOR STATE
As shown in Eq. (2), the extensive tensor diagrams only differ by the location of tensor B and a position dependent phase factor. It is insightful to make the following observation: for a given tensor B, the excitation ansatz Eq. (2) can be expressed as: where the tensor on the j-th site in |G Φ (λ) is given by represented by blue squares. Here, to simplify the notation, we have suppressed the dependence of |G Φ (λ) on tensor A, B and momentum k, keeping only λ-dependence explicitly. Expanding |G Φ (λ) into power series of λ, we find that the ground (excited) state |Ψ(A) (|Φ k (B) ) is contained in the zeroth (first) order term, both of which lie in the tangent space of the MPS manifold, while higher order terms are outside of the tangent space due to nonlinearity in tensor B [33]. Thus, we can eliminate the tensor diagram summation in |Φ k (B) by computing the first order derivative of |G Φ (λ) . Interestingly, Eq. (3) bears a similarity with the generating functional in QFT, where the parameter λ plays the role of source field in the latter. Note that, although |Ψ(A) and |Φ k (B) are translationally invariant, the generating function |G Φ (λ) is in general not invariant under one-site translation, except at momentum k = 0.
With |G Φ (λ) , the norm square of excited state can be expressed as . Using translation invariance of |Φ k (B) , we can lower the order of derivative with the following generating function for the excited state norm: with which, the norm square can be obtained as Here, the local tensor on site-j of the ket layer is A j (λ), the same as appearing in Eq. (3). Before proceeding further, let us discuss how to use generating functions in practice. Taking Eq. (4) as an example, we first compute G ||Φ|| (λ = 0) by contracting the network in a conventional manner with computational complexity O(D 5 ), followed by a back-propagation using AD, with which the first order derivative at λ = 0 is obtained automatically, hence the norm square [26]. Since the computational complexity of AD grows with the order of derivative, it is advisable to utilize a generating function with which low order derivative suffices.
With Eq. (4), it is straightforward to find the generating function for the norm matrix N: where we simply omit tensor B in the bra layer of G ||Φ|| (λ). From Eq. (5), the norm matrix N can be obtained as N = . Unlike the scalar G ||Φ|| (λ) where the derivative is taken with respect to λ at λ = 0, G N (λ, B) is a vector and the derivative is taken with respect to B at λ = 1, B = 0, which can also be computed with AD [34].

IV. GENERATING FUNCTION FOR OPERATOR
Another type of TN diagrammatic summation originates from global observables. E.g., when computing structure factor, one needs to take into account correlations of local operator at all distances. Moreover, when applying TN methods to long-range interacting models which appear naturally in Rydberg quantum gas and quantum chemistry systems, efficient and faithful encoding of the long-range interaction is one of the main issues [35]. Although approaches based on MPO and projected entangled-pair operator (PEPO) have been proposed, the bond dimension can be quite large, making them less appealing [36][37][38][39]. Here we provide generating functions for global observables, without involving MPO whenever possible.
Let us first consider spin structure factor. Here, the summations arise from the Fourier transform of onsite spin operator A generating function can be introduced as follows: which is a tensor product of onsite operators, withŜ α j (λ) = I + λe −ikjŜα j , λ ∈ R, acting on the j-th site. One can then . Taking two layers of Eq. (6) with independent parameters λ, λ , followed by a second order derivative with respect to λ and λ , one arrives at the static structure factor. Higher order moments of local operator can be similarly obtained. Note that, with translation symmetry, one can combineĜ S α (λ) and local spin operator to reduce the order of derivative. Interestingly, by removing the λ-dependence in the first operator, and replacing the phase factor in the others with corresponding distance-dependent coefficient, Eq. (6) can be used to efficiently encode long-range power-law decaying interactions. Moreover, this approach can be easily generalized to higher dimension, where PEPO technique for long-range interactions is imperfect [37]. A simple counting suggests the number of terms for expectation value of long-range interaction is N 2 , while with generating function Eq. (6), it is only of order N , where each term differs by the position of the first operator. Exploiting translation symmetry, the latter case can be further reduced to a single term. Without translation symmetry, a combination of Eq. (6) and MPO technique is possible, which may lower the bond dimension of MPO. While Eq. (6) relies on the fact thatŜ α is an onsite operator, the scheme can be generalized to the case where the local operator has support on more than one site, e.g., local Hamiltonian operator. Similar to Trotter gates, we introduce the following generating function for Hamiltonian operator under periodic boundary conditionĤ = N j=1ĥ j,j+1 (assuming N even and nearest-neighbor interaction): where the operator acting on the (j, j + 1) bond is given bŷ h j,j+1 (λ) = I + λĥ j,j+1 , λ ∈ R, represented by a red rectangle. It is then straightforward to verifyĤ = d dλĜ H (λ) λ=0 . For the effective Hamiltonian in the generalized eigenvalue equation, one could take the generating function for norm matrix Eq. (5), and insert a MPO representation of the Hamiltonian, or the generating function Eq. (7) with an additional derivative. The latter could be advantageous if the MPO operator has a relatively large bond dimension [40]. Removing the boundary term,Ĝ H (λ) can also be applied for systems with open boundary. Furthermore, taking two layers ofĜ H (λ) one can efficiently compute the energy variance.
A close comparison between generating functions for TN state, operator, and partition function in QFT suggests that, for excited state, it is constructed by introducing source field to the ground state tensor, while for operators, identity operator is taken as the reference. This formalism also reminds us of Lie algebra, whose generators can be obtained by expanding group elements around the identity. In that case, the generating function has an exponential form. Indeed, generating functions for moments and cumulants have been introduced in Ref. [41], which took exponential forms with possible Trotter error and a finite difference method was then employed. In contrast, our approach does not have Trotter error, nor finite difference error, and the generating functions for the state and operator can be unified in a simple way.

A. Variational excited states and dynamical structure factor
We now present the numerical results using generating functions. The model we use for benchmark is the spin-1 Heisenberg chain with Hamiltonian:Ĥ = J N j=1 S j · S j+1 , where S j = (Ŝ x j ,Ŝ y j ,Ŝ z j ) are spin-1 operators on site-j and J = 1 is taken as the energy unit. We use periodic uniform MPS to approximate the ground state, and compute the low-energy spectrum with excitation ansatz, where generating functions introduced above are used. For small system size (N = 16), we compare the variational energy spectrum with that obtained from exact diagonalization (ED). As shown in Fig. 1(a), with bond dimension D = 24, the variational result agrees remarkably well with ED, with maximal relative deviation 6 × 10 −4 for high energy levels, and the level degeneracy is also recovered to a high precision. Moreover, although the low-energy excitations around k = 0 are two-magnon continuum [42,43], it is nevertheless reproduced well by oneparticle excitation ansatz. With excited states available, dynamical properties of the system can be investigated in a straightforward manner by constraining the total Hilbert space to subspace spanned by variational ground and excited states [44]. One of the key observables is the dynamical spin structure factor (DSF), which reveals properties of quasiparticles and is accessible in neutron scattering experiments. The DSF is defined as: | and E 0 (E k n ) being energy of ground state |Ψ(A) (n-th excited state |Φ k (B n ) with momentum k). The delta function is replaced by a normalized Gaussian with broadening width σ = 0.05, and we further takê S z as the operator in DSF. The DSF can be efficiently computed with generating functions Eqs. (3), (6), and the results for k = 3π/4, 7π/8, π are shown in Fig. 1(b)-(d) with system size N = 16. Comparing with ED, the line shapes are reproduced well by excitation ansatz.
For larger system size (N = 60), we quantify the quality of the variational result with the energy variance Var(E) = Ĥ 2 − Ĥ 2 , computed using generating functions (see Eqs. (4), (7)). Here the bond dimension is D = 30. As shown in Fig. 1(e), the energy variance of the variational ground state and one-magnon excited states remains small (on the order of 1 × 10 −3 ), considering the large system size, while the variational multi-magnon excitations appear to have larger energy variances. The Haldane gap can be read off as ∆ = 0.4105, in agreement with Refs. [29,43]. Further evaluating DSF using excited states (shown in Fig. 1(f)), we find a strong peak appearing at (k, ω) = (π, ∆), where the first excited state is located. This confirms that the elementary excitation is magnon with momentum k = π, in agreement with the variational spectrum where the first excited state is a triplet at k = π. Vanishing DSF around k = 0 is also consistent with excitations being two-magnon continuum in that region [43].

B. Rényi entropy of excited states
Apart from variational energy spectrum and DSF, the generating functions also allow us to further study properties of excited states in great detail. Here we use them to investigate entanglement properties of excited states, e.g., Rényi entropy, which has received considerable interest in recent study [45][46][47][48]. It is well known that for 1D gapped systems, the entanglement entropy of ground state saturates with increasing subsystem size [49]. However, much less is known for excited states. Traditionally, computing Rényi entropy for excited states requires multiple summations and thus is hard to achieve. Here we explore this question using generating functions without any summations.
For a normalized excited state |Φ with bipartition of the system L and L, reduced density matrix (RDM) of subsystem L (with size l) is given by ρ L = Tr L |Φ Φ|, and the Rényi entropy is then defined as S (n) = 1 1−n lnTr L ρ n L . Here we will focus on the n = 2 case, which can be computed by taking two copies of |G Φ (λ) and G Φ (λ)|, each with an independent parameter λ. Through AD of a single diagram, a fourth order derivative (one for each layer) at zero point gives directly the Rényi entropy.
The Rényi entropy with n = 2 for spin-1 Heisenberg model is shown in Fig. 2, using excited states in Fig. 1(e). The saturation of ground state Rényi entropy with increasing l is evident in the inset of Fig. 2, as expected. However, for excited states, we find that generically the Rényi entropy does not saturate with increasing l. To further quantify the effect of quasiparticles in excited states, we consider the ratio between excited and ground state traces: is RDM of excited (ground) state with subsystem size l. Theoretical studies have shown that, under the assumption of large momentum and energy gap, F (2) takes a universal form [47]: 2 , with x ≡ l/N . Indeed, the k = π/2 result agrees with the theoretical prediction, although our model is neither integrable nor free. However, clear deviation is also observed for other momenta, which could be ascribed to the small energy gap or momentum and deserves further study.

VI. DISCUSSION
Above we have shown that the generating functions can be used to compute one-particle excitation spectrum of a quantum spin chain with finite length. The results are encouraging, as the low-energy content is fully captured by MPS excitations, finding direct applications to investigate edge properties of two dimensional (2D) system [3]. The generating functions can also be adapted for tangent space based excitation ansatz without translation symmetry [50]. There are several directions for further extension. To go beyond one-particle case, one option would be modifying the ground state tensor A with A + i λ i B i , similar to introducing several source fields in partition function of QFT. For infinite size systems, one can combine generating functions with fixed-point methods, which will enable computing excitation spectrum of 2D system with projected entangled-pair states, including anyonic excitations in topological phases of matter [3].

VII. CONCLUSION AND OUTLOOK
In this work, we have introduced a set of generating functions for both TN states and operators, thereby eliminating extensive diagrammatic summations in modern applications of TN methods. Using generating functions, we have shown that excitation spectrum of quantum spin chain can be computed efficiently and accurately with periodic uniform MPS, and the procedure is formally the same as ground state search using TNs. Moreover, the generating functions allow us to investigate dynamical structure factor of the system and entanglement property of excited states in a convenient way, the later of which is beyond the capability of traditional methods. We envision the generating functions introduced here will be powerful in the next generation of tensor network algorithms and applications.
Note added -During the preparation of this manuscript, we became aware of an independent work by Boris Ponsioen and Philippe Corboz [51].