Continuous Tensor Network States for Quantum Fields

We introduce a new class of states for bosonic quantum fields which extend tensor network states to the continuum and generalize continuous matrix product states (cMPS) to spatial dimensions $d\geq 2$. By construction, they are Euclidean invariant, and are genuine continuum limits of discrete tensor network states. Admitting both a functional integral and an operator representation, they share the important properties of their discrete counterparts: expressiveness, invariance under gauge transformations, simple renormalization group flow, and compact expressions for the $N$-point functions of local observables. While we discuss mostly the continuous tensor network states extending Projected Entangled Pair States (PEPS), we propose a generalization containing a continuum Multi-scale Entanglement Renormalization Ansatz (MERA) and more exotic states


I. INTRODUCTION
Tensor Network States (TNS) provide an efficient parameterization of physically relevant many-body wavefunctions on the lattice [1,2]. Obtained from a contraction of low-rank tensors on so-called virtual indices, they economically approximate the states of systems with local interactions in thermal equilibrium. Their number of parameters scales only polynomially with the lattice size [3,4], circumventing the exponential growth of the Hilbert space dimension. TNS have led to powerful numerical methods to compute the physical properties of complex system [5][6][7], most notably in one spatial dimension d = 1, where Matrix Product States (MPS) [8], the simplest incarnation of TNS, are at the basis of what is arguably the most successful method to describe strongly correlated systems [9][10][11]. In higher dimensions d ≥ 2, accurate results [12] have also been obtained using Projected Entangled-Pair states (PEPS) [13], a natural generalization of MPS. Another family of TNS, Multi-Scale Renormalization Ansatz (MERA) [14], has proved well suited to describe scale invariant states [15,16] appearing in critical phenomena.
Beyond numerical computations, TNS provide important insights into the nature of many-body quantum systems, and have helped describe and classify their physical properties. By design, their entanglement obeys the area law [17][18][19], which is a fundamental property of low energy states of systems with local interactions. They enable a succinct classification of symmetry protected [20][21][22][23] and topological phases of matter [24,25]. TNS also have a built-in bulk-boundary correspondence [26], which makes close connections to physical phenomena appearing in exotic materials [27,28]. Finally, they can be used to build toy models illustrating the holographic principle and the celebrated AdS/CFT correspondence [29][30][31].
For regular spin lattices, PEPS assign a tensor to each lattice site, with 2z virtual and one physical (spin) indices, * antoine.tilloy@mpq.mpg.de † ignacio.cirac@mpq.mpg.de where z is the coordination number. The virtual indices are contracted according to the lattice geometry, yielding a wavefunction for the spin degrees of freedom. This description is particularly useful in translationally invariant systems, as this symmetry may simply be imposed by choosing the same tensor on each site. For MERA [14,32], a tree-like structure of two types of tensors is used. In both cases, the whole many-body wavefunction is determined by one or few tensors, which encode all the physical properties. An important challenge in the theory of TNS is the generalization from lattice to continuous systems. Such an extension would allow the direct study of Quantum Field Theories, without the need for a prior breaking of spatial symmetries with a discretization. Further, the continuum provides a whole range of exact and approximate analytic techniques (such as exact Gaussian functional integrals, saddle-point approximations, or diagrammatic expansions) that have no obvious discrete counterparts and that could provide useful additions to the TNS toolbox.
A natural way to carry out such a program is to simply take the continuum limit of a TNS, by letting the lattice spacing tend to zero while appropriately rescaling the tensors. In fact, this has been done in one spatial dimension, d = 1, where it yields continuous matrix product states (cMPS) [33,34]. In higher dimensions, however, the task does not seem trivial. Naive extensions of cMPS have a preferred spatial direction and break Euclidean symmetries [35]. In [35], a proposal for cPEPS was put forward to overcome such a limitation, but the resulting state was no longer obtained from the continuum limit of a TNS. Thus, so far there seems to be no fully satisfactory way of extending TNS to the continuum in d ≥ 2.
In this article we propose a definition of continuous tensor network states (cTNS) that naturally extends TNS to the continuum. We obtain them as a genuine continuum limit of TNS, but manage to preserve Euclidean invariance. As in previous works [35][36][37], we exploit the similarity between a tensor contraction over the indices lying on the links of a tensor network and a functional integral over a field living on the continuum limit of this mesh. The key difference lies in the way the continuum limit is taken in higher dimensions: As we shall argue, the d = 1 case of cMPS is too peculiar to be directly extended.
The first definition of cTNS we will propose in section II takes the form of a functional integral over auxiliary scalar fields as advertised. From this definition, which makes local Euclidean invariance manifest, we will derive an operator representation similar to the one used for cMPS. Importantly, we will show in section III how this ansatz can be obtained from a continuum limit of a discrete TNS. We will then study some of its properties reminiscent of the discrete: its ability to approximate (possibly inefficiently) all states (IV A), its redundancy under some so called Gauge transformations (IV B), which play a crucial role for PEPS, its flow under scaling transformations (IV C), and its cMPS approximation in some carefully chosen limit (IV D). We will then propose various methods to carry computations with Gaussian and non-Gaussian cTNS (V). While most of our approach is aimed at the continuum limit of PEPS, we will finally generalize it to MERA-like states (and more exotic TNS) in arbitrary dimensions by including a metric and restricting physical fields to a boundary (VI).

II. CONTINUOUS TENSOR NETWORK STATES
We start by giving two equivalent definitions of continuous tensor network states, leveraging a functional integral and an operator representation. Our objective at this stage is only to provide a definition of a class of states for bosonic quantum fields, with only a crude intuition for why such an object could indeed be a good definition of a cTNS. We forgo the derivation of this cTNS from a class of discrete tensor networks to the following section.

State definition
We begin with the functional integral representation. It will be the most direct to derive from the discrete and makes Euclidean symmetries manifest. Definition 1 (Functional integral formulation). A continuous tensor network state (cTNS) of a bosonic quantum field on a domain Ω ⊂ R d with boundary ∂Ω, is a state |V, B, α parameterized by 2 functions V and α: R D → C, and a boundary functional B: L 2 (∂Ω) → C defined by the functional integral on an auxiliary D-component field φ: where |0 is the physical Fock vacuum state, The functions α and V may depend explicitly on position.
The auxiliary D-component field φ, which is integrated over in the functional integral, is the continuous equivalent of the auxiliary bond indices that are contracted in tensor network states (see Fig. 1). This intuition will be made more precise in the next section. For this reason, we call D the bond field dimension. As we shall see in section IV, the bond field dimension D bears similarities with the bond dimension χ of discrete tensor network states.
If Ω is R d or a torus with periodic boundary conditions, B can simply be set to 1. In that case, the state and its associated properties depend only on V and α, in the same way a TNS depends only on local tensors. If V and α do not depend explicitly on x, the cTNS describes a translationally invariant state. More generally, if ∂Ω = ∅, the boundary functional could induce e.g. : 1. Dirichlet boundary conditions: B(φ| ∂Ω ) ∼ δ(φ| ∂Ω ) fixing φ| ∂Ω = 0 in the functional integral, 2. Neumann boundary conditions: B(φ| ∂Ω ) ∼ δ(∇φ·n) fixing ∇φ · n| ∂Ω = 0 where n is normal to ∂Ω, 3. Something more general, given e.g. by a quasi local functional: where L is a function from R (d+1)D to C. This latter option will be generated naturally when we discuss gauge invariance in IV B.
We may rewrite expression (1) more explicitly as a sum over unnormalized field coherent states. Introducing the massless free field probability measure dµ(φ) for the auxiliary field: and a complex amplitude A V (φ): yields: where is an unnormalized field coherent state. Hence, just like cMPS in dimension 1, cTNS are a generalization of field coherent states. The latter are obtained e.g. if dµ(φ) is only nonzero for a given φ (for an infinitely deep V ), or, in the homogeneous case, if α is constant.

N-particle Wave function
A generic state |Ψ in the bosonic Fock space F[L 2 (R d , C)] can be expanded into a sum of n particle wave functions ϕ n : (6) where ϕ n is a completely symmetric function of its coordinates. Simply expanding the exponential of eq. (1) gives, for the cTNS |V, B, α : It provides an equivalent definition of the cTNS.

Correlation functions
A state is also fully characterized by its (equal time) correlation functions. To compute them, we first introduce the generating functionals for real sources j , j: They generate the normal ordered and anti normal ordered correlation functions respectively. For example, it is straightforward to verify that: Using the formula for the overlap of (unnormalized) field coherent states, and writing N = V, B, α|V, B, α we get: We observe an important fact which is that the function α appears squared, hence an α quadratic in the field already brings non-Gaussianities. To compute Z j ,j , one applies the Baker-Campbell-Hausdorff (BCH) formula to (8) to push annihilation operators to the right and get back to a computation of field coherent state overlaps. We obtain: hence the same as (12) but for the removal of the product j · j , which is responsible for the divergent equal point contributions upon functional differentiation.

State definition
We now provide an equivalent operator representation of cTNS. For simplicity we restrict ourselves to domains of R d that can be split into Cartesian products Ω = [−T /2, T /2] × S where ∂S = ∅. This is not strictly necessary but substantially simplifies the definition. [−T /2, T /2] and x ∈ S. A cTNS is then defined as: where T is the τ -ordering operator,φ k (x) andπ k (x) are k independent canonically conjugated pairs of (auxiliary) field operators: i.e. D copies of a bosonic Fock space on a d − 1 dimensional space. The trace is taken over this auxiliary Hilbert space. As before, V and α may depend on x and τ .
The operatorB acts on H aux and fixes the boundary conditions, e.g.B = 1 encodes periodic boundary conditions on the coordinate τ . Another natural option is to takeB = |in out|, which corresponds to the situation of Fig. 2.
Definitions 1 and 2 are equivalent for this subclass of domains. The proof is straightforward. One just applies the techniques of standard QFT textbooks to go from operator to functional integral representations with τ = it (see e.g. [38,39]). Mainly, one discretizes the τ -ordered product in (14) into a finite product of terms. One then inserts resolutions of the identity in the field basis |φ at every time step ∆τ and writes each resulting overlap in the conjugate momentum basis |π . Going back to the continuum limit yields a phase space functional integral which reduces to the formula of equation (1) upon Gaussian integration of the conjugate momenta π.
The boundary operator of (14) is related to the boundary functional of (1) by: where |φ out and |φ in are eigenstates of the auxiliary field operatorsφ(x). The auxiliary field φ decomposes into φ = φ in + φ out where φ in (resp. φ out ) has support on ∂Ω in (resp. ∂Ω out ) with ∂Ω = ∂Ω in ∪ ∂Ω out .
As before we may propose a repackaging of formula (14). Introducing the Hamiltonian density H( yields: This is a straightforward extension of the cMPS definition [33]

N-particle wave function
The N-particle wave-function ϕ n defined in (6) can also be computed in the operator representation. For −T /2 < τ 1 < · · · < τ n < T /2, we get, expanding the τ -ordered exponential (14) into an infinite product: As for cMPS [33], we may interpretĜ as a propagator andα as a scattering matrix creating a particle. It is the very specific form of H and hence ofĜ that is responsible for the Euclidean symmetries of the resulting state. Generalizing cMPS starting directly from (17) would make it hard to guess an appropriate expression forĜ.
Note that (17) amounts to taking as wave function a correlation function of auxiliary quantum fields. This is similar in spirit with the Moore-Read states [40] used for Hall physics or with the infinite matrix product states [41] used for critical spin chains.

Correlation functions
Finally, we may provide an expression for the generating functionals Z j ,j Z j ,j in the operator representation. Exploiting the operator definition of the cTNS (14), expanding the τ -ordered exponential into an infinite product of infinitesimal exponentials, we get: with the transfer matrix (with sources): Using as before the BCH formula yields: The functional derivatives can then be carried explicitly and one obtains e.g. , for −T /2 < τ 2 < τ 1 < T /2: and the transfer matrix T := T 00 . More generally, correlation functions are given by the trace of a succession of propagators M followed by operator insertions of α ⊗ 1 (respectively 1 ⊗ α * ) in the positions corresponding to ψ (respectively ψ † ).

A. (Discrete) tensor network states
We start with a very brief reminder on tensor network states (TNS), recalling only their elementary definition. For an understanding of their efficiency in representing quantum systems of physical interest, we direct the reader to the relevant literature (e.g. [42,43] and references therein).
TNS are variational ansatz for many-body wave functions that take the form of a contraction of local tensors. The simplest example, in spatial dimension d = 1, is provided by matrix product states (MPS). For a translation invariant quantum spin 1/2 chain with N sites, a generic state reads: where the wave function c i1,··· ,i N contains 2 N complex parameters. A MPS is a economical ansatz for this wave function: where A −1 and A 1 are two χ×χ matrices. These matrices contain the parameters which allow to vary the state. Their size χ, called the bond-dimension, encodes the depth of the variational class and upper bounds the amount of spatial entanglement that can be carried by the state. The two matrices can be be collected into a 3-index written graphically: where i is the physical index and k, are so called bond indices. Graphically, the corresponding wave function c can be written: where joint legs of the tensor A denote a summation on the corresponding index, associated with the matrix multiplication and subsequent trace in (23). This representation makes it natural to generalize matrix product states to projected entangled pair states (PEPS) [13] in arbitrary dimensions, e.g. in d = 2 for N × N sites: This latter object, in general d, is what the cTNS defined in (1) or (14) aims to extend to the continuum.

B. Constructing cTNS
Our objective is to show how cTNS can be obtained from a limit of a discrete TNS. To motivate this discrete ansatz, we first provide heuristics for why its main characteristics seem to be required. Mainly, we aim to justify: 1. why infinite bond dimension is needed, and 2. why the trivial tensor around which we expand is of the form we postulate.
The first point is a scaling argument. We ask for a strong notion of continuum limit: we require the discrete tensor to be approximately stable by fine graining to the UV. Namely, the discrete ansatz needs to be (at least approximately) expressible as a contraction of tensors with the same form but different parameters. Each blocking multiplies the physical dimension by 2 d , and this is why in the continuum limit one obtains a field theory on the physical degrees of freedom. But each blocking also multiplies the bond dimension by 2 d−1 (see Fig. 3). Hence, for d > 1, the bond dimension is increased when zooming out and decreased when zooming in. The only way to make the class of states considered approximately stable is for the bond dimension to be infinite. Notice in this respect that the d = 1 case allows finite bond dimensions even in the continuum limit [33,34]. It will be important to see if this peculiarity can be recovered in some appropriate limit in section IV D. Note that our argument in favor of infinite bond dimension does not imply that a discrete tensor network state with finite bond dimension could not behave, at distances sufficiently large compared to the lattice spacing, like a cTNS. Rather, our argument shows that any simple space discretization of a cTNS in d ≥ 2 into a TNS will have infinite bond dimension, even for arbitrarily small lattice spacing. As in [35], it could also be that a proper choice of boundary conditions would constrain the tensor contraction on a finite dimensional subspace, despite an apparent infinite bond dimension.
We now need to discuss more precisely the form of the elementary tensor. For notational simplicity, we now assume that an elementary tensorT to be contracted: is a vector in its bond indices but an operator acting on the vacuum in the physical space, namely: To obtain a continuum limit, the crucial choice lies in the elementary "trivial tensor", acting as the identity on the vacuum, and around which to expand: Indeed, it is natural to want the tensor corresponding to zero particle in an elementary cell of the physical space to dominate. It seems that any other choice would preclude the existence of a continuum limit. In d = 1 dimension, there is only one natural option which is to take the tensor corresponding to the identity: But in the same way as for the bond dimension, the situation is a little too trivial in d = 1 to give a precise hint for higher dimensions. In d > 1, there are several seemingly natural options which we have to inspect. We will discuss the d = 2 case but the reasoning holds for any d ≥ 2. We do not aim to prove that the tensor we will ultimately expand around is the only option, but rather that other seemingly simpler options bring difficulties. 1. A naive option is to generalize the identity on the auxiliary bond space by takingT (0) ijkl = δ ijkl , that is to take an elementary tensor corresponding to a Greenberger-Horne-Zeilinger (GHZ) state on the bond indices:T As will later be manifest, this choice is too brutal and would yield a state with a trivial spatial structure.
2. Another simple option is to take the identity along a diagonal, e.g. :T This is the choice that stays the closest in spirit with the d = 1 case. The problem of such a choice (made e.g. in [35]) is that it picks a prefered a direction and thus makes Euclidean invariance impossible to obtain directly (that is, without analytic continuation).
3. We may combine the 2 d−1 identity operator along diagonals in a sum, as an attempt to recover the Euclidean invariance lost with the previous choice: An issue is then that the corresponding tensor contraction contains loops: (35) which yield divergent terms as tr [1] = +∞ for infinite bond dimension.
None of the natural options seems to provide a simple Euclidean invariant continuum limit. Our proposal will consist in taking a regularized version of the first possibility, in the form of a "soft" delta: C. Discrete Ansatz in d = 2 The first lesson from the previous sections is that infinite bond dimensions seem to be required. We thus write the bond indices of the elementary tensor as D real numbers. In d = 2, this means that an elementary tensor has 4 bond indices φ(1), φ(2), φ (3), and φ(4) ∈ R D : As we mentioned before, the heart of the problem of the continuum limit lies in defining the proper trivial tensor around which to expand. We choose a "soft" delta: This ansatz forces the bond indices to remain close to each other, and contributes to the generation of the gradient squared term in the action. Being Gaussian, we also naturally expect its form to be stable. To this "trivial" part, we add local corrections of order ε d = ε 2 : In this expression, φ denotes whatever combination of the bond field indices φ(1), φ(2), φ (3), and φ (4). The simplest possibility is to take φ as the average of the bond indices but it does not matter for the continuum limit. The operator ψ † (x) anticipates the continuum and has commutation relations [ψ(x), ψ † (y)] = 1 ε 2 δ x,y δ 2 (x − y). Ignoring the boundary conditions for now, the contraction of the tensors amounts to integrate over all the bond indices: and we see that the differences in eq. (39) yield: We recognize the (rotation invariant) gradient square term of the continuum definition 1. Defining the path integral "measure" as: finaly yields the continuous tensor network state of equation (1) up to boundary conditions. To get a state on the physical Hilbert space, the auxiliary fields on the boundary just have to be contracted (or integrated) against a boundary functional which we wrote B in (1).

D. Discrete Ansatz in general d
For d ≥ 3 the derivation is carried along the same way as before. We just note that there is a small peculiarity in the d = 2 case because the auxiliary fields are adimensional. To generalize equation (38) to higher spatial dimensions d > 2, one naturally extends the prescription of summing all the differences of the squares of the nearest bond indices φ(1), · · · φ(2d). But, importantly, to obtain the continuum limit, one needs to multiply this expression by ε d−2 where ε is the length of the unit cell.
to obtain the integral of a gradient squared in the continuum limit. In d = 2, the ε d of the integration measure and ε −2 from the gradient square cancel each other and this scaling factor does not appear.
In retrospect, it is clear why deriving the continuum limit by perturbing around the GHZ tensor (32) would have given a trivial continuum. It would have corresponded to putting an infinitely large constant instead of ε d−2 in (43), an infinite "rigidity" that could not be compensated by locally small terms.

IV. PROPERTIES
We now explore the properties of cTNS that are analogous to those of their discrete counterparts.

A. Stability and expressiveness
It is first natural to wonder how "big" the class of cTNS is. It could be, for example, that even for arbitrary large D and arbitrary V , B and α, cTNS only spanned a small sector of the Fock space. As in the discrete, can any state be approximated, even if inefficiently, by a cTNS?
Let us first consider the stability of the cTNS class. The sum of two cTNS is still a cTNS, provided we are willing to accept singular potentials V . More precisely, let |V 1 , B 1 , α 1 and |V 2 , B 2 , α 2 be two cTNS with bond-field dimension D 1 and D 2 . Then we can easily rewrite their sum as a cTNS with bond-field dimension D 1 + D 2 + 1 (although it may in general require fewer auxiliary fields). For example, defining the cTNS |W Λ , C, β with: where θ is the Heaviside function, we indeed have |W ∞ , C, β ∝ |V 1 , B 1 , α 1 + |V 2 , B 2 , α 2 . Indeed, when Λ is sent to infinity, the auxiliary field φ becomes a "bit" taking values ±1 digitally splitting the functional integral into two contributions D φ φ≡±1 where each term of the sum gives the two initial states.
The expressiveness of cTNS is then easy to assess, following the same technique as for cMPS [35]. Taking α(x, φ(x)) = f (x) and V = a/Vol(Ω) we obtain any field coherent state with any complex weight e −a |f . Using the stability result, one can construct arbitrary linear combinations of such field coherent states which are dense in Fock space, hence one can get arbitrarily close to any state in the Fock space. With this construction, the bond field-dimension grows at each addition of coherent states.
Actually, using larger bond field dimensions is only a convenience and, provided V and α are arbitrary, a cTNS can approximate any state in the Fock space with D = 1. Let us consider a sum of coherent states |Ψ = m j=1 e −ai |f i . This sum can be approximated by the cTNS |V Λ , 1, α with: where 1 A [x] = 1 if x ∈ A and 0 otherwise. Indeed, when Λ → +∞, the auxiliary field is forced to sit on one of the m minima of the potential, which each have a complex weight e −aj . To each of these m possible values of the field, the term in α associates a different coherent state. Hence one can approximate |Ψ with arbitrary precision and hence all states in the Fock space.
Allowing larger values of D remains useful if V and α are restricted in some way, e.g. to being polynomials with a fixed degree. In that case, being able to take a larger bond field dimensions D substantially increases the expressiveness of a cTNS subclass. Gaussian cTNS (see V A) will provide such an illustration.

B. Gauge transformation
Different choices of V , B, and α can generate the same state. This is to be expected: In the discrete, the map between an elementary tensor and a many-body wavefunction is not injective either. Understanding the transformations between tensors generating the same state is fundamental in the theory of TNS, especially for the classification of symmetry protected and topological phases. It is thus natural to ask the same question for cTNS following the discrete construction.
where = and = gives the same contracted state up to new boundary terms (that vanish on a torus). Such transformations have proved central to classify topological phases of matter with discrete tensor networks. We would thus like to find an analog in the continuum.
For infinite bond dimension, the equivalent of an invertible linear transformation acting on discrete indices is a linear operator G acting on functions of D real variables (the auxiliary field): A G that is too generic will typically destroy the continuum limit when the corresponding gauge transformation (49) is applied on the elementary tensor. The main difficulty is to know what subset of operators to look at. For simplicity, we restrict ourselves to the case D = 1 from which the general case is easily deduced. A first option is to consider the subset of diagonal transformations. Let F and G be two operators acting diagonally: Acting on an elementary discrete tensor in the same way as in (49) with F and G simply changes the integration measure: dφ(1) · · · dφ(4) → f (φ(1))g(φ(4)) f (φ(3))g(φ(2)) dφ(1) · · · dφ(4) (53) If this change of measure is too general, there will be no continuum limit. A natural choice, preserving the continuum, is to take: In the continuum limit such a choice put in (53) yields: hence this adds a pure divergence term into the cTNS definition which can be transformed into a boundary term thanks to Stokes' theorem. This is exactly what a gauge transformation should do.  (49) has a non trivial result on the boundary only. In the continuum (right), the transformation of the elementary tensor is equivalent to the addition of a pure divergence term for the auxiliary fields in the bulk, which can then be integrated into a boundary condition.
This very special choice of operators does not exhaust the infinitesimal transformations compatible with the existence of a continuum limit. However, we conjecture that all discrete gauge transformations of the form (49) that preserve the continuum limit ultimately give rise to pure divergence terms as well. In any case, this discussion of the discrete setting is but a motivation for the introduction of (some) continuous gauge transformations of cTNS.

Continuum description
The previous inquiries motivate the following proposition which is at the same time a definition of a certain class of gauge transformations for cTNS.

Proposition 1 (Gauge transformation). Let F [x, φ(x), ∇φ(x)] be an arbitrary vector field in Ω.
If Ω has no boundary, the cTNS |V, α is left unchanged by the gauge transformation: The proof is trivial and is just a direct application of Stokes' theorem. More generally, if Ω has a boundary ∂Ω, the gauge transformation (57) adds a boundary term to the measure: where n(x) is the unit vector normal to ∂Ω in x. Gauge transformations of cTNS thus have a straightforward geometric interpretation.

C. Tensor rescaling
Our objective is now to relate different tensor network descriptions of the same state at different scales [14,44]. More precisely, considering a correlation function for a state parameterized by a tensor T (1) in the thermodynamic limit: C(x 1 , · · · , x n ) = T (1)|O(x 1 ) · · · O(x n )|T (1) , (59) the objective is to find a tensor T (λ) of new parameters such that: Naturally, in the discrete, this relation is at best approximate.
For the cTNS of definition 1, we can write the flow V (λ), α(λ) exactly, following the rather standard dimensional analysis of ordinary QFT. As such, we introduce new creation and annihilation operators ψ † (x) = λ d/2 ψ † (xλ) and ψ(x) = λ d/2 ψ(xλ). They indeed verify the standard commutation relations [ ψ(x), ψ † (x)] = δ(x − y). These new operators relate correlation functions at different scales: where the λ nd/2 factor just comes from the fact that the operators we introduced have a dimension. We just have to rewrite |V, α as a function of the new creation and annihilation operators to relate the different scales. To achieve this, we change of position variable introducing u = x/λ. The free field measure now reads: with φ(u) = λ d−2 2 φ(λu). This gives: This allows to discuss the IR behavior of cTNS in terms of relevant, irrelevant, and marginal couplings. To this end, we informally expand V and α in powers p of the fields φ and analyze the terms of each degree separately. The corresponding coupling dimensionality is ∆ = d + 2−d 2 p for terms in V and ∆ = d 2 + 2−d 2 p for α. Consequently: -For d = 2, All powers of the field in V and α yield relevant couplings. There are no irrelevant couplings.
-For d = 3, the powers p = 1, 2, 3, 4, 5 of the field in V yield relevant ∆ > 0 couplings. The power p = 6 is marginal in V . For α, the powers p = 1, 2 are relevant and p = 3 is marginal. All other powers are irrelevant.
Relevant powers will dominate the behavior of cTNS correlation functions in the IR and actually be the only ones allowed if the cTNS description is aimed to hold at all scales in the non Gaussian case (see V B).

D. Recovering continuous matrix product states
Continuous TNS should reduce to cMPS in an appropriate limit. This is an important property to check to demonstrate that our ansatz is a natural extension of cMPS to d ≥ 2.

Compactification
A cMPS of a quantum field defined on a space interval of length T is parametrized by 3 (χ × χ) matricesQ,B,R and is defined [33,34] as: |0 .
(68) To obtain a cMPS, we may directly instantiate our cTNS ansatz with d = 1, e.g. using its functional integral form (1). However, as we mentioned before, the d = 1 case is quite peculiar compared to other dimensions and so it is nice to see it can also be immediately obtained from a general d case where all dimensions but one are taken to be very small.
Indeed, consider a domain of the form Ω = [−T /2, T /2] × S where S is a d − 1 dimensional torus of length in all d − 1 directions. Expanding the auxiliary fields φ of (1) in Fourier modes on this d − 1 torus S and taking the limit → 0, yields a functional integral in which only the field zero mode on S survives. Hence, one obtains a functional integral of the form: where the 1 dimensional auxiliary field φ is the zero mode (on the shrunk torus S) of the initial d dimensional auxiliary field. It is what one would have obtained fixing immediately d = 1 in (1). In operator form, this yields: whereP k andX k are canonically conjugated pairs (D zero dimensional quantum fields). This is already a cMPS . However, the bond Hilbert space is now that of D particles in 1 dimension or 1 particle in D dimensions, hence χ = +∞.

Bond dimension quantization
To obtain a genuine cMPS (with finite bond dimension) from a d = 1 cTNS defined by (70), we need to choose a specific potential effectively reducing the Hilbert space dimensionality. The intuition is quite clear: take a potential with deep minima.
Let us take a potential with D deep minima m k on the vertices closest to 0 of an hypercube, i.e. m 1 = (1, 0, · · · , 0), m 2 = (0, 1, 0, · · · , 0), ..., and m D = (0, · · · , 0, 1). The effective dynamics is now restricted to a D-dimensional Hilbert space spanned by |m k corresponding to wave packets localized around each minima. In this reduced Hilbert space, the minima are coupled by tunneling. Because of the geometrical configuration we have considered, the minima can all be connected by independent saddle points, hence the effective coupling between the minima can be chosen freely. This means that we can obtain any D × D complex matrix Q of standard cMPS by adjusting the value of the D 2 saddle points of V .
The R matrix is fixed in the same way. The values of α[X] on the minima |m k of the potential fix the diagonal coefficients of R and the value on the saddle point connecting |m k1 and |m k2 fix the non diagonal terms R k1,k2 .
Hence, not only can cTNS reduce to cMPS when d = 1 (or when d−1 dimensions are small) for a specific choice of potential, but actually all (bosonic) cMPS can be obtained this way. In this context, the bond field dimension D reduces to the usual bond dimension χ.

V. COMPUTATIONS
To carry computations with cTNS, one could of course rediscretize them and use the standard TNS algorithms. We now mention techniques relying only on the continuum limit.

A. Gaussian states
There exists a subclass of cTNS for which all quantities of interest can be computed exactly: Gaussian cTNS.

Definition 3 (Gaussian cTNS). A cTNS is said to be
Gaussian if the functions V and α are respectively at most quadratic and affine in the auxiliary field: Naturally, a Gaussian cTNS is also a Gaussian state in the usual sense of the term. More precisely, for a Gaussian cTNS, Z j ,j and Z j ,j are manifestly Gaussian functionals. Let us compute Z j ,j in the translation invariant case Ω = R d .
Inserting definition 3 into equation (13) yields: Carrying the Gaussian integration we then obtain: where: Because of translation invariance, K(x, y) = K(x − y) which can be written in Fourier space: Inserting this expression into equation (75) and integrating over the variable u = (x − y) yields: which is difficult to make more explicit but could be computed exactly for given α and V . To get a intuition of the behavior of the two point functions, we may instantiate this expression on a simple example where the bond field-dimension D equals 1 and α (0) = V (1) = 0 for simplicity. In that case we have: Using (9) this gives the correlation function: Importantly here, the correlation function in momentum space C(p) ∝ p −4 when p → +∞. Hence, the integral is not UV divergent for x = y so long as d ≤ 3 in which case the particle density ψ † (x)ψ(x) is finite.

B. Non-Gaussian states
For a non-Gaussian cTNS, it is no longer possible to compute the correlation functions exactly in general. Fur-ther, the definitions we provided for the cTNS in (1) or (14) are generically divergent. Nonetheless one can use approximations or numerical techniques coming from the quantum field theory and tensor network toolboxes.

Regularization and renormalization
In the general case, the ansatz we put forward suffers from the same UV divergences that plague quantum field theories. As in QFT, these divergences are in a way inevitable: the gradient squared (∇φ) 2 in the path integral insufficiently penalizes high momenta in d ≥ 1 (note that, again, the d = 1 case is trivial). On the other hand, the locality of the underlying tensor network forbids higher derivatives. Hence, as in QFT, the divergences are tied to the very property (locality) that we require.
Given this state of affairs, there are essentially 3 options to deal with divergences, depending on what one needs the state for.
The first option is simply to regularize the state with a momentum cutoff Λ, either directly in the path integral -which will break locality and destroy the operator representation-or in the operator representation -which will generically break Euclidean invariance-. In both cases, the scale Λ will be reminiscent of the inverse lattice spacing of discrete tensor networks. The parameters appearing in the expansion of V and α will then be the equivalent of the bare parameters in QFT Lagrangians. As long as the state is used as a variational ansatz e.g. to minimize the energy of an anyway regularized QFT Hamiltonian, this is unproblematic. Indeed, carrying an optimization on bare or renormalized parameters will be equivalent, and the fact that some properties break above a cutoff momentum is anyhow imposed by the physical QFT being approximated. In this approach, there is no restriction on the powers of the auxiliary field appearing in V and α.
One may also be interested in the class of cTNS for their properties, and not necessarily to approximate the ground state of a given system. In that case, going beyond regularization and renormalizing the state with proper counterterms and renormalization conditions seems necessary to preserve the locality of the underlying tensor network. In the general case, this is equivalent to renormalizing a relativistic open quantum field theory, a problem which has received interest recently [45]. At the level of dimensional analysis, this restricts powers of the auxiliary field in V and α to renormalizable interactions, hence to the relevant and marginally relevant powers obtained in IV C. In d = 3, this restricts the parameters to a finite number of tensors appearing in the finite polynomial expansion of V and α. Allowing for more auxiliary fields is thus necessary to make the cTNS class arbitrarily large and expressive in d = 3.
Finally, a natural regularization may be provided by restricting the class of quantum states in d − 1 on which the transfer matrix T acts. As we will see, for special cases of T, one can indeed recover finite results in d = 2.
In that case, the state itself is implicitly defined by the approximate method used to contract it.

Dimensional reduction
We now discuss the last option. To compute physical correlation functions in the general case, one can exploit their operator expression (20) given by the exponential of a transfer matrix acting in a space of one dimension less. In the d = 2 case, the theory one needs to solve is thus simply a 1 dimensional QFT. The latter is solvable with cMPS (i.e. cTNS in 1 dimension less) which, as a bonus, have a built-in UV regulator [46] and bring the computation back to a 0 dimensional problem [47]. We outline the steps of such a computation on a simple example.
We consider a cTNS on a torus (τ, x), τ ∈ [0, T ], x ∈ [0, L] (as in Fig. 2), which reads, in the operator representation: Correlation functions for this state at a fixed τ take a particularly simple form, with the propagator exp(τ T) appearing only once. For example, the 2-point function reads: Such N -point functions at equal τ contain useful information about the state in the thermodynamic limit T → +∞, L → +∞. Indeed, because of Euclidean invariance, they give access to all correlation functions of aligned points, and a fortiori to all possible 2-point functions. This is sufficient to compute the expectation values of most homogeneous and isotropic quasi-local Hamiltonians.
To simplify the discussion, we now consider the special case of an Hermitian T (obtained e.g. when all the coefficients of V and α are real). In the T → +∞ limit, exp(T T) will be dominated by the projector on the eigenvector |ss of T with the largest eigenvalue. The correlation functions then simplify, e.g. : The right hand side is the correlation function for (2D copies of) a 1 dimensional bosonic field theory, which motivates the use of a cMPS. More precisely, we may use a cMPS defined on two copies of the auxiliary quantum fields to approximate the dominant eigenvector |ss .
Assuming only D = 1 auxiliary field, we can write: where Q, R 1 , and R 2 are (χ × χ) matrices and 1] are the creation operators associated to each copy of the Fock space on which T acts, |0 is the Fock vacuum of these two copies, and the trace is taken over the matrices. This is nothing but a translation invariant cMPS for two species of bosons. The dominant eigenvector can then be approximated by choosing The right hand side of (85) can be computed explicitly as a function of Q, R 1 , R 2 . Indeed, in the same way as we computed the normal ordered correlation functions for a cTNS in (20), one can compute the normal ordered correlation functions for a cMPS, replacing −H(x) by Q and α(x) by R 1 , R 2 [34], e.g. for x ≥ y: with the (0-dimensional) transfer matrix: One then just has to express T as a function of ψ 1 and ψ 2 instead of the field and conjugate momenta, which requires a choice, e.g. : for some Λ 0 . Taking the expectation value of products of local operators on the cMPS yields divergent contributions. They can be removed e.g. by normal ordering H and α [48] in the operator representation of (14), or by adding a counter term in the Hamiltonian as in [49]. In the end, the expectation value to maximize can be written: where M (Q, R 1 , R 2 ) is some polynomial of Q, R 1 , and R 2 explicitly calculable from V and α. This expression can be simplified in the thermodynamic limit and then be maximized e.g. by gradient ascent [34]. In practice, for transfer matrices with relativistic H like the ones we consider, this maximization has to be carried over matrices Q, R 1 , R 2 with a fixed maximum norm (or with a soft penalization of large norms). This is necessary to prevent the cMPS and its finite entanglement from capturing only the UV features of the stationary state [46]. For sufficiently large bond dimension χ, and taking into account this subtlety, we expect to get a good estimate of the stationary state. Once Q, R 1 , and R 2 are fixed this way, physical correlation functions can be computed analytically using equation (83). For example, if α is linear α(φ) ∝ φ we get, for x ≥ y: More complicated cases could be treated in a similar way. Through 2 successive dimensional reduction, we can thus compute certain correlation functions of a cTNS in d = 2 with an expression involving only matrices with a finite number of entries (d = 0). There is a priori no objection in principle to contract a d = 3 cTNS this way, but each additional dimensional reduction is done at the price of a variational optimization. For numerical purposes, the optimization of the cMPS is the crucial step. While current methods [46,[49][50][51] can be used, the prospect to use cMPS to solve field theories in more than 1 spatial dimension provides a strong additional motivation to make them more efficient.

Perturbation theory
Given that it is possible to compute correlation functions for Gaussian cTNS, it is natural to compute correlation functions for more general states by carrying a perturbative expansion around Gaussian states. One simply Dyson expands the non-Gaussian part of the exponential in the expression for the generating functional (13). It generically yields an expansion in terms of Feynman diagrams, similar to that of QFT.
For example, if α[φ] is linear in φ with a correction ∝ λ k φ k φ , the expansion will contain diagrams composed of vertices with 3 and 4 legs, corresponding to the α[φ]α * [φ ] term in (13), connected by Gaussian propagators. As previously mentioned, unless cancellations between different auxiliary fields occur, loop diagrams will be UV divergent and a regularization will be needed. We leave the derivation of the general Feynman rules, including a renormalization scheme, to future work. Notice that, in this approach, it is not the state itself that is defined through a perturbative expansion, but rather the correlation functions computed with it.

Others
There are of course many other ways one could compute correlation functions. As we mentioned before, one could rediscretize the cTNS to go back to a tensor network description, truncate the bond dimension, and use existing algorithms to contract it. However this would seem to partially defeat the purpose of introducing the continuum in the first place. An interesting avenue is to explore known approximations or tools of quantum field theory (besides perturbation theory) that would not be obvious in the discrete, like saddle point approximations, large D limits, or functional renormalization. Finally, direct Monte-Carlo sampling of the auxiliary field, although it will yield oscillating terms harming convergence in the general case, is a last resort option.

A. General metric and anisotropy
The main difficulty to overcome in order to construct cTNS lay in preserving local Euclidean symmetries. We may now wish to relax this constraint by allowing a general metric and anisotropic terms in the functional integral definition (1). Namely, it is natural to consider the following generalization.

Definition 4 (General functional integral formulation).
A continuous tensor network state (cTNS) of a bosonic quantum field on a smooth Riemanian manifold M with boundary ∂M and metric g, is a state |V, B, α parameterized by 2 functions V and α: R D+dD → C, and a boundary functional B: L 2 (∂M ) → C defined by the functional integral on an auxiliary D-component field φ: where all functions depend explicitly on position and summation on k is assumed.

B. Specialization: cMERA
A natural specialization of the previous generalization consists in having an auxiliary field living on an hyperbolic manifold M coupled to a physical field restricted to the boundary ∂M (see Fig. 5). Namely, we have in mind a state of the form: where the integral on the boundary may have to be taken as some appropriately rescaled limit of a bulk integral.
FIG. 5. Physical states on a boundary -In the discrete (left), the MERA is a (special case of) tensor network state with a hierarchical structure, with physical indices only at the boundary. Tensor network states with such a structure can also be extended to the continuum with our cTNS ansatz, by restricting the physical field to the boundary and choosing an appropriate metric (here hyperbolic) for the bulk auxiliary fields.
Such states could provide a natural generalization of the multi-scale entanglement renormalization ansatz (MERA) [32] in the continuum and for an arbitrary number of physical dimensions. They could provide a natural continuum versions of tensor network toy models of the AdS/CFT correspondence [29][30][31]. The form (93) is also reminiscent of field theory toy models of the AdS/CFT correspondence [52], where scalar field theories on a fixed AdS background are related to conformal field theories on the boundary. Note that this is approach is different in spirit from that of the standard entanglement renormalization approach to quantum fields [53][54][55][56][57], constructed as a unitary transformation applied on a QFT ground state. In the proposal (93), there is a straightforward lattice discretization and a natural "bulk" description in terms of auxiliary fields. However, the isometry property, characteristic of the MERA, is less straightforward to implement.

C. Fermions
We have defined our ansatz for bosonic quantum fields, because functional integrals and field coherent states are more natural in this context. To extend our proposal to fermions, one would have to introduce quite peculiar Grassmanian integrals with even kinetic and potential V terms but a Grassman odd term α in front of the creation operator ψ † . For fermions, it may be more convenient to start with an operator representation like that of (14), where Euclidean invariance is less natural, to subsequently derive the functional integral formulation.

D. Conformal field theory
We defined cTNS with the help of D auxiliary free massless scalar fields of measure dµ. A natural generalization would be to consider more general conformal field theories (CFT) for the auxiliary space, in the spirit of what has been proposed in the context of matrix product states with infinite bond dimensions [41]. Admittedly, some non-trivial measures can already effectively be emulated by tuning the real part of the potential V in (1). However in the general case, it may be more convenient to use the CFT machinery directly, for example on the wave function representation (7,17) or on correlation functions.

VII. DISCUSSION
We have put forward a new class of states for quantum fields that is obtained as a continuum limit of tensor network states and thus carries the same fundamental properties.
Although we have shown a number of interesting properties of our class of states, many interesting questions are so far open. Is it possible to find a quasi local parent Hamiltonian for such states? Can the transfer matrix T used to compute correlation functions in the operator representation be put in canonical form? Are there important gauge transformations our discussion in IV B ignores? How do V and α encode topological order and (local and global) gauge symmetries? Can this approach be combined with techniques developed on the lattice to study Gauge theories with tensor networks [58][59][60][61]? Are there non-trivial non-Gaussian cTNS for which correlation functions can be computed exactly? Do (possibly regularized) cTNS generically obey the area law like their discrete counterparts? Can cTNS be used to construct interesting toy models of the AdS/CFT correspondence? To what extent does the bond field dimension D quantify entanglement for (possibly only some) classes of cTNS?
Tackling these questions is an important goal for future work, to fully extend the success of tensor networks from the lattice to the continuum.