Transfer Matrices and Excitations with Matrix Product States

We investigate the relation between static correlation functions in the ground state of local quantum many-body Hamiltonians and the dispersion relations of the corresponding low energy excitations using the formalism of tensor network states. In particular, we show that the Matrix Product State Transfer Matrix (MPS-TM) - a central object in the computation of static correlation functions - provides important information about the location and magnitude of the minima of the low energy dispersion relation(s) and present supporting numerical data for one-dimensional lattice and continuum models as well as two-dimensional lattice models on a cylinder. We elaborate on the peculiar structure of the MPS-TM's eigenspectrum and give several arguments for the close relation between the structure of the low energy spectrum of the system and the form of static correlation functions. Finally, we discuss how the MPS-TM connects to the exact Quantum Transfer Matrix (QTM) of the model at zero temperature. We present a renormalization group argument for obtaining finite bond dimension approximations of MPS, which allows to reinterpret variational MPS techniques (such as the Density Matrix Renormalization Group) as an application of Wilson's Numerical Renormalization Group along the virtual (imaginary time) dimension of the system.


I. INTRODUCTION
Determining the vacuum of an interacting field theory or the ground state of a strongly interacting quantum system described by a local translational invariant Hamiltonian is one of the most fundamental and challenging tasks in quantum many-body physics. Once obtainedpossibly in some variational way -how much information about the Hamiltonian is then encoded within the ground state? We will demonstrate that it is possible to extract many low-energy features of the Hamiltonian by just having access to the ground state. This is possible due to the Hamiltonian being a sum of (quasi-) local terms; this locality is the key to uncovering the mysteries of quantum many-body systems, such as the presence of a finite group velocity in quantum lattice systems, known as the Lieb-Robinson bound, and the relation between the spectral gap and correlation length [1][2][3]. The latter result connects a single characteristic of the static correlation functions of the ground state to one particular excitation energy.
This work continues along this line by investigating to what extent information about the full dispersion relations of the different elementary excitations of the model is encoded within the ground state and its correlations. Throughout the paper we assume translation-invariant Hamiltonians, such that excited states can always be characterized by momentum. Any statement regarding the spectrum of a Hamiltonian is to be interpreted up to an overall energy scaling and a constant energy shift. The shift is typically chosen such that the ground state energy E 0 = 0. The overall energy scale is represented by a characteristic velocity (e.g. the Lieb-Robinson velocity related to the norm of the Hamiltonian terms, or some spin-wave velocity) in the system.
In theory, the full dispersion relation can be reproduced from the ground state if the map between a local Hamiltonian and its corresponding ground state is bijective. For strictly n-local Hamiltonians, i.e. Hamiltonians for which every term acts only on a finite number n of neighboring sites, such a bijective relation is generically obtained. There the n-site reduced density matrices (RDMs) of ground states represent extreme points in the convex set of all possible n-site RDMs. The Hamiltonian can then be represented as a hyperplane in the space of such RDMs, and the energy will necessarily be minimized for an extreme point in this set. Each of these points uniquely determines an n-local parent Hamiltonian as the tangent space to the boundary at this point [6]. This argument is however of very limited practical use as it is computationally infeasible to characterize this convex set [7]. Also, the uniqueness is only obtained by restricting to a class of n-local Hamiltonians and there might exist other (n+k)-local (with k ≥ 1) or quasi-local Hamiltonians for which this is the exact ground state. One of the main goals of this paper is thus to identify which features of all those Hamiltonians can be captured in the ground state and its correlations.
We follow a more practical approach based on local information contained within the ground state, which is naturally accessible through a tensor network representations of the same. A central local object arising in tensor network simulations is the tensor network transfer matrix (defined in Section Sec. II). Indeed, the main motivation for this work originates from numerical results obtained from tensor network simulations of ground states of strictly local translation invariant Hamiltonians in the thermodynamic limit. There it is observed that the spectrum of the transfer matrix exhibits a very peculiar structure, from which certain information about the low energy excitation spectrum of the underlying Hamiltonian can be extracted. These results are presented and discussed in Sec. III for a set of prototypical quantum models on lattices in one and two dimensions, as well as (1+1)-dimensional field theories.
We provide several arguments for explaining these observations in Sec. IV. There we argue how the structure of the eigenvalue spectrum of the transfer matrix allows to reproduce the expected form of correlation functions in gapped quantum states and use the single mode approximation to relate these eigenvalues to excited states of the Hamiltonian. We also show the converse, i.e. that the excited states of the Hamiltonian affect the static correlations functions, either by employing arguments from relativistic theories or by using momentum filtering to refine the celebrated proof of Hastings in Ref. 3 for the relation between gap and correlation length.
Sec. V follows an alternative approach by directly connecting the transfer matrix in the context of tensor network states to the exact quantum transfer matrix (QTM) [8,9] at zero temperature, which appears in path-integral formulations of partition functions or ground states of quantum systems. Tensor network methods for studying such transfer matrices have been successful since the invention of the Transfer Matrix Renormalization Group (TMRG) method to simulate classical models in two [10] and higher dimensions [11]. Invoking a quantumto-classical mapping, this method has later been generalized and used to simulate one-dimensional quantum models at finite temperature [12][13][14] and recently to also include real time evolution [15,16]. In these methods, the object which is approximated by a tensor network is the (quantum) transfer matrix itself. In this work however, we investigate the transfer matrix at zero temperature generated by a tensor network approximation of the ground state. We also explain how the renormalization group (RG) allows to interpret the MPS-TM as a compressed version of the QTM. More specifically, in Sec. V C we demonstrate how Wilson's Numerical Renormalization Group (NRG) for impurity systems -or its recent reformulation using the Multi-scale Entanglement Renormalization Ansatz (MERA) -allows to build an MPS approximation of the ground state with finite bond dimension D from the QTM. This construction yields a novel connection between tensor network states and RG methods.

II. TENSOR NETWORK TRANSFER MATRICES
In this section we define the regular and mixed Transfer Matrix (TM) for Matrix Product States (MPS) [17][18][19][20] on one dimensional lattice systems and continu-ous Matrix Product States (cMPS) [21][22][23] on (1+1)dimensional field theories respectively. In the context of higher-dimensional lattice systems described by Projected Entangled Pair States (PEPS) [19,24] we consider two-dimensional lattice systems on cylinders. There we obtain an effective one-dimensional lattice system by blocking sites on rings around the cylinder as described in Sec. II A.
As we are interested in bulk properties of quantum systems, we will generally work in the thermodynamic limit, where for gapped one-dimensional quantum lattice systems, a good approximation of the ground state can be obtained by using a uniform MPS ansatz with finite bond dimension D where A sj is a set of d matrices ∈ C D×D containing all variational parameters defining the state, s j labels states within the d-dimensional local Hilbert space on each site and j labels sites on the lattice. v L and v R are boundary vectors which have no effect on bulk properties. An optimal MPS representation of the ground state can readily be calculated using variational uniform MPS techniques [25,26]. For ground states of higher-dimensional lattice systems similar techniques can be used for uniform PEPS [27]. Equivalently, ground states of (1+1)-dimensional field theories in the thermodynamic limit can be well approximated by uniform cMPS, where, e.g., a one-flavor bosonic cMPS of finite bond dimension D is given by where again Q, R ∈ C D×D contain all variational parameters defining the state. Hereψ † (x) are bosonic creation operators, P is the path ordering operator, |Ω is the vacuum of the field theory and v L and v R are again boundary vectors having no effect on bulk properties. To obtain cMPS ground state approximations, the algorithm of Ref. 26 can be adapted accordingly.

A. Regular Transfer Matrix
It is well known (and reiterated in Sec. IV A) that static correlation functions with respect to a uniform MPS ground state are obtained using the regular MPS transfer matrix (MPS-TM) [19], which is given by withĀ s the complex conjugate of A s . To simplify notation we will generally omit the subscript A denoting the MPS matrix whenever no confusion can arise. For continuum results we define the (generator of the) cMPS transfer matrix as where againQ andR denote the complex conjugates of Q and R respectively. We will again omit subscripts whenever they are not necessary. The relation to the lattice transfer matrix is given by with the lattice spacing of an underlying lattice discretization and T the transfer matrix of the corresponding MPS defined on the discretized lattice. Finally, for two-dimensional systems studied using PEPS, we work in the setting of infinitely long cylinders. By blocking the PEPS tensors A s udlr on a ring along the (finite) transversal y-direction of the cylinder we can then interpret this contracted object as a uniform MPS along the (infinite) longitudinal x-direction of the cylinder and we define the longitudinal transfer matrix as in Eq. (3). For a square lattice geometry this MPS has bond dimension D Ny and physical dimension d Ny where N y is the number of sites along the circumference of the cylinder. A graphical representation of the obtained TM is given in Fig. II A.

B. Symmetries and the Mixed Transfer Matrix
If a uniform MPS defined by a set of matrices A s 1 is invariant under a local unitary symmetry operation u, one can show [28] that where A s 2 defines the transformed state. Here V is a unitary gauge transformation on the auxiliary space and e iθ is the dominant eigenvalue with magnitude one of the mixed transfer matrix In fact, the MPS is invariant under the local symmetry u if and only if the spectral radius of the mixed transfer matrix ρ(T A2 A1 ) is one, i.e. the fidelity per lattice site is one. In the case of higher dimensional lattice systems, a relation similar to Eq. (6) holds for PEPS [28].
For phases with a spontaneously broken symmetry, the ground state is degenerate and the variationally best ground state approximations within the manifold of MPS of bond dimension D are minimally entangled states which exhibit maximal symmetry breaking. Such states can be transformed into each other by applying symmetry operations of the broken symmetry.
In the ground state of a one-dimensional quantum system, continuous symmetries for which the order parameter does not commute with the Hamiltonian cannot be spontaneously broken [29,30]. Nevertheless, close to or within a gapless phase with a continuous symmetry, it is sometimes energetically favorable for a variational (c)MPS approximation of the ground state to break this symmetry and to approximate an excited state with very small excitation energy and much smaller entanglement instead. This pseudo symmetry breaking is purely an effect of finite bond dimension and also gives rise to a pseudo order parameter [31,32]. The symmetry is restored in the limit D → ∞.
In a phase with broken symmetry on a lattice, let A s 1 and A s 2 be MPS approximations of two ground states with the same variational energy but different order parameters and maximally broken symmetry. The orthogonality of these states requires that the fidelity per lattice site must be strictly smaller than one, i.e. the spectral radius of the mixed transfer matrix ρ(T A2 A1 ) < 1. Equivalently, for continuum systems we define the (generator of the) mixed cMPS transfer matrix as where Q 1 , R 1 and Q 2 , R 2 are two different cMPS representations. Similar to the lattice case, if Q 1 , R 1 and Q 2 , R 2 describe two equally good ground state approximations (i.e. having the same variational energy) with different order parameters and maximum symmetry breaking, the spectrum P Q2,R2 Q1,R1 has strictly negative real parts.
The degeneracy of the ground state in phases with broken symmetries gives rise to topologically nontrivial excitations (kinks or domain walls), which typically correspond to the elementary excitations of the model. The mixed (c)MPS-TM of type (7) or (8) of these symmetry broken ground states plays a crucial role in obtaining a variational approximation for such excitations, whereas the regular (c)MPS-TM of type (3) or (4) is the central object for topologically trivial excitations [32,33].

III. NUMERICAL RESULTS
This section illustrates and discusses typical spectra of the regular and mixed (c)MPS-TM of obtained (c)MPS ground state approximations and compares it to low energy excitations for several one-dimensional quantum models of interest. For the eigenvalues of the (c)MPS-TM we write λ j = e −εj +iφj (9) where ε j = − log |λ j | and φ j = arg λ j is the complex argument. This form already suggests that the ε j will be related to some characteristic energies of the model, as motivated throughout the remainder of this paper. Low lying excitation energies for one-dimensional models are obtained by means of both topologically trivial and nontrivial variational uniform (c)MPS ansatzes [32,33] and -if applicable -are shown together with exact solutions.
For two-dimensional models we exploit the observed relation between eigenvalues of the transfer matrix and location and magnitude of energy dispersion relations to give a first estimate of the dispersion of elementary excitations.

A. One-dimensional Lattice Models
We will first focus on three prototypical onedimensional lattice models. We start with the spin-1/2 XY model in an external magnetic field which can be solved exactly in terms of free spinless fermions [34][35][36]. Here S α j denote spin-1/2 operators defined on site j. We consider the gapped ferromagnetic regime 0 < γ < 1 and 0 < g < 1, where the system is in a symmetry broken phase and the ground state is twofold degenerate with local order parameter m x = S x j . Here the elementary excitations are domainwall-like and therefore well approximated by a topologically nontrivial MPS ansatz. Specifically, we consider the incommensurate phase γ 2 + g 2 < 1, where correlations oscillate with arbitrary wave vectors.
As a second example we consider the spin-1/2 XXZ model in an external magnetic field This model is solvable as well and the ground state and elementary excitations in the thermodynamic limit can be obtained via Bethe ansatz [37][38][39]. Here we consider the antiferromagnetic gapless incommensurate phase specified by −1 < ∆ < 0 and 0 < |h| < 1 − ∆, where there are gapless excitations at multiples of the Fermi-momentum k F = ( 1 2 − m z )π with m z = S z j the ground state magnetization. In this phase there is no spontaneous symmetry breaking, however due to criticality the finite D MPS ground state approximation breaks the continuous rotational symmetry in the XY plane (c.f. Sec. II B). This makes it possible to use a topologically non-trivial variational MPS ansatz for excitations.
The last example being studied is the (extended) Fermi-Hubbard model where c † σ , c σ denote creation and annihilation operators of fermions with spin σ, n σ = c † σ c σ and n = n ↑ + n ↓ . For V = 0 this model is non-integrable. We consider the repulsive regime, where U, V > 0, away from half filling (µ = 0), which again corresponds to a gapless incommensurate phase. There is no spontaneous symmetry breaking in this phase and we consider topologically trivial excitations only.
In On the left we plot the eigenvalues λ j = e −εj +iφj of the regular MPS-TM on the complex plane within the unit circle, whereas on the right we plot ε j = − log |λ j | vs. complex argument φ j , along with the 4 lowest variational excitation energies obtained from a topologically trivial variational MPS ansatz. We do not plot the dominant eigenvalue λ 0 = 1. In Fig. 3 we show results for the mixed MPS-TM and topologically non-trivial excitations for the XY model and XXZ model only, with the same parameters as above. If available we also plot the exact dispersion of the elementary excitations as well as the lower boundaries of multi-particle continua for reference. From Fig. 3 it is apparent that the topologically nontrivial variational ansatz captures elementary excitations with high accuracy. The two particle continuum consists of combinations of two elementary excitations and is thus partially captured by a topologically trivial ansatz (c.f. Fig. 2), which is consistent with results in Ref. 32. The full continuum can be recovered by using a variational MPS ansatz including scattering states of elementary excitations [40]. Consequently, low lying odd particle states are partially captured by a topologically non-trivial ansatz, whereas low lying even particle states are partially captured by a topologically trivial ansatz. In the case where there are no topologically non-trivial excitations, there is no such distinction.
Concerning the eigenvalues of the MPS-TM we can now make the following remarkable observations. In the plots on the left of Fig. 2 and Fig. 3 we see that most of the eigenvalues arrange themselves along lines with approximately constant complex argument φ j . This fact is reflected in the arrangement of the ε j in columns in the plots on the right, where the low lying ε j correspond to eigenvalues λ j close to the unit circle.
We can also observe that the complex arguments φ j of the low lying ε j coincide very precisely with the momenta k min of the minima in the dispersion of excitation energies. The corresponding values of ε j are related to the minimum excitation energies E min via some characteristic velocities v α and can serve as a first approximation for this energy if the velocity can be estimated. It appears that this velocity, which determines the energy scale, can also vary among excitation minima within each respective model in the examples shown in Fig. 2 and Fig. 3. Indeed, for the XY model, the momenta k min of the  minima of the elementary excitations and the three particle continuum are well reproduced within less then 1% accuracy by the eigenvalues of the mixed MPS-TM with largest magnitude. Consequently the same holds for the minima of the two particle continuum and the regular MPS-TM. For the elementary excitations we estimate the characteristic velocity v 1 relating the lowest excitation energy E min and ε 1 = − log |λ 1 |, where λ 1 is the eigenvalue with second largest magnitude, as v 1 = Emin ε1 ≈ 0.9409 (c.f. Fig. 3 inset), where we have extrapolated the value of ε 1 for D → ∞. Towards the end of Sec. IV C we show how the value of this velocity can be estimated from assuming a Lorentz-invariant low energy behavior. There we obtain an estimate for v 1 which agrees with the value obtained above within 1% accuracy.
For the XXZ model, the momenta of the gapless excitations at multiples of k F are again very precisely reproduced by the arguments φ j of the eigenvalues of the regular and mixed MPS-TM with magnitude close to one (i.e. ε j close to zero) within less than 1% accuracy. Notice that in the limit D → ∞ we expect ε j → 0, i.e. the spectral radius of the mixed MPS-TM also becomes one and the rotational symmetry in the XY plane is restored (c.f. Sec. II B).
For the extended Hubbard model, the star like structure of the eigenvalue spectrum in Fig. 2 is very pronounced. In the right plot the ratios of the variational dispersion minima and the lowest ε j at k = φ j ≈ 0.4π and k = φ j ≈ π might suggest a characteristic velocity v α > 1. However, one would expect to have gapless excitations at these momenta, as well as at k ≈ 0.2π, too, suggesting that the corresponding variational energies are not converged. It is instructive to either use larger bond dimension or enhance by using an ansatz including scattering, which we have not performed here. It is however interesting to note for these minima that the accuracy of the variational energies and the low lying ε j appears to be roughly on the same level.
The above observations are truly remarkable. The mere knowledge of the ground state MPS-TM already yields important information about the excitation spectrum of the underlying Hamiltonian. More specifically, the momenta k min of the excitation energy minima in the respective particle sectors can be determined accurately and the corresponding energies can be estimated in a first approximation just from static ground state properties. An advantage of this approach over considering just static correlation functions is discussed in Sec. IV B.

B. (1+1)-dimensional Field Theories
We will now turn to continuous (1+1)-dimensional field theories and study the Lieb-Liniger model [41] using cMPS methods. The Hamiltonian is given by with repulsive interaction strength c > 0 and chemical potential µ > 0, where ψ and ψ † are bosonic field operators. The model depends only on a single parameter γ = c ρ , with ρ = ψ † ψ the ground state particle density and it is critical for all values of γ. Fig. 4 and Fig. 5 show results for the eigenvalues of the regular and mixed cMPS-TM similar to the lattice case for (a) γ ≈ 1.35 and (b) γ ≈ 311.5 and D = 64. Given the relation in Eq. (5), the right column now plots the real part of the eigenvalues −σ j of the generator P versus their imaginary part, which is now interpreted as momentum. In the left column the eigenvalues λ j of T = exp(P ) are plotted on the complex plane within the unit circle as in the lattice case. In the continuum setting, momentum space is no longer 2π-periodic and the definition of T is not fully justified, as it can come with any real power x ≥ 0 in continuum correlation functions, where only integer powers appear in lattice correlation functions. Nevertheless, it helps in illustrating that the spectrum of eigenvalues of the transfer matrix exhibits a similar structure. The fact that this structure in Fig. 4 and Fig. 5 is less outspoken than for some of the lattice models indicates a larger contribution of microscopic effects for this specific case. Next we will study the simplest Lorentz-invariant theory available, the free (1+1)-dimensional Klein-Gordon boson described by the Hamiltonian: where we have taken the speed of light to be unity. The field operators φ and π can be written in terms of the cMPS Fock space operators ψ and ψ † as: where an arbitrary scale ν is introduced [42]. The Hamiltonian (14) needs to be regularised, and this is achieved by adding the term to the Hamiltonian. We make the following observations regarding the eigenspectrum of the generator of the cMPS transfer matrix P (4) corresponding to H KG , as plotted in Fig. 6. The eigenvalues σ j of P are all real and negative, for all values of m and D. This reflects the fact that the relativistic dispersion relation has a single minimum at momentum zero. The eigenvalue with largest real part of P converges to −m as D → ∞. As is argued in Sec. V A, it follows from the Euclidean invariance of the quantum transfer matrix of a relativistic theory that in the limit D → ∞ the eigenspectrum of P should be the same as that of H KG (up to the sign), corresponding to a characteristic speed v = c = 1. The above observation provides numerical support for this using and extrapolating from finite D data.
We can also study the distribution of the eigenvalues of P as a function of the bond dimension. For any value of m, the (negative) eigenvalues σ j become dense in the region [m, +∞). In the gapped phase (m > 0), the density of eigenvalues diverges at m, i.e. the ratio of the nth largest and second largest eigenvalue of P converges to unity for low lying n. This is similar to the density of states in a gapped single particle excitation branch near the minimum of the dispersion relation. As the mass m is taken to zero the theory becomes critical, and the cMPS approximation enters the so-called "finite entanglement regime" [43]. The eigenvalues of P are still all real and converge to zero as D → ∞, but their ratios now converge to values larger than one, and are expected to encode universal data [42]. This is tantamount to the statement that the effect of finite bond dimension is to only introduce a single scale into the underlying conformal field theory, and implies that universal quantities can be extracted straightforwardly from cMPS data [42]. Let us exemplify this by spelling out the results obtained by taking m = 0.2 in (14) and using the modest range of bond dimensions up to D = 36. Scaling with 1/D and extrapolating to D → ∞ yields that the second largest eigenvalue of P tends to 0.201, thus reproducing the mass accurately. The ratio of the third and second eigenvalue of P is estimated to converge to 1.040 as D → ∞, and the ratio of the fourth and second to 1.086. We note that, since the theory is free, the deviation from unity can not be due to convergence to some bound state just above the lowest branch, and can only be an effect of numerical accuracy. The same ratios for m = 0 converge to approximately 2.51 and 3.1, respectively, and are related to properties of the excitation spectrum of the underlying conformal field theory [42].

C. Two-Dimensional Lattice Models
The observed connection between the eigenvalues of the transfer matrix and the minimum of the dispersion relation opens up a way to infer properties of the dispersion relation of two-dimensional systems, which are notoriously difficult to deal with. To this end, given a two-dimensional translation invariant Projected Entangled Pair State (PEPS) [24] on a square lattice cylin-der with periodic boundary conditions in y direction, we block all PEPS tensors in a ring around the cylinder (i.e., with the same x coordinate along the cylinder). We then consider the quasi-one-dimensional system along the cylinder obtained that way, described by blocked tensors A s udlr , and its transfer matrix T . As the original state was also translational invariant in y direction, we can label the eigenvectors of T by eigenvalues e iky of the action of the translation operator on the auxiliary degrees of freedom, as given by Eq. (6). On a hexagonal lattice cylinder, we additionally block two neighboring PEPS tensors in order to obtain a quasisquare lattice before further blocking all obtained tensors in a ring around the cylinder. Assuming that the observed connection between the leading eigenvalues of the transfer matrix and the minimum of the dispersion relation holds for each k y independently, we obtain the location and relative energy of the minima of the dispersion relation for each possible value of k y , which yields a cut through the dispersion relation. By closing the periodic boundaries in different ways, we can obtain this information along different symmetry axes, allowing us to reconstruct the overall form of the dispersion relation.
We now apply this strategy to the AKLT model [44] on the square and hexagonal lattice. Its ground state is constructed by placing spin-1 2 singlets on the edges of the lattice and projecting the spin-1 2 's at each vertex onto the sector with maximal spin (S phys = 2 for the square lattice and S phys = 3 2 for the hexagonal lattice); this construction yields the unique ground state of the SU(2) invariant Hamiltonian H = i,j Π i,j , where Π i,j is the projector onto the subspace with spin S = S phys on neighboring sites i and j. This construction corresponds exactly to a PEPS with bond dimension D = 2. Even though the exact ground state is known, little is known about the excited states, although recently an indirect method was proposed to estimate the gap by means of a Tensor Network Renormalization Group method [45].
We can use an iterative eigensolver to exactly determine the low-lying spectrum of T with very high precision on cylinders with sufficiently large circumference. Since the model possesses SU(2) symmetry, we can additionally label the eigenvectors of the transfer matrix by their spin (which corresponds to the spin of the excitation), thereby aiding the identification of different excitation branches. As the transfer matrix of the AKLT model is hermitian (up to a gauge transformation), its eigenvalues are real, and thus k x = 0, π. The pairs (k x , k y ) for all eigenvalues of T are therefore arranged along lines in the Brillouin zone; by properly closing the periodic boundaries, we can thus obtain data points along different symmetry axes. The results for the square and hexagonal lattice are shown in Fig. 7 and Fig. 8, respectively. In both cases, we find an isolated branch of antiferromagnetic spin-1 excitations, with a two-particle continuum starting at about twice the elementary quasienergy gap, in agreement with known results for onedimensional systems. In particular, for the square lattice of the AKLT model on the hexagonal lattice along the symmetry axes indicated in the inset. As for the square lattice, we find a branch of spin-1 excitations, whose minimum is now around the Γ point with momentum (0, 0), as the unit cell contains two spins. Again, there is a continuum of two-particle excitations at about twice the quasi-energy at Γ. The data has been obtained from cylinders with circumferences Ny = 6 √ 3 for Γ-M, and Ny = 12 for M-K and K-Γ in units of the lattice constant. Note that due to the choice of the unit cell, the Γ-M line is folded up to a smaller Brillouin zone; a possible unfolding is indicated by the empty circles. Also, note that the M-K line seems to correspond to a continuation of the Γ-K line into the third Brilloin zone, as indicated by the dashed blue line in the inset.
we find the minimum of the dispersion at momentum (k x , k y ) = (π, π), whereas for the hexagonal lattice the minimum is found at (k x , k y ) = (0, 0). For both lattices the minima appear on the isolated S = 1 branches.

IV. STATIC CORRELATION FUNCTIONS
In this section, we elaborate on the relation between the eigenvalues of the transfer matrix and static connected correlation functions and use this to provide several arguments for understanding the peculiar structure of the transfer matrix spectrum. Without loss of generality we consider the case of one-dimensional lattice systems and write the static connected correlation function for operators A 0 and B n acting on single sites 0 and n as where . . . denotes the expectation value with respect to the ground state. These arguments can readily be extended to operators acting on multiple sites, as well as higher dimensional systems and continuum systems. In Sec. IV A we explain how the clustering of the eigenvalues of the transfer matrix in branches allows to recover the typical Ornstein-Zernike form of correlations (to be defined below) in the limit D → ∞. Sec. IV B uses the single-mode approximation to relate these branches to minima in the dispersion relation and also discusses why generically the full spectrum of the transfer matrix can provide more information than selected typical static correlation functions.
We also investigate this connection in the other direction by showing how the low-energy excitations of the Hamiltonian affect the static correlation functions in the ground state. Sec. IV C assumes a Lorentz-invariant low energy description to recognize the structure of the spectrum of the transfer matrix as the finite D manifestation of the Källén-Lehmann representation of correlation functions. Finally, Sec. IV D uses momentum filtering to refine the proof of Ref. 3 for the relation between the correlation length and the energy gap of the system.

A. Recovering the Ornstein-Zernike Form
Let us first recall how the regular MPS-TM gives access to all static correlation functions in the corresponding MPS. For this we assume a complete eigendecomposition of the transfer matrix where we have dropped the subscript A for notational simplicity. Here |j) and (j| are the right and left eigenvectors of T respectively with (i|j) = δ ij and we again write λ j = e −εj +iφj for the eigenvalues. We demand the ground state MPS representation to be injective, such that there is a unique dominant eigenvalue λ 0 = 1 and |λ j>0 | < 1 for all other eigenvalues. Furthermore we define the operator transfer matrix (OTM) for operators O acting on n sites as where O i1...in j1...jn = i 1 . . . i n |O|j 1 . . . j n . With the above definitions it is well known that Eq. (17) can be written as (20) where we have defined the form factors As for a finite bond dimension Eq. (20) is a finite sum of exponentials, connected correlation functions for sufficiently large distances n must decay as a pure exponential C AB (n) ∼ exp(−n/ξ) where the correlation length ξ corresponds to the inverse of the smallest ε j with nonzero form factor. In contrast, typical correlation functions in gapped phases are expected to decay at large distances n as where there is an additional power law contribution to the decay with an exponent η, which in principle depends on the operators A and B. Close to a critical point, this form can be motivated from conformal field theory or from a general renormalization group argument. Approaching the critical point takes the correlation length ξ → ∞, and a pure power law decay of correlations remains, where the scaling exponent η can depend on the choice of operators A and B. Sufficiently deep in a gapped phase, on the other hand, Eq. (22) is known as the Ornstein-Zernike form and η typically depends on the number of spatial dimensions d as η = d/2 [46,47], i.e. for a one-dimensional quantum system for large distances n.
In the limit n → ∞, an MPS with finite bond dimension would correspond to η = 0. Nevertheless, for D → ∞ the correct form of the correlation functions should be restored, at least up to some length scale. We proceed by noting that the scaling form in Eq. (22)   2. On each of these lines the λ jα become sufficiently dense for D 1 and the corresponding ε jα will follow some dispersion. We will then reorder the indices j α , such that the ε jα are in ascending order and ∆ α = ε jα=0 is the minimum. In the complex plane this parameterization corresponds to going from the λ jα closest to the unit circle along the line with constant complex argument φ α towards the center. For small j α , ε jα changes smoothly as ε jα = ∆ α + g j κ α with some constants g, κ > 0, possibly depending on α.
3. The form factors f AB jα also vary smoothly and follow, to leading order, some power law f AB jα ∼ j ρ α where the exponent ρ depends on the operators A and B and possibly α. For ρ = 0 the leading order corresponds to a non-zero constant.
With the above assumptions and using the Euler-McLaurin formula to approximate the discrete sum with an integral, the contribution of one line α of eigenvalues in Eq. (20) becomes approximately where we have replaced j α with a continuous parameter z and integrated to ∞ for convenience, as e −ngz κ decays sufficiently fast with increasing z, even for moderate n.
In Fig. 9 we show corresponding numerical evidence for the XY model, defined in Eq. (10), at γ = 0.5 and g = 1.05, i.e. in the gapped paramagnetic phase, with bond dimension D = 32. In this phase all the eigenvalues of the transfer matrix are real, i.e. φ j = 0. We plot ε j along with a quadratic fit ε(j) = a j 2 + ∆, as well as the form factors f XX j , where X stands for S x , vs. index j for j < 7. Again, j labels eigenvalues for which the form factors f XX j are non-zero, in ascending order. The data confirms the expected values of the exponents as κ = 2 and ρ = 0, yielding the Ornstein-Zernike behavior of the correlation function expected for this model in this parameter regime [35].

B. Static Structure Factor and the Single Mode Approximation (SMA)
Whereas the previous subsection indicates why the peculiar structure of the eigenvalue spectrum of the MPS-TM allows to recover the typical form of static correlation functions for D → ∞, it makes no connection between the branches of eigenvalues appearing in this spectrum and the dispersion relations of the elementary excitations of the model. To make this connection more explicit, we now reiterate the well-known result that the Single Mode Approximation (SMA) produces dispersion relations which are strongly dependent on the static structure factor [48][49][50][51][52]. In particular, we will illustrate why the energies E(k) of the lowest energy-momentum eigenstates |E α k of a local translation invariant Hamiltonian H = n h n become very small at momenta k where the TM has eigenvalues λ j = e −εj +iφj with ε j approaching zero and ±k = φ j .
The generality of this subsection is based on the recent proof that elementary excitations on top of a gapped strongly correlated ground state tend to be localized or particle-like [53]. This means that there exists a representation of all lowest lying excited states by acting with the Fourier transform of a quasi-local operator O ( ) centered around site n = 0 and having support on sites n ∈ [− , ], on the ground state: where O ( ) n = U n O ( ) U † n with U n the lattice translation operator over n sites and V is the volume of the system. For the remainder of this subsection, we work in a finite system with periodic boundary conditions in order to be able to define normalizable energy-momentum states.
The unnormalized state in Eq. (27) becomes exponentially close to a true isolated energy-momentum eigenstate |E α k with increasing linear size of the support of O ( ) [53].
We can readily use this set of states in a variational ansatz for low-energy excited states and we will consider their energy expectation value which is a good approximation for the lowest excitation energies of H at momentum k. The operator O ( ) should have zero vacuum expectation value, and can be chosen hermitian [73]. Together with the assumption of parity invariance for both H and O ( ) , it is well known [51] that Eq. (28) can be rewritten as where O k = 1 √ V n e ikn O n and we have omitted the superscript ( ) for notational simplicity. Here F (k) is the double commutator expectation value known as the oscillator strength and S(k) = S OO (k) is the static structure factor, which is related to the static correlation function by a Fourier transform, i.e.
In Eq. (29), we have again discarded subscripts denoting the operators O for notational simplicity. The static structure factor S(k) can however become very large. This will happen when the momentum k corresponds to the period of an oscillating static correlation function C(n) with very large correlation length. For these momenta it should thus be possible to optimize over possible operators O such that the resulting energy expectation value is very small and thus, by virtue of the variational principle, there exist excitations with small energies.
Generically, there is a one-to-one correspondence between the momenta k where S(k) peaks and the complex arguments φ j of the transfer matrix eigenvalues. This can easily be shown for the case of matrix product states where we use the same notation as in subsection IV A. Due to translation invariance the static structure factor can be written as This sum can be split into a part |x| ≤ 2l where operators O 0 and O n overlap and a part |x| > 2l where they commute. The first part can be bounded as The remaining part can be written as where J O is the operator transfer matrix defined in Eq. (19). To perform the geometric sum in the last line, we define the projector Q = 1 − |0) (0| and E = We can then write T n = QE n + |0) (0| = E n Q + |0) (0|, where the second term will not contribute to S(k) due to the zero vacuum expectation value of O. Assuming that there are no other eigenvalues with magnitude one -which is guaranteed if the MPS is injective -we can now safely perform the geometric sum and obtain where the form factors f OO j are given by Eq. (21). It is apparent that the fraction can become very large with k approaching the argument ±φ j [74] of an eigenvalue of T with ε j close to zero, i.e. magnitude close to one. In the case of symmetry breaking and topologically non-trivial excitations, similar arguments lead to the same form of S(k) as in (35), where E is replaced with the mixed MPS-TM [75].
As a final point, we elaborate further on the relation between the static structure factor and the spectrum of the transfer matrix. For an injective MPS, it will always be possible to find an operator, whose OTM exactly excites one or more of the eigenvalues of the TM on a branch with fixed φ α , i.e. (j α | J O |0) = 0. This, however, could be an operator with very large support. Conversely, it is possible that for operators with small support, several eigenvalues are excited on branches with different arguments φ α which are close together. For these operators, the static structure factor might then have a maximum at a momentum k which does not exactly correspond to one of the arguments φ α .
An example for this is the transition from commensurate to incommensurate order in the bilinear-biquadratic S = 1 Heisenberg chain H BLBQ = n cos(θ) S n · S n+1 + sin(θ) ( S n · S n+1 ) 2 . (36) where Z stands for S z (purple line) and variationally obtained dispersion E(k) (blue line) vs. momentum k along with εj vs. φj. Whereas S(k) still has its maximum at momentum π in (b), the minimum kmin of the dispersion has already started shifting away from kmin = π. In the spectrum of the transfer matrix this is reflected by the fact that the eigenvalues with largest magnitude have finished aligning along φα ≈ 0.755π. The maximum of S(k) doesn't start shifting until around θ ≈ 0.1314π.
We consider the regime 0 < θ < 0.25π, where the oscillation period of static correlation functions changes from π to some incommensurate period exactly at the AKLTpoint θ VBS = arctan(1/3) ≈ 0.1024π [54]. However by looking at the static structure factor of a simple one-site operator such as S z , the peak stays at k = π until some significantly larger valueθ ≈ 0.1314π. Based on the relation between the structure factor and the minima of the dispersion relation, we have reason to expect that the dispersion relation starts shifting away from π at this latter value of θ. It appears, however, that the minimum of the dispersion relation starts shifting at θ ≈ 0.12π, a value in between θ VBS andθ. In Fig. 10 we observe from the full spectrum of the transfer matrix, that this happens when the eigenvalues with largest magnitude have finished aligning along the line of constant phase φ α ≈ 0.755π. At this point it appears that the support of the operator O generating the excitation, has shifted from the -by now very small -branch with φ α = π to the now fully aligned branch with φ α ≈ 0.755π. This value of φ α further shifts towards 2π/3 with θ → 0.25π, where the gap then closes at momenta k = 0, ±2π/3 [55][56][57].
From this example it is apparent that especially in the vicinity of such peculiar crossover points, it is worthwhile to look at the full spectrum of the transfer matrix. It generically contains more information than a simple static structure factor, as the transfer matrix is completely independent of the choice of operator. It illustrates that there is a crossover regime where the MPS-TM has to react by developing additional branches with constant φ α and the support of the operator O generating the excitation has to shift to this newly developed branch. It appears that this process is completed after the eigenvalues have fully aligned along the newly developed branch. It would be interesting to investigate this process further with the precise knowledge about the operator O generating the excitation and how the form factors (j α | J O |0) develop with θ.
For now we conclude with the observation that there are indeed also situations where the locations of the lowest lying TM eigenvalues do not precisely coincide with the minima of the dispersion. These are however very special cases like the discussed crossover from commensurate to incommensurate order in the bilinear-biquadratic S = 1 Heisenberg chain, where the TM has to adapt to changing conditions within some crossover regime. There the peculiar structure of the TM eigenvalue spectrum is however also different from normal situations as investigated in Sec. III, which can serve as an indicator for such exceptional situations.

C. Källén-Lehmann representation
In the previous two subsections, we have discussed two effects of the peculiar distribution of the eigenvalues of the MPS-TM. Firstly, the clustering of eigenvalues onto lines starting from the origin allows to recover the typical form of static correlation functions in gapped quantum ground states in the limit D → ∞. Secondly, this distribution causes peaks in the static structure factor, which can be related to minima in the dispersion relation of excitations using the single mode approximation. In retrospect, the first effect only requires a dense distribution of eigenvalues along a line, without any connection between the location of these lines and the dispersion relation of the excitations of the system. The second argument only requires the existence of a single dominant eigenvalue with a phase corresponding to the momentum of the minimum of the dispersion relation, and does not explain why there needs to be a dense distribution of eigenvalues. This leaves open the question whether a different structure of the eigenvalue distribution could give rise to similar effects. Put differently, we would like to answer the reverse question, i.e. to what extent are static correlation functions and the clustering of the eigenvalues of the MPS-TM determined by the excited states of the Hamiltonian.
For Lorentz-invariant field theories, the Källén-Lehmann representation of two-point correlation functions provides such a direct connection to the excitation spectrum of the Hamiltonian. Whereas the Källén-Lehmann representation exists for arbitrary dynamical correlation functions, we here specialize to the case of static correlation functions between two scalar operators A(x) and B(y), where it is given by where, contrary to standard field theory notation, x and y denote spatial vectors, k denotes momentum and d represents the number of spatial dimensions. The first integral is over all possible masses in the theory, and ρ(M 2 ) corresponds to the density of states. If the lowest lying excitations correspond to single particle excitations with discrete masses M α , then ρ(M 2 ) will contain a contribution α δ(M 2 − M 2 α ). The state |M 2 α , 0 corresponds to the presence of such an excitation with mass M α and momentum zero.
Let us again restrict to the case of d = 1. We are interested in the long-range behavior of correlation functions in a lattice model, which is clearly dictated by the lowenergy behavior of the model. If this low-energy behavior can be captured by a Lorentz-invariant theory with masses M α , a remnant of the Källén-Lehmann representation of correlation functions should exist in the lattice model. For the field theory, the minimum of all dispersion relations is at momentum zero. However, in many cases taking the continuum limit of a lattice theory requires that e.g. N sites are blocked, and a single lattice dispersion relation with several minima for momenta φ α = 2πα N with α = 0, . . . , N − 1 gives rise to N independent dispersion relations of the field theory. A prototypical example for N = 2 is the XX-model, which corresponds to the staggered fermion discretization of relativistic Dirac fermions, where the two components of the Dirac spinor are put on even and odd sites respectively [58].
Eq. (37) represents the static correlation function in real space as an integral over momentum space. The corresponding representation for the static structure factor in momentum space can therefore easily be identified. In the Lorentz-invariant case, the single particle excitations add a contribution to the static correlation function that is given by a constant times the inverse of the dispersion relation of that excitation. For typical operators, these will be the dominant contributions. If the low-energy behavior of a lattice model is Lorentz-invariant, we can thus expect that S(k) should also receive contributions of the form c α Here, c α is a constant depending on the choice of operators A and B, m α is the mass of the excitations in units of the inverse lattice spacing and v α is the characteristic velocity, which for a proper Lorentz-invariant low-energy behavior should be the same (and thus independent of α) for all excitations. Since this form is only expected to hold for k ≈ φ α we can instead choose a 2π-periodic version of the dispersion relation to get contributions of the form We can now transform from the 2π-periodic variable k to the complex variable z = e ik and we change the notation of the static structure factor as S(k) → S(z = e ik ) to write where the contour integral is over the unit circle C, which is where S(z) is originally defined (see Fig. 11).
If we now were to construct an analytic continuation of S(z), we expect that every contribution of the form of Eq. (39) produces a square root singularity at z = 0 and inverse square root singularities at the points with ∆ α = m α /v α . We can thus choose the branch cuts to go from z (−) α to 0 and from z + α to +∞. Assuming that there are no other singularities, we can then deform the integration contour as in Fig. 11. In the limit η, η → 0 the arc segments A α and A α,β do not contribute, as they correspond to square root and inverse square root Figure 11: Change of integration contour: the original contour C corresponds to the unit circle, and is mapped to a contour C consisting of line sigments Lα and L α at distance η from every branch cut between z (−) α and 0, as well as arc segments Aα rotating by π around z (−) α with radius η, and arc segments A α,β rotating by 2π/N around the origin at radius η , between the branch cuts corresponding to z where we have transformed to y = ze −iφα and the contributions to the correlation function C(n) are now written as integrals over 0 ≤ y ≤ e −∆α along the branch cuts. It is apparent that a finite D MPS-approximation of the static correlation function tries to reproduce this continuum form with a discrete sum over the eigenvalues of the transfer matrix. These eigenvalues cluster on the branch cuts of the static structure factor, which we have related to the single particle excitation spectrum by assuming a Lorentz-invariant low-energy behavior.
For large n, where the low-energy contributions are dominating and the approximations are valid, we can again use the saddle point approximation around the point y = e −∆α and change variables to y = e −(∆α+x) , to obtain for the contributions to C(n), up to some factor which reproduces correlations of the Ornstein-Zernike form.
It is interesting to compare this to the discussion of the previous section. The single mode approximation for the dispersion relation E(k) was given in Eq. (29) as E(k) = F (k)/[2S(k)], which allowed to conclude that peaks in the static structure factor S(k) produces minima in the corresponding dispersion relation. Assuming Lorentz-invariance, the Källén-Lehmann produces contributions to S(k) of the form c/E(k) with c a kindependent constant, which looks like the reverse relation: minima in the dispersion relation give rise to peaks in the static structure factor. However, the single mode approximation also allows for minima which are not caused by S(k) but rather by a small or vanishing value for the oscillator strength F (k). It would be interesting if one could show that such excitations are necessarily related to low-energy features which have an intrinsically non-relativistic description, such as the gapless spin waves with quadratic dispersion relation in the ferromagnet.
As a final justification for the argumention of this subsection, we apply the above results to the case of the XY-model in the incommensurate gapped phase investigated in Sec. III A. There it is observed that for the elementary excitations the smallest excitation energy E min (i.e. the energy gap) and the eigenvalue of the transfer matrix with second largest magnitude λ 1 are related by some characteristic velocity v 1 via E min = v 1 ε 1 , where ε 1 = − log |λ 1 | and ε 1 = 1 ξ with ξ the correlation length ξ as established in Sec. IV A. On the other hand, in Eq. (41) we have defined ∆ α = mα vα , where for the dominant αbranch ∆ 1 is the inverse of the correlation length and we can thus identify ∆ 1 = ε 1 . We can now approximate the dispersion of elementary excitations E(k) around its minimum E min = E(k min ) by a Lorentz-invariant form. By identifying the mass m 1 with E min and thus noting that v 1 = m1 ∆1 = Emin ε1 we can then approximate for k ≈ k min . This suggests that the characteristic velocity relating E min and ε 1 is exactly the velocity of the Lorentz-invariant dispersion and can be estimated from the second derivative of the exact dispersion as For the parameters γ = 0.3 and g = 0.2 considered in Sec. III A this yields an estimate of v 1 = 0.9306, which differs from the value v 1 = 0.9409 obtained from v 1 = Emin ε1 by only ≈ 1%.

D. Momentum Resolved Relation between Correlation Length and Gap
In this section we derive more rigorous statements connecting the decay of static connected correlation functions to the dispersion E(k) of low energy excitations of local translation invariant Hamiltonians in the thermodynamic limit. We generalize the seminal work of Hastings [3] in which it is proven that the inverse of the energy gap ∆ of a local translation invariant Hamiltonian times a constant serves as an upper bound for the correlation length ξ of connected static correlation functions. This implies that if the gap vanishes, these correlations may (and in most cases will) be long ranged with a diverging correlation length.
The proof gives a statement relating the smallest overall excitation energy and the largest correlation length in the system, but does not take into account momentum information. In this work we extend the results in Ref. 3 to derive bounds on the decay of momentumfiltered correlation functions and relate them to the dispersion E(k) of low energy excitations. Specifically, we show that the inverse of the energy gap E(k) at a specific momentum k times a constant serves as an upper bound for a momentum-resolved correlation length ξ k . Conversely, the existence of a finite correlation length ξ k thus implies an upper bound for the energy E(k).
The detailed derivation of the bound in the most general setting is given in Appendix A; here we present the result for one-dimensional lattice systems and operators acting on single sites for the sake of simplicity. We start from the static connected correlation function of operators A i=0 and B j= >0 , for which we assume zero vacuum expectation value. We now attempt to extract momentum-space information while retaining real space information by replacing operator B at site with a gaussian wave packet centered around site , defined as 2r e ikn B +n (45) where N r is a normalization constant. We define the momentum-filtered correlation function as which corresponds to a Fourier transform of the product of the static correlation function and a gaussian wave packet centered around site . In momentum space this yields the convolution of the static structure factor S(k) and another gaussian wave packet in momentum space. In Appendix A it is proven that by tuning r as a fraction of , the momentum-filtered correlation function C k ( ) can, for sufficiently large , be bounded by where c 1 and c 2 are some constants and the inverse of c 2 is an upper bound for the momentum-filtered correlation length ξ k , given by Here v LR is the characteristic Lieb-Robinson velocity [1] and E * (k, δ) = min |k−k |≤δ E(k ) is the minimum of the dispersion E(k) in an interval around k given by δ.
The constant δ is introduced in the proof and and can be tuned to obtain the sharpest bound. From Eq. (48) it is clear that there is tradeoff, since increasing δ generally leads to a decrease of E * (k, δ). If k however corresponds to a minimum in the dispersion E(k) then the function E * (k, δ) is largely insensitive to δ in some region around the minimum and we can choose δ as large as possible within this region. If k on the other hand corresponds to a regular point where dE dk = 0 then there is a direct effect from increasing δ to decreasing E * . An optimal choice of δ is thus dependent on the form of E(k).
Colloquially, the above result means that the decay of the momentum-filtered correlation functions is dictated by the corresponding low energy states around that given momentum. Equivalently, one can say that a large correlation length for a given momentum k implies a small excitation energy E(k).
In practice, this means that one can deduce an upper bound for the low-energy spectrum of the Hamiltonian by looking at the momentum dependence of momentumfiltered correlation functions, where one would have to optimize over the parameter δ [76].

V. THE QUANTUM TRANSFER MATRIX (QTM)
This section uses renormalization group arguments to establish a close relationship between the MPS-TM and the exact Quantum Transfer Matrix of the model, thus providing a direct connection between the MPS-TM and the spectral properties of the Hamiltonian.

A. Imaginary Time Evolution as Tensor Network
Consider a one-dimensional local lattice Hamiltonian H with translation invariance and a unique ground state |ψ 0 with ground state energy E 0 = 0. To avoid issues with the infrared orthogonality catastrophe, we work in a finite system of N sites with periodic boundary conditions for this section. The ground state can be obtained by an evolution in imaginary time β with |φ init some initial state which is non-orthogonal to |ψ 0 . Since H is a sum of local terms h n , we can break β into small imaginary-time steps δ and use a Suzuki-Trotter decomposition of e −δH ≈ n e −δhn [8,59]. There are several strategies to write this (or alternative) decomposition(s) as a 2-dimensional tensor network with translation invariance in the spatial direction [16,60,61]. Obviously, all the information of H is thus encoded into this tensor network. If the state |φ init is initially in the form of a translation invariant MPS, then we can also interpret |ψ 0 as a translation invariant MPS by grouping contractions along imaginary time. Each column of the tensor network is an MPS matrix A s , itself being a half-infinite matrix product operator (MPO). A s then exactly represents the ground state up to a Trotter error [77]. A graphical representation of this construction is given in Fig. 12   We immediately see that ψ 0 |ψ 0 = Tr T N A where the MPS transfer-matrix T A was defined in Eq. (3). For a system with a unique ground state, the boundary conditions at β = 0 are irrelevant and an equivalent network with periodic boundary conditions in the temporal direction would be obtained for the thermal partition function Z β = Tr e −βH in the limit β → ∞. The exact MPS-TM for the ground state |ψ 0 thus corresponds to the quantum transfer matrix (QTM) at zero temperature, defined in Refs. 8, 9. It is important to note that all the information about the QTM -in particular its eigenvalues -is thus contained within the ground state |ψ 0 and its exact MPS representation A s . Note, however, that this exact representation with exponentially diverging bond dimension differs from a finite D approximationÃ s that can be obtained for instance from some variational algorithm. In Sec. V C we present a construction how such an approximation, which only retains degrees of freedom relevant for the physical degree of freedom s, can be obtained from this exact MPS representation A s and thus from the true QTM. To understand how this relates the MPS-TM to the Hamiltonian of the system, we first need to discuss how the latter relates to the exact QTM, which is the topic of the next subsection.

B. The QTM and the Hamiltonian
To understand what information about the underlying Hamiltonian can be extracted from the knowledge about the QTM, we first consider the case of relativistic (1+1)dimensional field-theories. The vacuum or ground state can be expressed in terms of a path integral formulation very similar to Fig. 12, but where both real space and imaginary time are continuous. We can then identify a factor e −δH for some infinitesimally small δ with a narrow horizontal slice of the entire network. Due to relativistic invariance, however, such slices are invariant under Euclidean rotations between real space and imaginary time. This means in particular that a vertical slicecorresponding to a field theory analogue of T -is also equal to e −δH . Hence, knowing the spectrum of T immediately yields knowledge about the spectrum of H. This in turn implies that all the information about the eigenvalues of H is already contained within the ground state. Of course, Lorentz-invariance also strongly restricts the dispersion relations of the theory, such that they are completely determined by a single parameter, being the mass of the excitation.
The relevant question is thus how much of Euclidean invariance between real space and imaginary time remains in non-relativistic lattice systems. As a concrete example we consider the one-dimensional XYZ model on a chain with N sites and periodic boundary conditions By invoking a Suzuki-Trotter decomposition [8,59] the finite temperature partition function Z β of this model can be mapped onto a classical 2-dimensional eight-vertex model [39,62] of dimension N × M , where the statistical weights W of arrow configurations depend on the model parameters, inverse temperature β, the Trotter number M and the type of chosen decomposition.
Here we use a generalized Suzuki-Trotter decomposition e −δH = e −δ n hn ≈ n e −δhn , with δ = β M . This leads to a real space decomposition introduced by Suzuki [8], where the quantum transfer matrix T is translation invariant on a slanted lattice [9] and we assume M to be a multiple of N . The statistical weights are then given by W abcd = cd|e −δhn |ab [39]. The same network is obtained after a rotation of the lattice by 90 degrees and considering an equivalent effective model on a chain of M sites and Trotter number N at an effective inverse temperatureβ. The QTM of the original lattice can therefore be written in terms of the effective model Hamiltoniañ H as T = exp(−β NH ) and effective statistical weights W up to a Trotter error in N . The effective model parameters and the effective inverse temperatureβ can be obtained from the vertex weights after rotating, i.e. from W bdac = W abcd . For a graphical representation of this mapping see Fig. 13. From this relation the spectrum of the effective Hamiltonian can in principle be studied by looking at the QTM generated by the ground state of the original Hamiltonian.
Note that depending on the model and parameter regime, the effective parameters may also become complex, resulting in a non-hermitian effective Hamiltonian and thus accounting for non-hermitian transfer matrices. This particular fact was already discussed in the case of systems with Lorentz-invariance, as the continuum limit can often only be taken after blocking N sites. In that case, all eigenvalues of the transfer matrix would have phases φ α = 2πn/N with n = 0, . . . , N − 1, and by defin-ingH = − log(T N ) up to some energy scale, a hermitian effective Hamiltonian would be obtained. For systems with incommensurate order, where the eigenvalues can have arbitrary phases which are not fractions of 2π, this is no longer possible.

C. Truncation of the Virtual System
In this subsection we show how to obtain an MPS ap-proximationÃ s with finite bond dimension from the exact QTM constructed in the previous subsections. As was shown there, an exact MPS representation of the ground state can be constructed from imaginary time evolution, where the MPS matrices A s are given as a semi-infinite MPO with exponentially diverging virtual dimension (c.f. Fig. 12). This construction allows to identify the MPS-TM with the exact QTM at zero temperature.
In the following we assume that the QTM can be written in terms of an effective local HamiltonianH = nh n as T = e −H , for instance via the construction of the previous subsection. In this representation, the MPS-TM appears to be completely translation invariant in the imaginary time direction, with no special role being played by the physical index s, where the matrices A s and We can then interchange the roles of real space and imaginary time and give two distinct new interpretations to the expectation value Ψ|O(x)|Ψ . Firstly, we can interpret the new imaginary time direction as the evolution of a pure state of an infinite one-dimensional virtual system in the x direction according to the MPS-TM T or equivalently, the corresponding HamiltonianH. At certain "times" x, there is an insertion of some impurity O at the fixed coordinate τ = 0, which destroys the translation invariance of the virtual system. An alternative interpretation is given below while constructing an approximation with finite bond dimension of the exact ground state MPS A s .
To arrive at a finite virtual dimension D, it is necessary to restrict the exponentially diverging amount of virtual degrees of freedom (DOF) of A s to a finite subset that is relevant for the physical DOF s at τ = 0. We thus need to identify the relevant low-energy subspace of the virtual system to describe the evolution of the impurity at position τ = 0 as a function of x; the finite D approx-imationÃ s is then obtained from the exact MPS A s by projecting onto this relevant subspace.
The relevant low-energy subspace for an impurity problem can be obtained by applying real space RG transformations, as was first shown in the seminal work of Wilson using his numerical renormalization group (NRG) [63]. However, we here follow the more recent construction using the multiscale entanglement renormalization ansatz (MERA) [64][65][66]. This approach allows to identify the relevant degrees of freedom as those living at the causal cone of the impurity [67,68]. For completeness, we repeat this argument for our specific case.  The following construction starts from the assumption that there exists a sequence of real space RG transformations U r that renormalizesH onto its low-energy subspace, where r labels the layers of successively applied RG transformations. Equivalently, U r renormalizes T onto the subspace spanned by its dominant eigenvectors. For concreteness, we consider a real space RG procedure that coarse grains four neighboring sites into two renormalized sites, which can be realized e.g. by a modified binary MERA [66]. IfH is scale-invariant, the transformations become independent of the layer index r after some amount of initial layers r * . At this point all RG irrelevant terms have been removed and the renormalized Hamiltonian is a fixed point of the scale invariant RG transformation U. For non-scale-invariant Hamiltonians the series of RG transformations terminates after r max ≈ log(R) layers, where R is the dominant length scale of the Hamiltonian. At this point there are no relevant DOF left. These two cases can e.g. be represented by a scale-invariant or a finite-range MERA respectively [65], an example of which is shown in Fig. 14.
The MERA construction allows to conclude that any perturbation at position τ = 0 can only affect the degrees of freedom living at its causal cone [68]. This causal cone is shown as red solid lines in Fig. 14. In particular, let us focus on one side of this causal cone, e.g. the lower half of Fig. 14. It is apparent that only one index of the RG network protrudes the boundary of the causal cone for each layer. We can therefore interpret these legs as the sites of an effective lattice system L W defined along the boundary of the causal cone, which we call the Wilson chain. We label the sites along this chain by the layer index r. Note that every the r-th site on this chain is an effective renormalized description of N r ≈ 2 r sites of the original lattice, reminiscent of the logarithmic discretization introduced by Wilson, hence the chosen nomenclature [63,67,68].
Contracting the RG network outside the causal cone (non-shaded region in Fig. 14) allows to renormalize the MPS-TM T A into a new transfer matrix T W along the Wilson chain L W [67,68]. It is immediately clear that T W corresponds to the transfer matrix TÃ of a new MPS with matricesÃ s , which are obtained by projecting the exact MPS representation A s exactly onto this subspace of relevant DOF along the Wilson chain. For non-critical systems, the Wilson chain is of finite length L and we obtain a virtual system with finite dimension D ≈ χ rmax , where χ is the bond dimension of the RG network. For critical, scale-invariant systems however, the Wilson chain is still infinite and any truncation to a finite system introduces some error.
It is of course well known that ground states with finite correlation length ξ can be well approximated by a finite D MPS [6]. By virtue of the previous sections, the effective HamiltonianH is gapped, where the gap is exactly equal to ξ −1 . It is interesting to contrast this interpretation with the more usual considerations regarding the gap of the actual physical Hamiltonian H of the system. For Lorentz-invariant systems, H andH are equal and the correlation length is directly given by the inverse of the gap in units where the speed of light c = 1. Without Lorentz-invariance, there is no obvious relation between H and its Euclidean rotationH. The result of Ref. 3 allows to bound ξ (and thus the gap ofH) in terms of the gap ∆ of H. When the low-energy behavior of H still allows for an effective relativistic description, it can indeed be expected that ξ ∼ v/∆ where v replaces the speed of light with some characteristic speed of the system. However, there is also the possibility of an inherently nonrelativistic low energy behavior, resulting in e.g. gapless hamiltonians H for which the ground state is an exact MPS [69] so that the correspondingH is gapped.
Finally, we can obtain an alternative interpretation by applying the Jamio lkowski isomorphism to map the pure state defined on the infinite one-dimensional virtual system of the previous discussion to the density matrix of a half-infinite one-dimensional virtual system with a boundary at τ = 0. By again interchanging the roles of real space and imaginary time, this half-infinite system undergoes dissipative evolution in the new timedirection x, corresponding to a Hamiltonian containing the terms ofH with support on the region τ > 0, and additional Lindblad operators corresponding to the action ofH across the boundary and the additional action of an operator O inserted at certain "times" x. Hence, all the Lindblad operators are acting near the boundary, and the truncation of the bond dimension corresponds to selecting the relevant degrees of freedom to describe the boundary of the system. It is a virtue of the MERA that it naturally unifies the process of selecting the relevant degrees of freedom for boundaries and impurities [68]. The resulting virtual system becomes effectively zero-dimensional, and can be interpreted as providing a holographic description of the physical system [22].

VI. CONCLUSIONS
In this paper we have investigated how much information about the excitation spectrum of a local translation invariant Hamiltonian can be obtained from local information and static correlations in the ground state. We have approached this question using the formalism of tensor network states in particular, but have also established several general results not restricted to tensor network formulations.
We have started by defining the regular and mixed tensor network transfer matrix for lattice and continuum models in Sec. II A and Sec. II B. We have then obtained tensor network approximations for the ground states of various prototypical quantum models on a lattice in one and two dimensions and (1+1)-dimensional field theories and studied the spectrum of eigenvalues λ of the tensor network transfer matrix in Sec. III. We have observed that the the complex arguments φ of the dominant eigenvalues correspond to the momenta k min of the minima in the low energy dispersion of the system. Especially for critical models one can therefore easily determine the momenta for which there exist gapless excitations directly from ground state properties. We have also related the logarithm of the absolute values ε = − log |λ| to the minimum excitation energies E min by some characteristic velocity E min = v ε, which can be estimated by e.g. assuming a Lorentz-invariant low energy behavior (c.f. Sec. IV C).
These observations are of large practical importance in the context of simulating quantum many body systems using tensor network techniques: it implies that already fairly accurate information about the structure of the low energy excitation spectrum can be obtained just from a variational ground state calculation. We have demonstrated that this is especially useful for two-dimensional systems using the PEPS formalism, for which no other efficient methods are presently known to extract information about excited states beyond the value of the gap. In particular we have investigated the AKLT model on a square and hexagonal lattice cylinder, where we have obtained a first approximation of the dispersion of the elementary excitations, for which there currently exist no other competitive numerical methods.
In Sec. IV we have gathered several arguments to explain how the eigenvalue spectrum of the transfer matrix affect the low-energy excitations of the model and vice versa. We have explained how a clustering of eigenvalues along lines of constant complex phase allows to recover the Ornstein-Zernike form of correlations in gapped phases if the distribution of eigenvalues becomes sufficiently dense in the limit D → ∞ in Sec. IV A. Using the single mode approximation and a recent proof about the locality of elementary excitations [53], we have argued why the phase of the eigenvalues along these lines correspond to minima of the dispersion relation of the elementary excitations in those models in Sec. IV B. There we have also discussed that the full spectrum of the transfer matrix generally contains more information than static correlation functions of specific operators, as the transfer matrix is completely independent of the choice of operators. To approach this connection from the reverse direction in Sec. IV C, we have first called on the assumption of a Lorentz-invariant low-energy behavior to identify the spectrum of the transfer matrix with a discrete version of the Källén-Lehmann representation of correlation functions. Finally, in Sec. IV D we have introduced momentum-resolved correlation functions by defining gaussian wave-packets of operators centered around a certain momentum k in order to refine the celebrated result of Hastings [3] relating the correlation length in the system to the gap of the Hamiltonian.
Furthermore we have identified the tensor network transfer matrix of an exact tensor network representation of the ground state with the quantum transfer matrix, appearing in path integral formulations of partition functions or ground states of quantum many-body systems, in Sec. V A. We have then argued how the quantum transfer matrix is related to the original Hamiltonian for system with and without Lorentz-invariance in Sec. V B. We have demonstrated that for systems without Lorentz-invariance the quantum transfer matrix can be written in terms of an effective Hamiltonian with effective parameters which are related to the original Hamiltonian. These parameters can in principle also be complex, thus yielding non-hermitian effective Hamiltonians and non-hermitian transfer matrices. For systems with commensurate order, a hermitian effective Hamiltonian can be obtained by blocking several sites in constructing the transfer matrix.
As a final point we have demonstrated in Sec. V C how a tensor network approximation of the ground state with finite bond dimension D can be obtained from the exact quantum transfer matrix through a renormalization process where the physical system acts as an impurity. More specifically, the tensor network transfer matrix stemming from a finite D tensor network ground state approximation is a low energy representation of the exact quantum transfer matrix after applying several renormaliza-tion group transformations. Further details of this relation will be published elsewhere.
Appendix A: Derivation of the bound on the decay of momentum-filtered correlation functions This section contains the proof for the claim of Sec. IV D that the correlation length of a momentum-filtered correlation function at momentum k is upper bounded by a constant divided by the minimum of the dispersion around momentum k, plus another constant. The following derivation is general for any spatial dimension and lattice geometry.
Let us first introduce the relevant notations and conventions. Throughout this appendix, we assume to be working on a d-dimensional lattice Λ ⊂ R d generated by the primitive translation vectors a 1 , . . . , a d . The unit cell has a volume V cell = |det[a 1 |a 2 | · · · |a d ]|. Arbitrary lattice sites are denoted as x, y, . . . ∈ Λ; sets of sites are denoted as X, Y, . . . and the cardinality of a set X is denoted as |X|. To every lattice site x ∈ Λ we associate identical Hilbert spaces H x ; the Hilbert space of the whole system is H Λ = x∈Λ H x . This requires that we are working with a finite lattice, and we assume periodic boundary conditions with a period p i ∈ N in the direction of lattice vector a i , such that a site x and x + p i a i are identified for any i (no summation). The lattice is thus given by the set of points Λ = {n 1 a 1 + . . . + n d a d | n i = 0, 1, . . . , p i − 1, ∀i = 1, . . . , d}. We will however be interested in the thermodynamic limit p i → ∞, ∀ i = 1, . . . , d, since this scenario was used throughout the main text.
The reciprocal latticeΛ consists of all vectors K such that exp(iK · x) = 1, ∀x ∈ Λ. In particular, we can define the reciprocal basis vectors b i satisfying b i · a j = 2πδ i,j . The Fourier transform of a lattice function f : Λ → C is defined as and satisfies F (k + K) = F (k) for any K ∈Λ. Hence, we can restrict to momenta k ∈ B, where in the Brillouin zone B is the Wigner-Seitz unit cell ofΛ. Because of the periodic boundary conditions, momentum space is discretized and can be identified withΛ = {n 1 /p 1 b 1 +. . .+n d /p d b d |n 1 = 0, . . . , p i −1, ∀i = 1, . . . , d}. Anticipating the thermodynamic limit and in order to harmonize the notation with the main text, we nevertheless denote the inverse Fourier transformation as with V B = (2π) d /V cell the volume of the Brillouin zone. By using the Euclidean scalar product to define p · x, we can use the Euclidean distance as compatible lattice metric dist(x, y) = x − y ; the distance between two sets X, Y is defined as dist(X, Y ) = min x∈X,y∈Y dist(x, y) and the diameter of a set X is defined as diam(X) = max x,y∈X dist(x, y). We introduce a shift operator T x for all x ∈ Λ that shifts a state |Ψ ∈ H Λ over the lattice vector x. The Hamiltonian is given by H Λ = X⊂Λ H X where the terms H X are supported on a subset X, such that H Λ is translation invariant and local, i.e., there exist positive constants µ, s for which This allows to use Lieb-Robinson bounds [1,2] [A X (t), for two operators A X and B Y supported on disjoint sets X and Y . Furthermore, we assume that H Λ has a unique, translation invariant (i.e. momentum k = 0) ground state |Ψ 0 ∈ H Λ with ground state energy 0. All eigenstates of H with can be labeled by a momentum vector k ∈ B and an index α that labels all eigenstates within a given momentum sector. We denote these energy-momentum eigenstates as |Φ k,α , with eigenenergies E k,α . The lowest excitation energy at momentum k is given by E(k). Note that T x |Φ k = e −ik·x |Φ k for every vector |Φ k in the sector of momentum k. For two operators A X and B Y supported on disjoint finite subsets X, Y , we define the static connected correlation function as We now write A XBY (k) as [Ã X ,B Y (k)] + (A X −Ã X )B Y (k) +B Y (k)Ã X and again use the triangle inequality to bound |C(k)| as where we have defined a new correlatorC(k) = Ψ 0 |[Ã X ,B Y (k)]|Ψ 0 . For both the second and third term on the right hand side of Eq. (A14), we introduce a resolution of the identity, which we separate into two parts, one coming from momentum sectors with momentum k satisfying k − k ≤ δ (for the second term) or k + k ≤ δ (for the third term), and one coming from the rest. For the latter contribution, we use the Cauchy-Schwarz inequality to write for e.g. the third term To bound the new correlatorC(k), we replace the wave packetB Y (k) by a completely local versionB Y (k) defined asB and thus define yet another correlatorC(k) as The error in operator norm can be bounded by the triangle inequality as where an accurate determination of c Λ requires detailed knowledge about the structure of the lattice Λ. We can thus write and using the Cauchy-Schwarz inequality in the different terms and Eqs. (A12),(A13), we obtain |C(k)| ≤ |C(k)| + c Λ (2c erf + 1)e − 2 2r A X B Y .
Finally, to also boundC(k), we separate the time integral in the definition ofÃ X into two pieces. For t < T , we obtain 1 2π where we have also used For the contribution of |t| > T , we use the Gaussian in the integrand to bound The bound onC(k) is thus given by Finally, putting everything together, we obtain In this expression, we can tune the constants r, δ, q, and T as a function of dist(X, Y ) and the characteristics of E(k) around k, i.e. E * (k, δ). It is clear that we want to impose the restrictions: Let us start by setting = α dist(X, Y ) with α < 1 and T = dist(X, Y )/v with v > 2s µ (1−α) . In addition, we choose q = β dist(X, Y ) and r = γ dist(X, Y ) so that We can now fine tune the remaining constant so as to have a similar decay in all exponentials. Assuming the system has reflection invariance, i.e. E * (k, δ) = E * (−k, δ), we can therefore choose the constants such that If there is no reflection invariance, the smaller of both E * (k, δ) and E * (−k, δ) will determine the slowest exponential decay and should be used in the equations above. Clearly, we should set β = 1 vE * (k,δ) . We then obtain µ(1 − α) − 2s v = E * (k,δ) 2v from which we determine the optimal velocity v = 4s+E * (k,δ) µ(1−α) . From the second equality we can fix γ = α/δ. The remaining equation is therefore E * (k,δ) 2v = αδ 2 . Inserting the velocity, we obtain 2µ(1 − α)E * (k, δ) = αδ[4s + E * (k, δ)] which determines α as and only leaves δ to be determined. Note that the restriction α < 1 is satisfied. The fastest exponential decay is obtained by maximizing If we want to minimize the denominator, there is a clear tradeoff since increasing δ decreases the first term and increases the second. If k corresponds to a minimum of the dispersion relation E(k), then the function E * (k, δ) will be insensitive to δ in some region, and we can choose δ as large as possible within this region. However, if k corresponds to a regular point where dE(k) dk (k) = 0, then there is a direct effect from increasing k to decreasing E * (k, δ). A more intuitive result is obtained if we treat the term coming from the Lieb-Robinson bound separately, as the decay properties of this term are specific to the details of the Hamiltonian. Let us assume that we only know about the existence of some maximal velocity of propagation v LR , such that for any v LR |t| ≤ dist(X, Y ), we can write Clearly, choosing v LR larger results in a smaller ξ (a quicker exponential decay of the Lieb-Robinson bound) and vice versa. In order for the Lieb-Robinson bound to be the smallest error in Hastings proof of the exponential decay of correlations, we need to choose v LR at least large enough such that If we now use the new Lieb-Robinson bound of Eq. (A30) in the above, we would again choose = α dist(X, Y ) but we would need to fix T = (1 − α) dist(X, Y )/v LR . Hence, we use a fixed velocity of propagation and do not optimize over it. We will again try to have an equal decay in all exponentials, except for the one coming from the Lieb-Robinson bound, which we allow to decay faster. We thus obtain Setting γ = α δ and β = 1−α vLRE * (k,δ) reduces these equations down to Clearly, the first inequality is trivially satisfied, since E * (k, δ) ≥ ∆E and the fixed velocity v LR satisfies Eq. (A31). From the last equation, we obtain α = 1 1 + vLRδ E * (k,δ) (A34) and the rate of the exponential decay in dist(X, Y ) is given by .

(A35)
We can then optimize over δ to find an optimal decay.

Bounds on Fourier Transforms of Gaussians
We compute the discrete Fourier transform of a sampled Gaussian by inserting the inverse continuous Fourier transform withΛ the reciprocal lattice and V B the volume of the Brillouin zone (which is the unit cell ofΛ). Hence, the Fourier transform of the sampled Gaussian is a sum of Gaussians centered around the different lattice points of the reciprocal latticeΛ. We are only interested in the value of G(k) for k ∈ B, so if r is sufficiently large the contributions of the Gaussians around the points K = 0 will be very small. In general, there exists a constant c gauss so that we can bound G(k) by |G(k)| = N r r 2π