Exact bosonic Matrix Product States (and holography)

We derive an exact formula for a matrix product state (MPS) representation (or a PEPS in higher number of dimensions) of the ground state of translationally invariant bosonic lattice systems in terms of a single one-dimensional Euclidean quantum mechanical path integral with sources. We explicitly evaluate the general formula in the special case of the one-dimensional Klein-Gordon harmonic chain, being a spatial discretization of 1+1 dimensional free boson QFT, obtaining an exact MPS with an infinite dimensional bond space. We analytically diagonalize the transfer matrix obtaining two Fock spaces with continuous modes and check that the exact MPS construction reproduces the correct correlation functions. We also comment on possible holographic interpretations.


Introduction
Matrix Product States (MPS) and more general tensor network constructions arose from the DMRG (Density Matrix Renormalization Group) program and are an effective way of studying wavefunctions, spectra and correlation functions in various interacting lattice systems [1]. Due to the large dimensionality of the Hilbert space of lattice systems, generic wavefunctions would require an exponentially large number of parameters. Tensor networks like MPS, their higher dimensional analogs -Projected Entangled Pair States (PEPS) [2] and further refinements like Multiscale Entanglement Renormalization Ansatz (MERA) [3] provide a very effective variational ansatz with a much smaller number of parameters which captures properties of low lying states of systems described by hamiltonians with local interactions [4]. These ansatzae are build up from tensors with some legs in the physical local Hilbert space at each site and some legs in some auxiliary vector space of dimension D -the so-called bond space. The coefficients of these tensors are found by numerically minimizing the expectation value of the hamiltonian. When one increases D one gets better and better approximations, of course at the cost of increasing computation time and storage requirements.
Apart from these practical applications, a new source of interest in tensor networks appeared coming from the AdS/CFT correspondence [5] or holography [6] -the exact equivalence between nongravitational field theories 1 and theories in a higher dimensional spacetime with boundary which include gravity (predominantly superstring theories). Despite immense progress and overwhelming quantitative evidence, we do not have a microscopic understanding of the emergence of the dual higher dimensional theories directly from the nongravitational boundary field theory.
In [7], Swingle first noticed a remarkable direct analogy between holographic computation of entanglement entropy [8] and a corresponding computation within MERA. Indeed there is an appealing similarity between the structure of MERA and the geometry of Anti-de-Sitter spacetime relevant for holography for conformal field theories. In the case of continuous MERA (cMERA) [9], a proposal was made how to associate a concrete geometry (bulk metric) to the cMERA description of the boundary field theory [10]. Since then, there has been an immense research effort in this direction (see e.g. [11,12,13]). On a cautionary note, however, tensor network descriptions are also quite far from what we would like to have in holography in certain important respects. In particular, we do not have a natural bulk action nor any kind of field theoretical description (which could be potentially interpreted in terms of gravity coupled to a specific matter content).
In contrast to the practical applications of tensor networks as a computational tool for finding wavefunctions and computing correlation functions in concrete strongly interacting systems, when we want to investigate their potential basis for holography, we would like to have an exact description which generically necessitates an infinite dimensional auxiliary bond space. Although variational MPS ansatzae with an infinite dimensional bond space have already appeared in some contexts [14], the present construction has a different origin and aims at obtaining an exact description, albeit for a relatively simple system.
Apart from the case of cMERA, there are very few cases of exact tensor network wavefunctions for natural physical systems. A notable example is the exact MPS for the XY spin chain [15], which was the key source of inspiration for this work. One should also mention a construction of successive approximate MPS for representing CFT correlation functions [16].
The main aim of the present paper is to construct an exact Matrix Product State for a chain of bosons with harmonic interactions (aka Klein-Gordon chain), which is the standard discretization of a 1+1 free (massive or massless) boson. This physically interesting (and trivially exactly solvable) system unfortunately does not fall directly into the class of applicability of so-called Gaussian Matrix Product States introduced in [17] (see comments in section VIIG of that paper), hence this investigation. We hope that the exact Matrix Product State description will be useful for testing hypothesis on the link of holography and tensor networks.
As a byproduct of the construction, we provide expressions (see formulas (23)-(24) below) for the exact MPS (or higher dimensional PEPS) of bosonic lattice systems with a local self-interaction potential in terms of a 1-dimensional Euclidean path integral with sources. These formulas may hopefully be used as an alternative to numerical variational approach for such systems. As the focus of this paper is different, we do not perform any concrete experiments in this direction.
Let us finally comment on the rationale for choosing Matrix Product States and the specific physical system as ingredients for the exact tensor network construction.
On the one hand, we did not want to restrict to a purely conformal system, as scale invariance imposes various very powerful symmetry constraints and a lot of freedom in e.g. modifying the underlying metric as exploited in [12] (but see also [18] for going beyond conformality). Nonconformal systems are, in a sense, much more rigid. On the other hand, the Klein-Gordon chain has of course a gapless limit (and one can get a conformal massless boson in the continuum limit for an appropriate choice of parameters), but we would like to have a common description with gapped systems.
The motivation for treating a discrete system instead of working directly in the continuum despite computational complexities due to discreteness, was that in this context the standard tensor networks would be directly applicable. In addition, one would not have to deal with various infinities and their renormalization. Moreover the exact formulas may be potentially of interest for some condensed matter applications.
In the present paper we choose to work with Matrix Product States instead of MERA, which is typically the tensor network of choice in the context of potential holographic applications. Matrix Product States are conceptually simpler, and in the case of large D (bond dimension) describe well both gapped and gapless systems. For MPS it is easy to implement exact translation invariance even for a discrete system (and finite D). MERA, in contrast, introduces coarse-graining and the very structure of the tensor network does not respect translation invariance by one lattice site.
Last but not least, the structure of a MPS is quite intriguing from the point of view of holography. The boundary points (physical lattice sites) are directly disconnected and interact only through the intermediate auxiliary bond space. As for an exact MPS we expect the bond space to be infinite dimensional, it is tempting to speculate that the infinite dimensional auxiliary space could have a natural interpretation in terms of space of functions in some emergent dimension. The exploration of such a scenario was a direct motivation for this work.
The plan of this paper is as follows. In section 2 we give some further holographic motivations for the computations carried out later in the paper. In section 3 we review the formalism of Matrix Product States, the transfer matrix and the computation of correlation functions. Then in section 4 we give the general construction of exact MPS and PEPS for a wide class of bosonic systems. The following 3 sections are devoted to implementing this construction for the Klein-Gordon chain, diagonalizing the transfer matrix and computing correlation functions. In section 8 we return to the holo-graphic motivations and discuss the obtained results in view of the points raised in section 2. We close the paper with a summary and outlook.
The reader interested just in the exact MPS constructions can safely omit sections 2 and 8. The reader interested in potential links with holography can safely omit sections 6 and 7 and use section 5 just for setting the notation.

Holographic motivation
As described in the introduction, one of the main motivations of the current work is to try to identify some features characteristic of a holographic description in a completely solvable example system for which one can find an exact tensor network description.
A holographic description of a given d-dimensional physical system means that there exists a higher dimensional theory (in at least d+1 dimensions with a boundary) which is equivalent to the original one and thus allows us to compute observables like correlation functions, entanglement entropy etc. using the higher dimensional theory. Moreover, the computational holographic prescriptions typically have an interpretation of putting some objects/operator insertions/boundary conditions at the boundary which defines the particular quantity to compute in the higher dimensional theory.
A quintessential example is the original AdS/CFT correspondence linking N = 4 SYM and superstring theory in AdS 5 × S 5 spacetime. The AdS 5 geometry can be written in the following form where z ≥ 0 is the 'holographic' coordinate, and z = 0 is the boundary. The boundary is intuitively identified with the space where the original field theory lives. Note that the distance between any pair of points at the boundary is infinite, so that all observables are computed through passage into the bulk. Intuitively this is quite similar to the description in terms of Matrix Product States (see next section for an introduction), where there is no direct connection between local Hilbert spaces at different physical sites (i k below), and all interactions occur through the auxiliary bond space (the space of χ k 's) Figure 1: The correlation function of an operator with large dimension can be calculated from the geodesic in the bulk between the two boundary points.
A natural question which arises is whether the auxiliary bond space could have an interpretation of an emergent dimension. The aim of constructing an exact Matrix Product State description in this paper was directly motivated by the above question.
A simple example of a holographic prescription is the computation of a correlation function of an operator of large dimension shown in Fig. 1. According to the AdS/CFT correspondence, we associate to each local operator in the field theory a corresponding bulk field. The two-point correlation function of the field theory operator is roughly given 2 by the boundary limit of the bulk Green's function of the corresponding bulk field [22] O When the bulk field is heavy, one can use the geodesic approximation for the Green's function, and one obtains the picture in Fig. 1.
We will now specialize to the case of a 1+1 dimensional boundary and assume that the bulk field has a free field mode expansion The equal time correlation function would then have the following form with an unconstrained double integral reflecting 3-dimensional bulk modes.
Let us now contrast this with a direct evaluation of a 2-point correlation of a free boundary theory 3 , where the local operator O (x, t) is just the local field φ(x, t) with the mode expansion with µ(p) ∝ 1/ ω(p). The resulting 2-point correlation function is then given by a single integral Of course, for a theory on an integer lattice, the range of p integration would be restricted to the interval p ∈ [−π, π]. The aim of this paper is to construct exact Matrix Product States with an infinite dimensional auxiliary bond space for a bosonic lattice system and check whether some auxiliary modes depending on k appear in a natural way in the computation of correlation functions using the Matrix Product State formalism. In particular we would like to investigate whether the result can be naturally recast into a formula of the form (4).

Matrix Product States
Matrix Product States are a convenient ansatz for representing wavefunctions of 1D quantum mechanical systems which enables efficient computations using the variational method. Their importance stems from the fact that a generic wavefunction of a system with L sites has an exponentially large number of components (dim H) L , where H is the Hilbert space at each site. Even storing this data for large systems is problematic, let alone performing a numerical optimization. Matrix Product States assume instead a factorized structure of the wavefunction in terms of three-legged tensors, given below for a translationally invariant infinite system (see Fig. 2): The i k = 1 . . . dim H are physical indices, while χ k and χ k+1 run from 1 to D, the dimension of an auxillary "bond space". The key ingredient is the 3 legged tensor and it is clear that the number of parameters does not grow exponentially with the system size. As one increases D, one gets better and better approximation of the wavefunction. The approximation is indeed very good for local hamiltonians with a gapped spectrum. However, by increasing D one can also get relatively good descriptions of gapless systems for intermediate distances.
In higher number of dimensions, there is a direct analog of the MPS ansatz, called Projected Entangled Pair States (PEPS) where the basic building block (e.g. in 2D) is a 5 legged tensor The components of A i χχ (or A i χχ χχ ) are typically determined through numerical minimization of the expectation value of the hamiltonian for given values of the bond dimension D. In this work we will concentrate, however, on obtaining exact formulas for the MPS/PEPS representation of bosonic lattice systems. The key point is obtaining a natural identification of the infinite dimensional auxillary bond space. Before we proceed with the construction, let us briefly review the other components of the MPS framework: the transfer matrix and the computation of correlation functions.
In order to normalize the wave function and e.g. to compute correlation functions of local operators, we need to contract the physical index of A i χχ with the corresponding index of its complex conjugate. The resulting composite tensor with 4 legs in bond space plays the role of the transfer matrix: For correct normalizability on an infinite line, we require the existence of an eigenstate |E with eigenvalue equal to unity T F |E = |E (which is also the highest eigenvalue of T F ). This can be represented pictorially as With the above ingredients it is very simple to give a prescription for computing the correlation function of two local operators on a line. It has the following obvious pictorial representation: In order to evaluate this expression one has to diagonalize the transfer matrix T F and compute the form factors of the individual operators, marked in the figure by the dashed lines.
In the rest of the paper we will carry out the construction of all these ingredients. In section 4 we will provide general exact path integral expressions for the MPS and PEPS states (8) and (9). In section 5 we will specialize to the case of the Klein-Gordon harmonic chain and explicitly evaluate the MPS formula, then we will proceed in section 6 to evaluate and diagonalize the transfer matrix (10). As a cross check, we use the above results to evaluate the correlation function (12) in section 7. We then try to look at the obtained formulas for the correlation functions from a higher dimensional perspective in section 8.

General construction of exact bosonic MPS and PEPS
Let us consider a general bosonic lattice system whose Lagrangian can be written in the form or its say 2D version In this section we will allow the self interaction term V (φ) to be arbitrary, but in concrete calculations in the rest of the paper we will specialize to the quadratic case V (φ) = K 2 φ 2 . We would like to construct an exact MPS/PEPS representation of the ground state wavefunction. This is far from trivial even in the quadratic case where we have at our disposal an explicit formula for the exact wavefunction in position representation: where the (doubly infinite) matrix C kl is explicitly known 4 . Note that this form is very far from an MPS ansatz, as in the wavefunction (15) all physical sites are directly coupled to each other. The problem of finding a MPS representation corresponds to decoupling the physical sites through an intermediate auxillary bond space i.e. rewriting (15) as Let us now proceed with the construction in the generic case. The ground state wavefunction up to an overall normalization can be of course obtained in the standard way through an Euclidean path integral with We would now like to rewrite this path integral expression as a product of MPS's. Intuitively, we should partition the Euclidean path integral into strips (or more concretely in the present discrete case into lines ending at the physical sites) and reinterpret these ingredients as successive MPS's. However in order to do that we have to explicitly identify the bond variables. At first glance, it would seem that one could just use the Euclidean u n (τ ) in this role. However it is difficult, if not impossible, to go along with this choice. Firstly, the u n (τ ) are directly linked to the physical variables φ n through the boundary condition u n (T ) = φ n . In contrast, for discrete MPS's A i χχ the bond space indices χ, χ are not constrained in any way by the physical index i. Secondly, due to the direct coupling of u n (τ ) and u n±1 (τ ) in the path integral formula it is difficult to make a consistent choice (valid at each lattice site) which of the three above variables would be identified with the bond variables χ and χ .
In order to overcome these difficulties, we use a less direct construction of the bond variables. As a first step, we will rewrite the spatial interaction term in the action by introducing a functional Lagrange multiplier Now we will interchange the order of integration and first integrate out s n (τ ). Thus the above terms in the action will become 1 2D Plugging this back into the Euclidean path integral we see that we can now directly interpret the expression as a product of MPS's with the bond space being the space of the Lagrange multiplier functions χ n (τ ) defined on a half line τ ∈ [0, +∞). Thus the conventional finite dimensional approximate MPS is substituted by and in the product of MPS's, the summation over the intermediate bond space index becomes a functional integral . .

(22) The final formula for the exact MPS
where we distributed the quadratic terms 1 2D χ 2 n of (20) symmetrically into the two neighbouring MPS's. Clearly the same construction goes through in higher number of dimensions providing an exact representation for PEPS, the only difference being additional χ's in the other spatial directions.
2 +V (u)+i(χ−χ +χ−χ )u+ 1 4D (χ 2 +χ 2 +χ 2 +χ 2 ) (24) Of course even for lattices in higher number of dimensions, the above Euclidean path integral remains just one-dimensional. Note that the formulas (23)-(24) are valid for an arbitrary interaction potential V (u). The above formulas may thus be used as an alternative to the conventional variational approach by picking some finite dimensional basis in the functional space of χ(τ )'s and obtaining matrix elements in that basis by numerically evaluating the 1D Euclidean path integral. Alternatively one could perhaps employ various quantum mechanical approximation schemes for obtaining approximate analytical results. In particular, it should be feasible to apply perturbation theory in their evaluation.
In the present paper we will from now on concentrate on the case of quadratic potential V (u) = 1 2 Ku 2 which corresponds to a spatial discretization of Klein-Gordon theory in 1+1 dimensions. In this case one can evaluate explicitly the MPS (23), as well as diagonalize the resulting transfer matrix.

Exact MPS for the Klein-Gordon harmonic chain
The Klein-Gordon harmonic chain is given by the Lagrangian In the following we will denote ω 0 = K/M . The above system is a standard spatial discretization of the free boson when one takes the lattice spacing ε → 0 with the following identifications and m 0 is the mass in the continuum. However all our subsequent computations will be valid in the discrete case for any values of the parameters M , D and K. We will just comment about what happens in the continuum limit. For completeness let us recall that, in the general discrete case, the excitations of the harmonic chain have energies with p ∈ (−π, π). As a cross check of the MPS construction and the diagonalization of the transfer matrix, in this paper we will compute the correlation functions φ 0 φ m and π 0 π m , where π m = Mφ m . These are given in the discrete case by quite involved expressions 5 which we quote here from [23].
The formula for the MPS (23) can be evaluated using the well known driven harmonic oscillator path integral. In view of diagonalizing the transfer matrix in the following section, it is convenient to keep T in (23) finite as a regularization parameter. The result for and We will fix the overall normalization when we discuss the transfer matrix in the following section. We see that the dependence on φ is Gaussian while the dependence on the arguments of the MPS, i.e. on χ(τ ) and χ (τ ), although also Gaussian is nonlocal. Despite that, one can explicitly diagonalize the transfer matrix, which we will do in the next section.

The transfer matrix
Given the MPS state (32) we construct the transfer matrix by performing the Gaussian integral over φ in whereχ(τ ) andχ (τ ) denote the independent arguments of the lower MPS (with no relation to χ(τ ) and χ (τ ) -see Fig. 3). The * in the above denotes complex conjugation. This integral can be readily carried out giving T 0 χ 2 +χ 2 +χ 2 +χ 2 dτ (36) For ease of notation it is convenient to denote jointly the arguments of the upper and lower MPS entering the transfer matrix as X (τ ) = (χ(τ ),χ(τ )) X (τ ) = (χ (τ ),χ (τ )) (37) We will also use the notation The transfer matrix can then be written as where Diagonalization of the transfer matrix seems to be quite involved at first glance as we have to solve a functional eigenvalue equation involving a functional integral This can be done, however, in two steps. Firstly, let us assume that the quadratic form Σ(τ, τ ) has been diagonalized so that substituting where we introduced a normalization factor which should be chosen so that the MPS wavefunction is normalized to unity. This is translated into the requirement that the highest eigenvalue of the transfer matrix should be equal to 1.
Solving the functional eigenvalue equation (41) now factorizes into a product of 1D eigenvalue problems Each such 1D eigenvalue problem can be solved by harmonic oscillator wavefunctions with eigenvalue where we fixed the normalization factor in (44) to We thus see that we obtain in a natural way an underlying Fock space structure in the auxiliary bond space (which should not be confused with the physical excitations (27)). To complete the picture, we need to determine the eigenvalues Σ l and eigenfunctions f l (τ ) of the kernel Σ(τ, τ ) appearing in (42) and (43).
The Σ(τ, τ ) quadratic form takes the form The 2 × 2 structure can be disentangled by passing to Then the quadratic form becomes 1 2M We thus have two separate eigenvalue problems which we will consider in turn (we will now set M = 1 and reinstate it at the end).
The j + modes Let us first consider In order to solve it we use the fact that ∆(τ, τ ) is a Green's function and satisfies d 2 Then acting with this operator on (53) leads to The fact that ∆(0, τ ) = ∆(T, τ ) = 0 provides, again using (53), the boundary conditions for equation (55) j + (0) = j + (T ) = 0 (56) which leads to the following eigenfunctions and eigenvalues: The j − modes We now turn to the second slightly more complicated problem.
Since K(τ, τ ) is anihilated by d/dτ 2 − ω 2 0 , we will again get with j − (0) = 0 but now supplemented with a more involved boundary condition at τ = T : The first two conditions determine j − (τ ) ∝ sin kτ with Σ − = 1/(ω 2 0 + k 2 ). Then the nontrivial condition (60) simplifies to cos kT = 0. We thus get Transfer matrix eigenvalues Taking the T → ∞ limit, the spectrum becomes continous with k = (l+1/2)π T or k = lπ T , so in both cases we get We thus get two sets of modes parametrized by real positive k which form a Fock space of particles 6 . The eigenvalue of the transfer matrix is given by where n(k) is the integer particle number of mode k. λ(k) is the single particle eigenvalue (given by (46) with n = 1), where we explicitly indicated its dependence on k. It is convenient to write λ(k) in exponential form We will refer to the function E T F (k) introduced here, as the transfer matrix dispersion relation. Its explicit form is obtained by putting together (62), the definition of D in (45) and (46). It reads To summarize, the eigenvalue of the transfer matrix is given by where n(k) is the integer particle number of mode k. It is indeed quite surprising that the dispersion relation of the auxillary Fock space modes is so complicated for the simple bosonic chain. In the continuum limit (26) the dispersion relation (65) however significantly simplifies to

Correlation functions from the exact MPS construction
We will now use the exact MPS and the results on the diagonalization of the transfer matrix to compute the correlation functions φ 0 φ m and π 0 π m . Since their exact expressions are quite complicated in the discrete case for arbitrary separations (see (28)-(30)), this will be an important cross-check of our construction.
The wavefunction to the right of the rightmost operator φ is just the vacuum Since Similarly we have Putting together the above equations as in (12) leads to (79) Removing the regulator by taking the limit T → ∞ gives for the correlation function where E T F (k) is given by (65)-(66), and the form factor of the φ field is given by 7 with x(k) given by (66). It is interesting to contrast the form of the correlation function obtained here (80) with the conventional expression in terms of the physical dispersion relation The modes appearing in the MPS expression should not be confused with the physical modes in the first Brillouin zone p ∈ (−π, π). In addition they give an exponentially damped contribution to the correlation function instead of an oscillatory one as in the expression (83). The equivalence of the MPS expression (80) with the conventional one (83) follows from the considerations in section 8.
The π 0 π m correlation function The calculation of the correlation function π 0 π m goes along the same lines with just some small modifications which we underline here.
Instead of (69) we now have to compute The result is with J + defines as in (71). Its matrix elements are given by The difference w.r.t. the matrix elements of J − comes from the different boundary conditions of the f + l and f − l modes at τ = T . Now we can immediately write the direct analog of (79) and remove the regulator T → ∞ to obtain the final result We checked numerically using Mathematica for a variety of parameters that the obtained expressions (80) and (88) agree with the exact answers given by (28)-(30).

Comments on holographic interpretation
In the preceding sections we have constructed an exact MPS with an infinite dimensional bond space for the bosonic Klein-Gordon harmonic chain. Moreover we have analytically diagonalized the transfer matrix and uncovered a Fock space of two kinds of modes each parametrized by a continuous parameter k ∈ [0, +∞). The eigenvalue of the transfer matrix corresponding to both of these modes 8 is given by In holography, the fields in the dual higher dimensional geometry directly serve to compute correlation functions. This strongly suggests that the potential tensor network counterpart should not really be the MPS but rather the MPS contracted with its complex conjugate. Hence at each lattice position on the boundary we have effectively two Fock spaces with continuous modes parametrized by k ≥ 0.
In the above, the Fock space modes are associated to each lattice site -they are in position space w.r.t lattice. In order to establish contact with the considerations in section 2, it is natural to pass to momentum space w.r.t lattice site positions. Thus the new Fock space modes will be generated by From now on we just consider one set of these Fock spaces 9 . Let us suppose, as in section 2, that a bulk field at t = 0 and z = 0 has the following mode expansion: Then the 2-point equal time correlation function takes the form Let us compare this formula with the expression that we obtained from the MPS computation which we wrote now allowing for an arbitrary sign of n. Since the MPS computation was done in position space w.r.t. the lattice sites, we should try to rewrite the result in momentum space. Indeed we can use the formula π −π dp 2π sinh ε cosh ε − cos p e ipn = e −|n|ε (95) to find F (k, p) so that (93) and (94) exactly coincide. This will occur when Using the quite involved explicit formulas (81) and (89), it is encouraging to find that after many cancellations the final result is very simple and reads (97) 9 More precisely the one with nonvanishing form factor of the φ n field.
Some comments are in order here. On the one hand, we see no trace of a nontrivial background geometry. What we recover is just the flat (discretized along the boundary direction) bulk corresponding the Euclidean path integral (17) as indicated by the form of (97) which looks exactly like the corresponding propagator.
On the other hand, we have a double unconstrained integral which could be interpreted as coming from quantizing some 2+1 dimensional theory and reducing to a spatial timeslice as in (3)-(4) and not a single integral characteristic of a quantized 1+1 dimensional theory.
Most probably these two interpretations are not in contradiction but are rather two complementary points of view on the same expression. In order to reach a definite conclusion, we need, however, to explicitly study time evolution in the exact MPS framework as the picture must be more involved than just (3) in order to account for correlation functions of both φ and π fields at different times and correctly incorporate both sets of modes appearing in the diagonalization of the transfer matrix. We plan to address this issue in future work.
Finally, let us note that using the formulas in this section we can prove analytically the equality of the MPS and standard formulas for the correlation function φ 0 φ n . Indeed, performing first the integral over p in (93) with F 2 (k, p) given by the expression (97) reproduces, by construction, the MPS formula (94), while integrating first over k yields the standard formula (83) thus proving their equivalence.

Summary and outlook
The main results of this paper are the general formulas, valid for any interaction potential, for exact Matrix Product States or PEPS for a bosonic lattice system in terms of a 1-dimensional Euclidean path integral with sources, and their evaluation for the 1-dimensional bosonic harmonic chain. We hope that the general formulas may be of use as an alternative to the conventional variational approach for obtaining MPS and PEPS representations by evaluating the exact formulas for some appropriately choosen subset of functions from the bond space.
A very general feature of the above construction is that the auxiliary bond space becomes the space of functions on a half line. In the concrete case of the Klein-Gordon bosonic chain we can be much more explicit and an-alytically diagonalize the transfer matrix. Its eigenvalues and eigenfunctions form two Fock spaces of modes parametrized by a continuous real parameter k ∈ [0, +∞). The dispersion relation of these modes (understood here as minus the logarithm of the transfer matrix for a single particle state) and given by (65) is surprisingly complicated for such an apparently simple physical system. From the MPS computations of the 2-point equal time correlation functions we obtained quite surprising novel integral representations of these correlation functions which exactly reproduce the well known correlators of this bosonic chain. It would be interesting to employ the MPS construction to more nontrivial position space observables like entanglement entropy etc.
One of the main motivations for this work was that an exact MPS representation of the ground state wavefunction of the Klein-Gordon bosonic chain could be used as a testing ground for various hypothesis on the relation of holography and tensor network constructions. We hope that the result of this paper will be useful in this respect.
From the point of view of computing correlation functions, the most natural counterpart of the holographic (spatial) geometry is not a single MPS but rather the MPS contracted with its complex conjugate keeping the intermediate bond space states arbitrary (i.e. not summed over).
Along these lines, we observed some encouraging and some discouraging results. On the one hand, the infinite dimensional space of modes which diagonalize the transfer matrix indeed can be identified with a Fock space of modes in some emergent dimension. Also, the correlation functions can be written in a way reminiscent of a higher dimensional quantized theory (although this is not the only possible interpretation). On the other hand, we do not see a trace of some nontrivial background geometry. So the relation of the obtained results to holographic expectations is not entirely clear-cut. However, we believe that it is necessary to analyze in detail temporal evolution within this exact MPS framework to make any definite conclusions. This remains as an important problem for future research.
Apart from that, there are numerous other interesting directions for further research, like the development of approximate methods for treating the general formulas (23)-(24) in the interacting case, possible applications of tensor network renormalization [24] in the framework of the present construction, modifying the boundary metric, moving directly to the continuum and merging the above approach with some of the constructions in [25]. We plan to investigate some of these issues in future work.