Finite Correlation Length Scaling in Lorentz-Invariant Gapless iPEPS Wave Functions

It is an open question how well tensor network states in the form of an infinite projected entangled pair states (iPEPS) tensor network can approximate gapless quantum states of matter. Here we address this issue for two different physical scenarios: i) a conformally invariant $(2+1)d$ quantum critical point in the incarnation of the transverse field Ising model on the square lattice and ii) spontaneously broken continuous symmetries with gapless Goldstone modes exemplified by the $S=1/2$ antiferromagnetic Heisenberg and XY models on the square lattice. We find that the energetically best wave functions display {\em finite} correlation lengths and we introduce a powerful finite correlation length scaling framework for the analysis of such finite-$D$ iPEPS states. The framework is important i) to understand the mild limitations of the finite-$D$ iPEPS manifold in representing Lorentz-invariant, gapless many body quantum states and ii) to put forward a practical scheme in which the finite correlation length $\xi(D)$ combined with field theory inspired formulae can be used to extrapolate the data to infinite correlation length, i.e. to the thermodynamic limit. The finite correlation length scaling framework opens the way for further exploration of quantum matter with an (expected) Lorentz-invariant, massless low-energy description, with many applications ranging from condensed matter to high-energy physics.


I. INTRODUCTION
The study of interacting quantum matter is of enormous current interest, with questions ranging from quantum spin liquids, topological matter, correlated electrons in solids, ultracold atoms in optical lattices to strongly coupled quantum field theories.
In this context numerical approaches play a very important role, with tensor networks being a central player. For problems in one spatial dimension methods such as the density matrix renormalization group (DMRG) [1,2], (infinite) matrix product states ((i)MPS) [3][4][5], and the multiscale entanglement renormalization ansatz (MERA) [6] have proven to be very powerful, both for gapped states and for quantum critical states with a low-energy conformal field theory description [7][8][9][10].
In two spatial dimensions infinite projected entangled pair states (iPEPS) [11] have become a competitive numerical approach with successful applications to many problems in the field of quantum magnetism and for strongly correlated fermions [12][13][14][15][16][17]. Furthermore theoretical work has established how different forms of topological order can be represented and understood in iPEPS wave functions with finite bond dimension [18,19].
It is however an open question how well tensor network states in the form of an iPEPS tensor network can approximate gapless quantum states of matter. Here we address this issue for two different physical scenarios: i) a conformally invariant (2 + 1)d quantum critical point in the incarnation of the transverse field Ising model on the square lattice and ii) spontaneously broken continuous symmetries with gapless Goldstone modes exemplified by the S = 1/2 antiferromagnetic Heisenberg and XY models on the square lattice. We find that the energetically best wave functions display finite correlation lengths and we introduce a powerful finite correlation length scaling framework for the analysis of such finite-D iPEPS states.
The outline of this paper is as follows: We start by providing a short introduction to the iPEPS framework and the energy optimization strategies in Sec. II. In Sec. III we study the critical properties of the (2 + 1)d transverse field Ising model as an example of a quantum critical point in the (2+1)d Ising universality class. In Sec. IV we present results for the S = 1/2 antiferromagnetic Heisenberg model and the S = 1/2 XY model, as examples for continuous symmetry breaking. In Sec. V we provide an extensive discussion and interpretation of the results obtained and we conclude in Sec. VI.

II. INFINITE (SYSTEM) PEPS
Considering a two-dimensional quantum many-body system consisting of p-level particles, placed on an infinite square lattice, one can make an ansatz for a wave function describing the system, |ψ = σ1,σ2,...
which is commonly known as infinite system projected entangled pair state (iPEPS) [11]. The graph in Eq. (1) is called a tensor network diagram, where nodes (edges) represent tensors (their corresponding indices). Edges connecting two tensors indicate summation indices and are of dimension D, which we call bond dimension of the iPEPS. The open indices σ 1 , σ 2 , . . . are of dimension p. In this work we only consider iPEPS with a single-site unit cell, i. e. iPEPS where the same tensor is used for each site, but all the techniques described here can be generalized to arbitrary unit cells.
iPEPS are a straight-forward generalization of infinite system matrix product states (iMPS) for two spatial dimensions and obey an area law for the entanglement entropy as well [11,20]. In contrast to iMPS, iPEPS can be constructed to have infinite correlation lengths, already for the simplest nontrivial case, D = 2 [21].

A. Contraction
The key challenge in all iPEPS algorithms is the so-called contraction of the state. For instance, to evaluate single-site observables, one needs the single-site density matrix, which consists of an infinite sum of double-layer tensors, The idea of a contraction is to find an approximation with controllable error for this infinitely large tensor network. There are several ways to go: Finding an approximate environment in the form of a boundary matrix product state (bMPS) [22] or a corner transfer matrix (CTM) environment [23] or by directly applying renormalization group (RG) schemes such as tensor RG (TRG) [24], tensor-entanglement RG (TERG) [25], second RG (SRG) [26], higher-order TRG (HOTRG) [27] or tensor network renormalization (TNR) [28]. Formally, a bMPS is nothing but an iMPS, with bond dimension χ, that is an approximation for the dominant eigenvector of a transfer matrix of double layer iPEPS tensors, with corresponding eigenvalue 1. Note that the quality of the approximation can be improved by increasing χ. From the singular values, ζ j , of the matrix C (Rényi) entropies, can be computed. In addition to the bMPS tensors one has to determine a horizontal dominant eigenvector, with corresponding eigenvalue 1 to be able to write the singlesite density matrix as The state-of-the-art method to find such a bMPS is described in Ref. [29] and has also been used in this work. Another way to contract a state, is to find corner transfer matrices (CTMs), and half-line transfer tensors (HLTTs), where again the contraction dimension χ is used to control the error of the approximation. With the CTMs and HLTTs, the single-site density matrix can be written as A powerful numerical tool to find CTMs and HLTTs is the socalled CTM renormalization group (CTMRG) [23]. The specific CTMRG procedure introduced in Ref. [30] has been used in this work, as it is a particularly stable variant of CTMRG. For this work both bMPS and CTMRG contractions have been implemented and we observe that for equal χ both methods give almost identical results.

B. Energy Optimization
As iPEPS are especially well-suited to describe ground states, energy optimization algorithms for iPEPS are of particular interest. For almost a decade there was only a single strategy for this task: imaginary time evolution [22] in various variants. Recently Refs. [31] and [32] introduced new direct variational approaches, which achieved lower energies than imaginary time evolutions, even in the limit of vanishing Trotter step size. Both variational methods rely on interaction contractions, similar to the norm contractions described in the previous subsection, but including all interaction terms of the Hamiltonian.
The first method [31] makes use of CTMRG interaction contractions to formulate a generalized eigenvalue problem (GEVP) for a given iPEPS, where the eigenvector corresponding to the lowest eigenvalue is used to propose the next iPEPS tensor in the minimization run.
However, it should be noted that also CTMRG interaction contractions can be used to obtain energy gradients and vice versa bMPS interaction contractions to obtain the optimization GEVP [37].
For iPEPS energy optimizations in this work bMPS interaction contractions in combination with the BFGS algorithm were primarily used, as this method turned out to be the most stable one. Also, the BFGS algorithm seems to be more stable close to the minimum compared to CG methods. Some states also have been optimized using a brute force minimization method, i. e. by computing finite difference energy gradients. All iPEPS optimized in this work have a single-site unit cell with a complex tensor which was forced to be invariant under spatial rotations. No symmetries at the virtual level were imposed. It turned out that starting with several random states and small contraction dimensions (16 χ 32) is the most economic way to bootstrap energy minimizations. The contraction dimension is then successively increased -in our optimizations up to values of χ = 512. For the transverse field Ising model we observed that it can also be beneficial to use intermediate minimization results to initialize minimizations for nearby values of the transverse field.

C. Correlation Length
A correlation function, c(r), of two arbitrary observables can be written as using a bMPS contraction, where the tensors u and v contain the specific observables. It should be noted that the same correlation function can be expressed by replacing the bMPS tensors with HLTTs from a CTMRG contraction. We write the eigendecomposition of the transfer matrix as with |λ 0 | ≥ |λ 1 | ≥ . . . and w. l. o. g. we assume λ 0 = 1. Then the dominant correlation length of the system, ξ, can be extracted as With this method the dominant correlation length can be extracted without even knowing which observables lead to the corresponding correlation function.
It should be noted that in contrast to local observables, the correlation length requires huge contraction dimensions χ to converge. Incorporating the ideas of Ref. [38] we observe the functional behavior which enables to extract a converged value for the correlation length, ξ(∞), precisely already for moderate values of χ. In cases where λ 1 λ 2 (due to degeneracy, as observed e.g. in the S = 1/2 Heisenberg model) one should use the largest eigenvalue different from λ 1 instead of λ 2 for the scaling in Eq. (14). In a first application we study the critical behavior of our variationally optimized iPEPS tensors for a quantum critical point in the (2 + 1)d Ising universality class. In the following we call this universality class described by a (2 + 1)d dimensional, Lorentz-invariant conformal field theory (CFT) the 3d Ising CFT. Note that this critical behavior is distinct from the one observed in the so called Ising PEPS [21], which amounts to promoting the thermal partition function of the 2d classical Ising model into a two-dimensional quantum many body wave function in PEPS form with bond dimension D = 2. This wave function can be parametrized by the temperature T entering the partition function, and at T = T c describes a critical wave function with algebraically decaying correlation functions. However the critical properties of this wave function are described by the 2d Ising CFT.

A. Overview
The Hamiltonian studied in this section is the transverse field Ising model on an infinite square lattice with an addi-  [39]. The lower six rows denote the additional critical exponents probed in this work, which are derived from the scaling dimensions above.
tional longitudinal magnetic field, where J = 1 denotes the strength of the ferromagnetic Ising interactions and sets the energy scale, h denotes the transverse field while h z parametrizes the longitudinal magnetic field. The phase diagram of this model is sketched in Fig. 1.
For h z = 0 the model has a Z 2 spin inversion symmetry and features a quantum phase transition at h/J = h c = 3.04438(2) [40], which separates a paramagnetic phase with m z ≡ σ z i = 0 for h > h c from a symmetry broken phase with m z = 0 for h < h c . The quantum critical point at h = h c , h z = 0 is described by the 3d Ising CFT. For all finite h z = 0 there is no critical behavior and the z-magnetization m z is finite as a response to the finite h z . In the entire phase diagram with h = 0 the transverse x-magnetization, m x ≡ σ x i , is finite. In the following we will explore the physics in the vicinity of the critical point h = h c , h z = 0 using iPEPS simulations. In order to provide field theoretical predictions to compare with we briefly review the properties of the 3d Ising CFT for our purposes. The 3d Ising CFT has two relevant perturbations, the O σ and the O field. Their scaling dimensions ∆ σ and ∆ are given in Tab. I. They clearly differ in 2d and 3d, which leads to distinct critical behavior. In the transverse field Ising model it is expected that in the continuum limit we can Since the perturbation is relevant, i.e. ∆ φ < d, it opens a mass gap proportional to λ 1 d−∆ φ , respectively it leads to a finite cor- In the transverse field Ising model it is understood that the microscopic cou- The 3d Ising CFT has two relevant perturbations, which translate into the two correlation length exponents, ν = 1/(d−∆ ) and ν σ = 1/(d−∆ σ ), which appear in the scaling We are in a situation where the perturbed theory can display magnetizations, and we thus define critical exponents α φ , β φ which describe how the magnetizations grow as a function of the perturbing coupling for a perturbation φ. α φ is the exponent we use when measuring m The subscript σ, denotes the perturbing field as for the correlation length exponents before. For the transverse magnetization we expect |m x,c − m x | ∼ |λ| α φ , while for the longitudinal magnetization we expect m z ∼ |λ| β φ . The definitions and values for these exponents can be found in Tab. I. Some of the exponents we have introduced here for clarity reasons are also more commonly known as β ≡ β , δ ≡ 1/β σ and ν ≡ ν in the statistical mechanics literature. Note that the d = 4 critical exponents are equivalent to the mean field exponents. This is because d = 4 is the upper critical dimension for the Ising criticality, i.e. the dimension where mean field theory becomes exact.

B. iPEPS Results for the Transverse Field Ising Model
The transverse field Ising model (with h z = 0) has been studied extensively with both finite-size PEPS and iPEPS approaches in the past [22,23,32,[41][42][43][44][45][46]. A common feature of all the these simulations is that for h below a bond dimension D dependent value, h (D) c , the system shows a finite zmagnetization m z . However in the past the values of h (D) c and the functional behavior of m z in its vicinity did depend significantly on the tensor optimization methodology. We believe the newest generation of optimization algorithms put forward in Refs. [31,32] do not suffer from these shortcomings anymore, so that a detailed analysis of the intrinsic iPEPS finite-D behavior close to the 3d Ising CFT is finally possible.

Magnetizations in the vicinity of the critical point
We start by presenting in Fig. 2 [40]. In an earlier study based on one-dimensional iMPS states for the (1 + 1)d transverse field Ising model, a bond dimension D > 10 was required to reach a comparable accuracy [43]. We have also tried to optimize D = 4 tensors, but albeit technically possible, it turns out to be extremely difficult to obtain energies which are lower than our best D = 3 results so far. We will come back to this observation later in this section. Finally in panels (c) and (d) we display m z and m x along an orthogonal cut at fixed h = h c with varying h z > 0, i.e. along the violet axis in Fig. 1.
While the plots in Fig. 2 seem to suggest large differences between D = 2 and D = 3 it should be noted that the h and h z ranges shown are already quite small. Shown on a scale h ∈ [0, 4] it would be difficult to visually spot the differences between the two D values.

Critical exponents
In a next step we explore the critical behavior contained in the presented data. In Fig. 3 line with the expected slope β as a guide to the eye. In the corresponding inset we numerically calculate the derivative and find a collapse between the D = 2 and D = 3 data at larger distances from the critical point. The D = 2 running estimate for the critical exponent reaches a maximum of ∼ 0.29 and then drops to zero as h c − h goes to zero. The D = 3 running estimate rises to ∼ 0.32 before it also drops to zero as h c − h goes to zero. We expect that D > 3 would get even closer to the expected β ≈ 0.326 418 (2), before dropping to zero as h c − h goes to zero. The drop to zero is clearly due to the finite-D remnant m z at the thermodynamic value h c . Let us therefore investigate what happens when we analyze the data as a function of h Fig. 3 displays the corresponding data. The D = 2 data clearly shows a limiting mean-field behavior β − h (in the inset), as observed previously in Ref. [43]. The D = 3 data shows some hint of an intermediate plateau around the true 3d Ising CFT value for β before also crossing over to the mean-field value β So we learn that when the location of the critical point is known beforehand, the critical behavior can be determined quite accurately already with a surprisingly small bond dimension of D = 3. In the vicinity of the finite D critical points we however observe mean field behavior and due to the crossover it is more difficult to extract the genuine critical behavior.

Correlation lengths
After having analyzed the critical behavior of local observables as a function of perturbing couplings, we now investigate the correlation lengths in our optimized iPEPS wave functions in the vicinity of the critical point. In the vicinity and at a quantum critical point in (1 + 1)d represented with a finite bond dimension iMPS we know that only a finite correlation length can appear. Since iPEPS is in principle able to represent wave functions with algebraically decaying correlations, i.e. states with infinite correlation lengths [21], it is not a priori clear what to expect in our optimized iPEPS wave functions. Let us note first that the correlation lengths for D = 1 (product states) vanish identically, even though the spontaneous magnetization shows critical mean-field behavior at h , h z = 0. We are quite confident that this short correlation length is not an artefact of incomplete optimization, but is a genuine feature of optimized, translationally invariant, finite-D iPEPS wave functions for Lorentz-invariant quantum critical systems with a 3d space-time description. For D = 3 we also observe finite correlation lengths, but now the maximum is substantially larger: ξ . So both D values seem to indicate that our variational optima feature a finite correlation length. This is one of the key results of this paper, whose possible origin we are going to discuss later. We will however show in the following, that the finite correlation length is also a blessing, as it helps us to understand and organize the finite D effects in field theoretically motivated formulas based on ξ(D).
Before doing this, let us investigate the functional behavior of the correlation length in the vicinity of the critical point. Depending on the cut in parameter space we expect ξ ∼ |h c − h| −ν or ξ ∼ |h z | −νσ , with the values of ν and ν σ given in Tab. I. Indeed the data in panels (b) and (d) of Fig. 4 shown on a log-log scale are (roughly) compatible with the expected correlation length exponents in some intermediate window of the couplings. This is expected since far away from the critical point we are outside the quantum critical regime, while very close to the critical point ξ saturates. Still the agreement for the h z -detuning is much better than the |h c − h|-detuning.
Finally we plot the expectation values of the two magnetizations m z and |m x,c − m x | as a function of the measured correlation length ξ (for both parameter cuts) in panels (e) and (f) of Fig. 4. As discussed earlier we expect this relation to be governed by the scaling dimensions ∆ σ and ∆ respectively. While the D = 2 results do not match too well, the D = 3 results for both parameter cuts and both observables are in good agreement with the expected scaling dimensions.
Even though it seems that the critical exponents and scaling dimensions can be obtained more precisely based on the observables as a function of the couplings than of the correlation lengths, it is nevertheless rewarding to observe that the correlation lengths are also following the expected quantum critical behavior with increasing D, within the stated limitations.

Critical energy convergence
We now investigate the energy convergence of the TFI model at its critical point h = h c , h z = 0 for increasing bond dimension D. It is one of the open problems in practical iPEPS calculations to understand the convergence of energies as a function of D. Here we advocate that the variational energy of an optimized iPEPS tensor at the critical point h = h x , h z = 0 can be understood as a particular type of a Casimir effect controlled by the correlation length ξ. It is well known that the ground state energy density e = E/N sites of a 3d quantum critical system in a torus geometry with modular parameter τ [48] is given as where L denotes the linear extent of the torus, v is the "speed of light", e.g. the critical spin wave velocity in a TFI model, and α QCP τ is a τ dependent Casimir amplitude which otherwise depends only on the universality class of the quantum critical point (QCP), for example for the 3d Ising CFT and a square torus with periodic boundary conditions α (3d Ising CFT) τ =i = +0.35 (2) according to Ref. [49]. We now postulate that our iPEPS setup can be considered as a new, distinct geometry with its own Casimir amplitude α (3d Ising CFT) iPEPS , where however the length of the torus is replaced by the correlation length ξ, such that We stress that this Ansatz is in agreement with the expected scaling behavior of the one-point function of the stress-energy tensor T , whose scaling dimension in d = 3 is ∆ T = 3. is itself three orders of magnitude smaller than the square torus amplitude. So multiplying these two factors we are led to conjecture that the extrapolated iPEPS energies are accurate to about 10 −6 already at D = 3! This result might explain why D = 4 simulations are so demanding, as the expected remaining energy gains are tiny and are accompagnied with wave functions bestowed with large correlation lengths, which are even harder to contract accurately.
In order to corroborate the advocated scaling Ansatz in Eq. (17)

IV. CONTINUOUS SYMMETRY BREAKING
In this section we want to explore the properties of iPEPS wave functions as they represent or approximate quantum many body states which exhibit continuous symmetry breaking in (2+1) space-time dimensions. This is a rather ubiquitous phenomenon, ranging from magnetic order in O(N ) symmetric quantum magnets with N ≥ 2, to bosonic and fermionic superfluids, to superconductors and Goldstone phases in high-energy physics. In these systems the continuous symmetry is spontaneously broken, accompanied by the the appearance of a finite order parameter |M| > 0. Another hallmark is the required presence of gapless excitations (so called Goldstone modes), which are the soft long-wavelength modes of order parameter variations. We will study the S = 1/2 Heisenberg antiferromagnet and the S = 1/2 XY model, both on the square lattice, as paradigmatic examples for O (3) and O(2) continuous symmetry breaking with a three-respectively two-component vector order parameter.

A. Overview
The field theoretical description of collinear magnetic order in O(N ) quantum magnets relies often on a quantum nonlinear sigma model (NLSM) formulation, or on a description in terms of the symmetry breaking phase in an N -component interacting φ 4 theory. The gapless Goldstone modes are also known as spin waves in the magnetic context and taken in isolation they behave as a collection of free massless scalar fields. These gapless modes are responsible e.g. for the algebraic decay of spin-spin correlations to their limiting value and the finite size corrections of the energy or the order parameter. These linearly dispersing modes require a (2+1)d Lorentz-invariant description at low energies, similar to the Ising CFT discussed before. The quantum non-linear O(N ) sigma model is described in detail in Ref. [51]. For our purpose it is sufficient to know that this is basically a hydrodynamic theory of quantum magnets with collinear order. Its description, being hydrodynamic, relies only on a handful of effective parameters entering the description: the spin stiffness, ρ s , the spin wave velocity v, the transverse susceptibility, χ ⊥ , and the squared order parameter in the thermodynamic limit, m 2 (∞). The first three parameters are actually related via v 2 = ρ s /χ ⊥ . Similar to the finite size corrections to the ground state energy discussed in the quantum critical context, the finite size corrections to the ground state energy and the order parameter have been derived for the quantum non-linear O(N ) sigma model in Refs. [52][53][54][55][56]. The finite size corrections to the ground state energy e = E/N sites in d = 3 are as follows: For a square torus with periodic boundary conditions α NLSM τ =i ≈ 1.437745 has been obtained. The finite size correction for the magnetic order parameter squared are as follows: (20) For a square torus with periodic boundary conditions µ NLSM τ =i ≈ 0.62075 has been found. As we will see below, our variationally optimized iPEPS wave functions have a finite correlation length ξ, which depends on the model and on the bond dimension D. We now conjecture the following d = 3, finite ξ corrections for the ground state energy density e and the magnetic order parame-ter squared m 2 , and The potential power of these formulas lies in the fact that one can extrapolate the results at finite ξ(D) to the limit ξ → ∞ based on iPEPS data fits to the above formulas. Furthermore it is possible to predict the finite size corrections for other microscopic models, once the "universal" values of α NLSM iPEPS and µ NLSM iPEPS are determined. The knowledge of N , v and v/ρ s allows then to quantitatively predict the slope of the finite size corrections. Conversely it might become possible to estimate v and v/ρ s based on precise iPEPS data for a model with a known value for N .

B. S = 1/2 antiferromagnetic Heisenberg model
The S = 1/2 antiferromagnetic Heisenberg model has also been studied frequently using finite size PEPS and iPEPS approaches in the past [31,32,41,44,[57][58][59]. The ground state of this model has antiferromagnetic Néel order, which breaks the continuous O(N =3) rotation symmetry spontaneously down to a residual O(2) symmetry. The presence of N − 1 = 2 Goldstone modes leads to an algebraic decay of the two-spin correlation function.
The Hamiltonian is defined as with J = 1 and the S α i are spin 1/2 operators. In order to be able to work with a single-site unit cell we perform a spin rotation on one Néel sublattice, which negates the sign of the y and z parts of the interactions in the actual calculations.
We proceed optimizing the variational energies of iPEPS tensors for D = 2, 3, 4, 5. Then we measure the order parameter squared m 2 = | σ i | 2 (i.e. the maximum possible for the order parameter squared amounts to one), as well as the correlation length ξ. These correlation lengths are finite and range from ξ(D = 2) ≈ 0.7 to ξ(D = 5) ≈ 4.5. The correlation lengths are thus substantially smaller than those of the critical transverse field Ising model at D = 3.
In Fig. 6 we present the energy per site e as function of 1/ξ 3 in the left panel and m 2 as a function of 1/ξ in the right panel. This is the conjectured ξ scaling form of Eqs. (21) and (22). It is striking that for both observables a linear fit leads to accurate extrapolations to the limit ξ → ∞, when compared to high precision QMC reference values [56,60]. We fit the largest three D values for the energy e and all the D values for m 2 and obtain the following iPEPS ξ → ∞ estimates: e HB ≈ −0.669417 and m 2 HB ≈ 0.380. Using the iPEPS fit slopes, the value N = 3 and the known QMC values of the hydrodynamic parameters v and v/ρ s [56,61], we can then proceed to determine and µ NLSM iPEPS ≈ +0.045.
Note that similar to the Ising Casimir amplitude, the NLSM energy Casimir amplitude is three orders of magnitude smaller than the square torus one, highlighting that an iPEPS calculation at a certain ξ should be considered three orders of magnitude more accurate than a square torus with L ∼ ξ. The iPEPS order parameter amplitude is however only one order of magnitude smaller than the square torus result.
The S = 1/2 XY model has also been investigated with PEPS and iPEPS in the past [32,57,62]. The ground state of this model has ferromagnetic order in the x−y spin plane, which breaks the continuous O(N =2) in-plane rotation symmetry spontaneously down to a residual discrete Z 2 symmetry. The presence of N − 1 = 1 Goldstone mode leads to an algebraic decay of the two-spin correlation function.
The Hamiltonian is defined as, with J = 1. Note that for this model the two choices of the sign of J can be mapped into each other. Since we want to use a single site iPEPS unit cell we adopt the ferromagnetic sign convention. We proceed optimizing the variational energies of iPEPS tensors for D = 2, 3, 4, 5. Then we measure the order parameter squared m 2 = | σ i | 2 (i.e. the maximum possible for the order parameter squared amounts to one), as well as the correlation length ξ. These correlation lengths are again finite and range from ξ(D = 2) ≈ 1.3 to ξ(D = 5) ≈ 8.4. The correlation lengths for each D are roughly a factor two larger than for the Heisenberg model. So it seems that the smaller number of Goldstone modes has a beneficial effect on the growth of the correlation lengths.
In Fig. 7 we present the energy per site e as function of 1/ξ 3 in the left panel and m 2 as a function of 1/ξ in the right panel. This is the conjectured ξ scaling form of Eqs. (21) and (22). We then fit the largest two D values for the energy e and all the D values for m 2 and obtain the following iPEPS ξ → ∞ estimates: e XY ≈ −0.548822 and m 2 XY ≈ 0.757. These values compare well with the QMC results of Ref. [63]. Since the QMC results are relatively antiquated, it is not inconceiv- able that the iPEPS results are actually more precise than the QMC results.
We are now also in a position to corroborate the O(N ) universality of the conjectured NLSM finite ξ corrections. Using the hydrodynamic parameters v and v/ρ s from Ref. [63], the iPEPS amplitudes from Eqs. (24) and (25) and inserting N = 2 we arrive at the NLSM predictions, which are shown by solid lines in both panels of Fig. 7. The nice agreement between the linear fits and the NLSM predictions for the slopes provide further support for the validity and therefore power of the field theoretically inspired finite-ξ correction formulae.

V. DISCUSSION AND INTERPRETATION
After having studied the three different models we are confronted with the fact that in all cases the correlation length ξ(D) was finite. While we have developed a powerful finite ξ scaling framework, where many observables, including energies and order parameters, can be analyzed and extrapolated to ξ → ∞, we are still left in the dark both regarding the underlying origin of the finite ξ(D) in the first place, and regarding the functional dependence of ξ on D for a given model or universality class.
While we don't yet have compelling answers to both questions so far, we can at least try to shed as much light as possible based on our numerical data. In Fig. 8 we have assembled the correlation lengths observed in the critical transverse field Ising model, which is described by one of the simplest nontrivial 3d CFTs, together with two instances of O(N ) continuous symmetry breaking phenomena with N = 3 for the S = 1/2 antiferromagnetic Heisenberg model and N = 2 for the S = 1/2 XY model. In the inset we also show the D ≥ 2 data in a log-log plot, yielding some rough estimates for a putative power-law relation ξ(D) ∼ D κ (it is not clear that such a law holds). With the two points for the critical transverse field Ising model we obtain κ (3d Ising CFT) ≈ 5, while both continuous symmetry breaking cases seem to share the same κ NLSM ≈ 2. The latter two cases however differ by a factor two in the prefactor, which incidentally is also the inverse ratio of the number of Goldstone modes. This could mean that the XY model has twice as large a correlation length as the Heisenberg model because it only has half the number of Goldstone modes. It would be interesting to explore whether the known additive logarithmic correction to the entanglement area law in continuous symmetry breaking states [64] might be at the origin of this behavior. There the prefactor of the log contribution is proportional to the number of Goldstone modes. This would also explain why the speculative values of κ are so different between the 3d Ising CFT and the symmetry breaking cases. In terms of their low-energy degrees of freedom, the 3d Ising CFT and the XY model both contain a single real scalar field each -an interacting field in the Ising CFT case and a massless free field in the XY case. One could thus have expected the values of κ to be roughly similar. The crucial difference therefore seems to come from the broken continuous symmetry in the XY model, which is absent in the Ising CFT.
The most pressing question remains however as to why we only find finite correlation lengths in the variationally optimized iPEPS wave functions for massless Lorentz-invariant (2 + 1)d scenarios. We believe that our observations are actually the generic result, and that the previously known examples of iPEPS states with algebraic correlations are fine-tuned and non-generic. As shown as a cartoon in the left panel of Fig. 9, we think of our (2+1)d iPEPS states as wave functions whose correlation functions are represented by a path integral with a finite ξ τ (D) extent in the (real or imaginary) time direction. The Lorentz-invariance (or Euclidian invariance after a Wick-rotation) of the fixed point we try to approximate then forces the spatial correlation lengths to be finite as well. The well known iPEPS states with algebraic spatial correlations at finite D can actually be represented by a purely in-plane path integral, where the temporal extent ξ τ is basically zero (right panel). This is certainly true for the 2d classical partition function Ising PEPS [21], quantum dimer Rokhsar-Kivelson states [18,65] and certain quantum Lifshitz theories [66,67], where there is a built-in space-time asymmetry. Some further evidence supporting this picture might be obtained from the entanglement entropies S bMPS of the boundary MPS resulting from the contractions of the iPEPS, as it seems plausible that this entanglement entropy is amplified if the correlation volume extends into the τ direction. The data shown in Fig. 10 can indeed be interpreted that the entanglement entropy grows more rapidly with the correlation length ξ in the genuine 3d space-time cases, compared to the quantitatively well understood logarithmic scaling of the (2 + 0)d wavefunctions, here exemplified by the Ising PEPS at various temperatures [68]. It would be interesting to study these boundary MPS entropies more systematically, as they roughly dictate how large the boundary χ value in the contraction schemes has to be.
Coming back to the reason for the finite ξ τ (D) in the first place: We speculate that finding the dominant eigenvector of a plane-to-plane transfer operator along the temporal direction combined with a projection to a finite D iPEPS leads invariably to a finite correlation length in the temporal direction. This view is also supported by a more formal argument stating that an entropic area law does not automatically imply an efficient iPEPS representation [69].
To a first approximation we understand the finite-D iPEPS to mimic a ground state wave function of the field theoretical fixed point perturbed by a certain amount of the most relevant allowed perturbation given the symmetry constraints imposed on the iPEPS. Such states are driven away from the gapless situation, leading to a finite correlation length ξ(D). In the examples studied here this always corresponds to a perturbation by coupling to the magnetic order parameter. Imposing symmetry constraints on the tensor might change the nature of the allowed relevant perturbation, and could also affect the values of some of the finite ξ correction amplitudes introduced in this work.

VI. CONCLUSIONS
In this paper we have introduced a powerful finite correlation length scaling framework for the analysis of finite-D iPEPS which have been variationally optimized for Lorentzinvariant (2 + 1)d quantum critical or continuous symmetry breaking Hamiltonians. This framework is important i) to understand the mild limitations of the finite-D iPEPS manifold in representing Lorentz-invariant, gapless many body quantum states and ii) to put forward a practical scheme in which the finite correlation length ξ(D) and field theory inspired formulae can be employed to extrapolate the data to infinite correlation length, i.e. to the thermodynamic limit. We have shown that some of the amplitudes entering the equations have a field theoretical interpretation and are therefore universal to some degree. We have demonstrated the power of the method for the energy convergence in all three considered models, and for order parameter extrapolations for the continuous symmetry breaking models, where the previously employed 1/D extrapolation schemes performed poorly. We believe that an-other advantage of this approach is also that the variational quality of an iPEPS tensor and the resulting correlation length (and other observables) seem to go hand in hand in the vicinity of a local optimum, in such a way that different data points still lie on the same finite ξ extrapolation curve.
We have also carefully analyzed the critical behavior of the transverse field Ising model in the vicinity of its critical point by measurements of local observables as a function of the two magnetic fields and we were able to obtain critical exponents within a few percent of the most recent conformal bootstrap results already at bond dimension D = 3.
The finite correlation length scaling framework opens the way for further exploration of quantum matter with an (expected) Lorentz-invariant, massless low-energy description. Within the quantum critical or CFT related questions one could explore the Wilson-Fisher O(N ) fixed points beyond the N = 1 case studied here, Gross-Neveu-Yukawa universality classes arising in interacting Dirac fermion systems, or QED 3 related behavior of gapless quantum spin liquids or deconfined criticality. In the context of continuous symmetry breaking various superfluid and superconducting phases in bosonic and fermionic systems should be described by the O(2) symmetry breaking results established in this work. Further directions are non-collinear magnetic order or SU(N ) continuous symmetry breaking. Ultimately one should also explore the occurence of finite correlation lengths and their scaling in interacting systems with a Fermi surface.
On the fundamental level it will be important to develop an understanding of the ξ(D) relation, that we started to explore here. Our preliminary results indicate a ξ ∼ D κ scaling with values of κ which could be universal. In the (1 + 1)d iMPS context, the values of κ depend only on the central charge c [8]. It would be interesting to understand in (2+1)d whether κ also depends only on some universal data, such as in the Ftheorem [70,71].