Computing the renormalization group flow of two-dimensional $\phi^4$ theory with tensor networks

We study the renormalization group flow of $\phi^4$ theory in two dimensions. Regularizing space into a fine-grained lattice and discretizing the scalar field in a controlled way, we rewrite the partition function of the theory as a tensor network. Combining local truncations and a standard coarse-graining scheme, we obtain the renormalization group flow of the theory as a map in a space of tensors. Aside from qualitative insights, we verify the scaling dimensions at criticality and extrapolate the critical coupling constant $f_{\rm c} = \lambda / \mu ^2$ to the continuum to find $f^{\rm cont.}_{\rm c} = 11.0861(90)$, which favorably compares with alternative methods.


Sec. I | Introduction
Solving interacting quantum field theories (QFTs) out of the perturbative regime is difficult, and is arguably one of the most pressing computational problems in theoretical physics nowadays.Among all non-trivial quantum field theories, the simplest one is perhaps the self-interacting scalar field in two spacetime dimensions, a.k.a.φ 4  2 .This theory provides a good case study because it is reasonably simple to define (it is super-renormalizable), and was actually one of the first interacting QFTs to be rigorously constructed [1,2].However, it does not bypass the real difficulties of QFT, namely it is neither free, nor integrable, nor supersymmetric.Furthermore, its physical interpretation is rather simple as it does not display multiple instantons or topological complications.Nevertheless, it does contain a second order phase transition, which is ubiquitous in physics.Finally, and perhaps most importantly, there are still some questions for which we do not possess exact answers.
Recently, there has been a renewal of interest for deterministic methods that aim at probing the deep nonperturbative regime of QFTs, and of φ 4  2 in particular.One can mention Hamiltonian truncation techniques as well as their renormalized refinements [3][4][5], Borel resummation [6], or tensor network based approaches [7][8][9].Ultimately, one common objective of all these methods is to outperform lattice Monte-Carlo techniques, which are primarily used in non-perturbative high-energy computations, e.g. in quantum chromodynamics [10][11][12].Although some of these new techniques temporarily outperformed the previous state of the art in the case of φ 4  2 , Monte Carlo methods [13,14] were quickly improved and now seem back on a par [15].As part of this effort, we propose in this manuscript a computation of the renormalization group (RG) flow of φ 4  2 using tensor network renormalization techniques.
Tensor network methods, which have proven a very potent tool in the context of quantum information and condensed matter theory, can be divided into two distinct families.In the first, one writes an ansatz for a lowenergy state of a given Hamiltonian as a tensor network, and optimizes its parameters variationally so as to get closer to the true ground state [16][17][18][19][20][21][22][23][24].In the second, one decomposes a Euclidean partition function directly as a tensor network, and the challenge then becomes to contract it approximately [25][26][27][28][29][30][31], with the aid of techniques akin to Kadanoff's blocking scheme [32].In both cases, and as the name suggests, tensor network methods have been primarily developed for lattice systems.To deal with QFTs, two options are available: Either one extends the tensor networks themselves to the continuum, as has been tried in recent years [33][34][35][36], or one discretizes the theory on a lattice, making the well-established tensor network techniques available.Although the first option is a promising initiative, we shall explore the second one in this manuscript.In short, we wish to study φ 4  2 by discretizing its functional integral on a very fine-grained lattice, and solving it with tensor network renormalization techniques.
The best-known method for approximately contracting tensor networks that represent partition functions is the so-called Tensor Renormalization Group (TRG) algorithm.It was applied to φ 4  2 in [37], and more recently in [9], where the authors report an estimate of the critical coupling of the theory in the continuum limit.Although it is a simple and fairly efficient algorithm, it has been known for several years that, despite its name, TRG does not in fact yield a proper renormalization group flow.Indeed, it was shown to violate a fundamental principle of RG, namely that all ultraviolet (UV) details relative to a scale below the one considered should be integrated away.Concretely, this means that the fixed points of the flow depend on non-universal features.This makes the identification of critical points more difficult, and the scale invariance taking place at conformal fixed points nearly impossible to witness.
Over the years, several more advanced algorithms that address the shortcomings of TRG, have been proposed [28][29][30][31].In this paper, we use the algorithm that was last introduced, namely Gilt-TNR, which is based on the notion of graph independent local truncations.In addition to being physically sounder than TRG, it yields notably more accurate observables.Another technical contribution of our study is the introduction of a more transparent and controlled field discretization.
Our main result is the computation of the critical value of the coupling in the continuum limit.As a benchmark, we also compute the scaling dimensions of the theory, which are known analytically from conformal field theory.Apart from the arguable gain of precision brought by tensor network renormalization techniques, we shall emphasize how these methods provide some precious qualitative insight about the RG flow of the theory.Furthermore, this study permits us to highlight some interesting specificities of the tensor network approach.In particular, we will show that they are rather complementary to other approaches: They work well at strong coupling, with excellent behavior in the infrared (IR), but rather suffer in the UV if the fixed point is the massless free boson, which is the case for super-renormalizable and asymptotically free bosonic QFTs.
In sec.II we recall the continuous and discrete formulations of the theory.The tensor network reformulation is explained in sec.III.In sec.IV we present the Gilt-TNR algorithm as a tool to perform renormalization group transformation for tensor networks, and the method pursued to identify the critical coupling.The algorithm is applied to φ 4  2 and the results are tabulated in sec.V. Finally, in sec.VI, we discuss the possible source of errors, as well as generalizations to higher dimensions.

Sec. II | The model
In this section, we define φ 4  2 theory on the lattice and in the continuum, and recall some of its basic properties.

II.A. From the continuum to the lattice
At a rather heuristic level, one can see (Euclidean) φ 4  2 theory as the theory of a random scalar field φ of probability measure with and Z = Dφ exp(−S E (φ)).Of course, taken at face value, this probability measure is ill-defined.One option forward is to regularize the theory by discretizing the continuum R 2 into a lattice, localize the field on its vertices, and write where a is the lattice unit length, i, j denotes a sum on nearest neighbors and Naively, if there were no quartic term, taking the continuum limit would simply amount to sending a → 0, while keeping the dimensionful parameter m 0 fixed.However, it is well understood that taking the continuum limit in the presence of non-quadratic terms requires a more subtle scaling of the parameters of the lattice theory.

II.B. Perturbative renormalization
To take the continuum limit, one has to take into account the fact that the leading effect of the quartic term is to shift the quadratic (mass) part by a term ∝ − log(a).Unless this effect is countered, the continuum limit is a trivial theory with an infinite mass.
One option is to add a mass counterterm and define where A(µ 2 ) corresponds to the (lattice) tadpole diagram responsible for the logarithmic divergence: .
It can be shown that this prescription is equivalent to the normal ordering used in the operator representation [38].The continuum limit is then obtained by sending a (or equivalently λ) to 0, while keeping f = λ/µ 2 fixed.In the case of φ 4 2 , this simple mass renormalization is the only requirement to make the continuum limit well defined.

II.C. Universal and non-universal properties
Let us now summarize what is known about the behavior of φ 4  2 , and distinguish its universal and non-universal properties.
First, the lattice theory has two phases separated by a second order phase transition along a line in (λ, µ) space.For a fixed λ, this phase transition is reached for a value µ c (λ) = 0. Whereas µ is a renormalized mass, since it contains a counter-term necessary to make the theory well defined in the continuum limit, it is not the physical mass associated with the spectral gap, which is a fully non-perturbative quantity.The continuum field theory also has a phase transition for f cont.c := lim λ→0 λ/µ 2 c (λ) ≡ g/m 2 c .Since φ 4 2 is a super-renormalizable QFT, defining its continuum limit is straightforward.It is a UV problem that is weakly coupled, and thus perturbatively solvableonly one loop is required.However, finding the position of the critical point, in the discrete setting as well as in the continuum theory, is a difficult problem.Indeed, it is an IR problem, and thus strongly coupled and nonperturbative [2].The most recent Monte-Carlo study [15] provides the estimate f cont.c 11.055 (20), which is more accurate than Borel-resumming the perturbation expansion to 8 th order in λ [6].
That being said, the IR properties of φ 4 2 exactly at criticality are known analytically from conformal field theory (CFT).The critical point is in the Ising universality class, and is the same on the whole (λ, µ c (λ))-line.Computing precisely the scaling dimensions on the critical line thus yields no new knowledge, but provides a good sanity check for any method aiming to solve φ 4  2 .

Sec. III | Tensor network reformulation
In this section, we explain how the partition function of the discretized theory can be computed as the contraction of a discrete tensor network.

III.A. Tensor network with field indices
Given the discretized action S a E , we consider the lattice partition function Assuming that we accept continuous indices, this partition function can be rewritten as a tensor network [9], which in the case of a 2 × 2 lattice reads where we introduced the tensors In these diagrams, a connected link means a "continuous contraction", i.e. an integral R dφ over real numbers.Equivalently, one could eliminate the Dirac deltas and define that each vertex corresponds to a field integral.

III.B. Field discretization
To make the previous tensor network with continuous indices amenable to standard techniques, it is necessary to discretize the field variable.One could simply digitize the field φ [37], or use Gauss-Hermite quadratures as was done in [9].We choose another way, which we feel is more transparent and gives a better control of the error.
It is not possible to carry out the field integrals independently of each other because of the exp(φ φ r ) term coming from the expansion of the kinetic term in F .We can simply solve this problem by expanding the exponential in Taylor series which allows to carry out all the integrals independently.This yields a factorization of F that can be pictorially represented as where This effectively replaces integrals by discrete (infinite) sums.Let us perform this factorization for every F in the tensor network: where we only represent a patch of the network, leaving the periodic boundary conditions implicit.Let us now define the elementary tensor T abcd with discrete indices: so that the previous tensor network expression reads For a 2 × 2 square lattice, the partition function can now be depicted as The indices (a, b, c, d) ∈ N 4 summed over still need to be restricted to a finite range K min , K max for numerics.The advantage of this expansion is that the truncation error is precisely controllable and decreases faster than exponentially in the cut-off K max (see app.A).Finally, the coefficients T abcd can be easily evaluated in terms of special functions, namely for n ∈ N where the confluent hypergeometric function of the second kind U (α, β, z) is defined as The latter has efficient implementations in most programming languages, and in particular in mpmath [39].
For half integers, h vanishes, and thus T abcd = 0 only if a + b + c + d is even, which is a consequence of the Z 2 symmetry of the problem.We exploit this feature in our algorithm by defining Z 2 symmetry preserving tensors as described in [40,41], which significantly reduces the computational time.

Sec. IV | Renormalization in the tensor network language
In this section, we explain how the renormalization group flow of φ 4 2 can be computed in the context of tensor networks using the Gilt-TNR algorithm.For a more general and detailed description of this algorithm, we encourage the reader to consult [31].

IV.A. Heuristics
We explained in the previous section how the partition function of the theory can be decomposed as a homogeneous tensor network, whose underlying graph is a square lattice with periodic boundary conditions.In this context, computing the partition function boils down to contracting the corresponding tensor network.In order to obtain an accurate prediction of the partition function in the thermodynamic limit, the network should be taken to be as large as possible.
Since we are dealing with a homogeneous network, it could be systematically contracted by first blocking a square of four tensors, and then iterate the process after replacing the original tensor with the result of the first blocking.Given an initial bond dimension χ along all the legs, one such round of contractions would result in a new tensor with bond dimension χ 2 , i.e.
The cost of such an operation would therefore grow exponentially in the lattice size, meaning only relatively small networks could be contracted in practice.
In order to deal with larger systems, some approximations need to be made, in such a way that the computational time of the contraction scheme only grows polynomially in the lattice size.If we are merely interested in the actual value of the partition function, then many options are available.However, our focus is rather on the renormalization group of the theory.This means that we want a method allowing us to define a flow in the space of tensors, obtained by defining effective tensors representing the same system at different length scales.
In order to obtain an effective tensor that represent the infrared length scale of the theory, the bond dimension must be prevented from growing at every iteration of the blocking scheme.In general, this is done by replacing a local patch of tensors by another one, for which the bond dimensions along some of the legs have been reduced.Naturally, such a replacement must only induce a small error, otherwise the benefit of dealing with a larger lattice size would not outweigh the penalty of performing approximations.This is guaranteed as long as some elements of the tensor account for short-range physics only, and can therefore be safely discarded as they become irrelevant at larger scales.

IV.B. Coarse-graining
The first algorithm to follow the philosophy described above was the so-called Tensor Renormalization Group (TRG) algorithm [25].This algorithm relies on replacing single tensors by their truncated singular value decomposition.Given a p × q matrix M , its singular value decomposition (SVD) is a factorization of the form M = U SV † , where U † U = 1 p×p , V † V = 1 q×q and S is p × q diagonal matrix.The entries S a := S ab δ ab ≥ 0 are known as the singular values of M .A truncated singular value decomposition is then obtained by discarding some of the smallest singular values in the factorization above, i.e. it is a low-rank approximation of the form where χ < rank(M ).This definition is easily generalized to tensors and can be schematically illustrated as , where the dotted lines indicate that there is an arbitrary number of legs on both side of the tensor A. By blocking {U, S 1 2 } on the left-hand side, and {S 1 2 , V } on the righthand side, we obtain a split of the tensor: where := and := .
Replacing the SVD by a truncated SVD in this process yields an approximate split of the tensor.The key ingredient of the TRG algorithm is the replacement of all the tensors by their approximate splits.More precisely, for a square lattice, the algorithm proceeds in two steps: First, the initial tensors are split via truncated singular value decompositions along two orthogonal directions: Second, four tensors resulting from the splittings are contracted together to define a new four-valent tensor: This results in a new homogeneous tensor network whose underlying graph is still a square lattice, but rotated by 45 • : Repeating the same two steps would yield a square lattice in the original orientation, with an effective lattice spacing twice as large as the original one.Since the truncated SVDs prevent the bond dimensions from growing uncontrollably high, this coarse-graining scheme can be iterated, hence yielding a flow in the space of tensors from UV to IR.

IV.C. Conceptual shortcomings
The TRG coarse-graining scheme is undoubtedly simple and has proven very successful to study two-dimensional classical systems and ground or thermal states of onedimensional quantum Hamiltonians. 1 Despite its success, the conceptual shortcomings of this simple approach have been known for several years.Indeed, it was emphasized in [28] that, in spite of its name, TRG does not yield a proper RG flow.This is seen easily by inspecting the fixed point structure of the algorithm.Indeed, a simple analysis concludes that the fixed points of the flow generated by TRG depend on the initial tensor, meaning that they display features that are non-universal.This goes against the spirit of RG which states that upon coarse-graining, the UV details below the scale considered should be integrated away.Translating this statement into our language, this signifies that all information contained within the tensors describing short-range physics at a given scale should be contracted during a coarse-graining step.But, perhaps counter-intuitively, the TRG algorithm contracts away only a fraction of short-range information.This is typically illustrated by means of a toy model for purely short-range physics known as the corner-doubleline (CDL) tensor.Given an arbitrary matrix M , a fourvalent CDL tensor is of the form When organized in a square lattice, these CDL tensors give rise to so-called loops of short-range correlations within each square plaquette.These loops indicate that any observable supported on a given leg of the network would be highly correlated with any other observable that has support around one of the plaquettes that include this leg.Conversely, these loops are completely uncorrelated with anything beyond their immediate surrounding.Upon a coarse-graining step of the TRG algorithm, half of the loops are traced away when blocking the tensors resulting from the splits.However, the other half remain untouched and become purely short-range information with respect to the new length scale: Consequently, although these CDL tensors are only nontrivial at the length scale of the lattice, they form fixed points of the flow yielded by TRG.
More generally, all the short-range physics details that are not removed at a given coarse-graining step of TRG accumulate, polluting the coarse-grained tensor and therefore the RG flow.As we already mentioned, this implies that the fixed point structure of the flow is not physically correct, making the identification of the genuine IR fixed points more difficult.Quantitatively, flowing irrelevant UV information in the tensor at each step reduces the space available to store the relevant IR information, and thus the accuracy of predictions at fixed bond dimension.
Crucially, these difficulties become particularly pregnant when studying critical points, making the scale invariance near impossible to witness, and the measurement of critical exponents less reliable [31].Note that by allowing for higher bond dimensions, the quantitative predictions of TRG can be naturally improved.However, because of its failure to contract away all the UV details at every step, it cannot yield a proper RG flow no matter the bond dimension.In many cases, these conceptual shortcomings do not drastically undermine the usefulness of this algorithm, but if we are interested in the RG flow, then it is crucial to use an algorithm that addresses them explicitly.

IV.D. Graph independent local truncations
In order to remove the short-range information that survive the TRG coarse-graining, more truncations need to be made by means of local replacements.But these additional local replacements should be such that the geometry and homogeneity of the network are preserved at every step, so that a flow in the space of tensors can still be defined.Over the years, several schemes that address the shortcomings of TRG and yield proper RG flows have been constructed, e.g.Tensor Network Renormalization TNR [28], Loop-TNR [29] or TNR+ [30].In this manuscript, we use a more recent approach based on an algorithm referred to as Gilt, which performs graph independent local truncations [31].
Gilt is a general-purpose algorithm that performs local truncations of bond dimensions in any tensor network without affecting its geometry.We use it to truncate the bond dimensions of carefully chosen individual legs of the network, thereby removing the short-range information present in some environment surrounding them.
Let us consider a plaquette of the network and pick a leg that is included in this plaquette.Let χ be the bond dimension of this leg.The Gilt algorithm consists in replacing the leg, which can be thought as being the identity matrix, by a non-trivial matrix.This non-trivial matrix is designed to be as low-rank as possible.That way, by splitting the matrix and contracting the components of this split with the neighbouring tensors, we effectively truncate the bond dimension of the original leg.This process can be summarized as follows: with χ < χ.Let us now explain how to design such low-rank matrix.We shall merely provide the solution and refer the reader to [31] for additional details and motivation.
Let us consider the network obtained by removing the leg whose bond dimension we wish to truncate, and perform the SVD illustrated below: The first identification indicates that we consider the blocking of the four tensors as a single tensor.We then perform the SVD of the matrix obtained by grouping the two middle upper most legs together, providing one index of the matrix, and all the other legs together, providing the other index.Notice that we do not truncate any bond dimension at this point.The tensor U provided by the SVD may be used to form a resolution of the identity, which we insert on the original leg of the network In the equation above, we introduced a special notation for the partial trace of U such that t i = tr U i .We are now ready to define the Gilt-matrix as where the vector t is defined in terms of t, the singular values of the SVD (20), and a parameter Gilt as which effectively provides a soft truncation of the singular values of S. This solution was found in [31] to be very good, although not necessarily optimal, to minimize the rank of the Gilt-matrix while maximizing the overlap between the original network and the one after inserting the Gilt-matrix.The hyper-parameter Gilt is chosen by the user depending on the variables of the model and the parameters of the simulation, e.g. the maximal bond dimension authorized.Intuitively, Gilt controls the trade-off between rank reduction and global accuracy.A large Gilt is likely to substantially reduce the rank of the Gilt-matrix and remove more irrelevant short-range information.However, it might introduce a larger local error.Conversely, a small Gilt is less likely to reduce the rank of the Gilt-matrix, leaving us in a situation akin to that of TRG, where irrelevant information propagates from a scale to the next, ultimately taking up all the parameters of our finite χ tensor.In app.B, we explain our procedure to choose Gilt .

IV.E. Renormalization scheme
Recall that the failure of TRG to yield a proper RG flow stems from its inability to remove the short-range information within half of the plaquettes.It was shown in [31] that by applying the Gilt algorithm along the edges surrounding the plaquettes at issue, the short-range information that survive the coarse-graining can be truncated away.More specifically, given such a plaquette, we iterate the Gilt algorithm on all the legs surrounding it, effectively removing the loop of correlations contained within: Notice that the Gilt matrices generally differ from one leg to another.Nevertheless, since the network is homogeneous and we only need to apply the Gilt algorithm on every other plaquette, only four matrices need to be computed throughout.By preceding every iteration of TRG with an iteration of Gilt, we obtain the renormalization scheme referred to Gilt-TNR.First, Gilt matrices are inserted on the legs surrounding every other plaquette, which are then split via truncated SVDs: Blocking four components of the splits around every tensor of the original network, we obtain two different tensors: At this point, we are in the initial configuration of (18) with the difference that we distinguish two types of tensors and half the loop of correlations have been removed thanks to Gilt.TRG can now be applied svd − − → .
It remains to block sets of four tensors as follows ≡ , and in the process remove the last loops of short-range information.The resulting network is still homogeneous so that the whole algorithm can be iterated until the desired infrared length scale is reached.FIG. 1. RG flow in the space of tensors for λ = 0.1 at five different values of µ 2 0 using TRG (upper panel) and Gilt-TNR (lower panel).For every linear system size (or equivalently number of RG steps), the 45 largest singular values of the coarse-grained tensors computed according to (24) are displayed.Each line corresponds to the flow of a given singular value.For both TRG and Gilt-TNR, the maximal bond dimension is χ = 110, and the threshold parameter was chosen to be Gilt = 5 • 10 −8 for Gilt-TNR.Whereas TRG yields non-universal fixed points, Gilt-TNR provides the right fixed point structure away from criticality and reveals a non-trivial fixed point at criticality that is characteristic of the underlying CFT.

IV.F. Extracting observables
As mentioned earlier, we are interested in the RG flow in the space of tensors generated by the Gilt-TNR algorithm.Concretely, the flow is analysed by computing the spectrum of the tensor after every iteration of the algorithm, where by spectrum we mean the largest singular values of the following decomposition: This spectrum provides a precise and basis independent characterization of the tensor.In particular, it allows for the immediate identification of the trivial fixed points of the RG flow of φ 4 2 , characterized by the number of nonvanishing degenerate singular values.Once we have identified the trivial fixed points, we can proceed by dichotomy to find the critical points at which phase transitions take place.This method is simple, robust and its accuracy only depends on the resolution of the input parameters.
As a benchmark, it is also useful to compute the scaling dimensions of the conformal theory at criticality.These observables are also straightforwardly obtained by diagonalizing the transfer matrix resulting from the contraction of two coarse-grained tensors as follows: which is related via a logarithmic conformal map to a discrete version of the dilation operator [44].

Sec. V | Results
In this section we present the results of applying the Gilt-TNR algorithm to φ 4 2 .

V.A. Renormalization group flow
Let us first analyse the RG flow of the lattice theory for a given λ that is not too small.In fig. 1, we present the RG flows generated by the Gilt-TNR and TRG algorithms for λ = 0.1.We display on these graphs the spectra as 2. RG flows in the space of tensors at criticality for four different values of λ using the Gilt-TNR algorithm at bond dimension χ = 110.After a number of iterations that grows as λ gets close to zero, the same non-trivial fixed point, which is characteristic of the underlying CFT, is observed.On the right-hand panel is depicted a heuristic picture of the RG flow, which illustrates the fact that as we go to the continuum limit, the QFT UV fixed point, with a continuum of eigenvalues, starts being visible.
defined in the previous section for coarse-grained tensors corresponding to different linear system sizes.Five pairs of panels are presented corresponding to five different values of µ 2 0 : the one at criticality, and two on both sides of the critical point.
We notice immediately the difference in fixed point structures.First of all, the Gilt-TNR algorithm yields the right trivial fixed points regardless of the initial parameters.These are characterized by a single or a pair of dominant singular values, that correspond to either the disordered phase or the symmetry broken phase, respectively.The TRG algorithm yields non-universal fixed points that depend on the initial parameters.While the non-universal features that correspond to IR-irrelevant information are seen to accumulate as we get closer to the critical point, they are systematically removed by Gilt, eventually leading to a collapse of the singular values, producing the trivial fixed points.
Strikingly, Gilt-TNR yields a non-trivial fixed point at critically that is characteristic of the underlying CFT.Indeed, we observe that as the system size grows, the spectrum remains approximately identical for many iterations, which is the hallmark of scale invariance.This cannot be observed using TRG because of the accumulation of short-range information that requires an unmanageable growth of the bond dimension.
When going to the continuum limit by sending λ to zero, the RG picture becomes more interesting.We naturally have the same IR phases and fixed points, but a new UV fixed point -corresponding to the short distance regime of the continuum QFT -appears.Indeed, we now start from an "ultra-UV" lattice theory, then flow close to the QFT UV fixed point (the massless boson) as irrelevant lattice artifacts decay, and then only reach the IR regime of the QFT, which is independent of the continuum limit.The lattice artifacts at the start of the RG are the same no matter the initial value of λ, but only with very small values of λ can the system flow close to the UV QFT before the relevant perturbations (mass and coupling) drive it away to the infrared (see fig. 2).The traces of the UV QFT are clear in the tensor spectrum.As λ goes to zero, we observe that the initial singular values are systematically larger.We then observe a decay of the first singular values that is more pronounced as λ is small, sign that we are near a fixed point with a continuum of eigenvalues, which is the case for the free boson.As we argue in sec.VI, this makes going to the continuum limit expensive.

V.B. Critical coupling in the continuum
Analysing the renormalization group flow of the model as described above, we can determine the critical coupling f c (λ) of the lattice theory for any values of λ.Since we are interested in the continuum field theory, we probe small values of λ, going as far as λ = 0.0025.We then extrapolate to the continuum limit.The values of f c (λ) for finite λ are shown in fig. 3.
The precision of our data points leaves little doubt    that the leading term in the fit to the continuum limit is not linear but ∝ λ log λ.More precisely, we fit the three following functions to extrapolate to the continuum limit: where χ 2 /d.o.f. is the "variance of unit weight" of the weighted least square fit3 .This gives = 11.0915(32) .
A rigorous analysis of the errors of our algorithm is difficult.Indeed, although we know the errors associated with every approximate replacement performed throughout the algorithm, we are unaware of any methods to accurately estimate the global error from them.As we explain in app.B, we crudely estimate the error bars from the fluctuations of our results as a function of the bond dimension.We notice that the uncertainty associated with the choice of fit is greater than the uncertainty within each fit.To be as conservative as possible, we average the final result over the three fits and quote as error the distance between the two extreme fits majored of their associated errors.We thus quote f cont.c = 11.0861(90), a value which favourably compares with previous studies (see table I).
The existence of a λ log λ has been an open question in the existing literature, with the most recent studies unable to know for sure if a logarithmic term was needed with the precision available [15].We do not know of a theoretical argument for this scaling for the particular case of φ 4  2 , but our results strongly suggest that it should be included.Note that the presence of a logarithmic term explains the slight underestimation of the continuum critical coupling reported in most previous studies that did not consider it.

V.C. Scaling dimensions
To illustrate the efficiency of the algorithm, it is interesting to present benchmark results for the scaling dimensions at criticality, which are known from conformal field theory.Recall that on the whole (λ, µ c (λ))-line, the critical point is in the Ising universality class, and that the Ising model corresponds to λ → ∞.We report in table II, the values we obtain for the first generations of scaling dimensions using the transfer matrix method reviewed briefly in sec.IV F. We do so for different values of λ to showcase the universality of the critical line.Note that the best accuracy is usually obtained for the coarse-grained tensors corresponding to the length scale where scale invariance starts being visible when inspecting the spectrum aforementioned, e.g. for a linear system size of 2 11 for λ = 0.1.We notice that the precision does not depend strongly on the value of λ, at odds with what we see for the critical coupling, where errors increase as we get closer to the continuum limit.This is because the critical exponents are universal and thus insensitive to microscopic details.This is where the fact that we do a proper RG is fundamental: we reach the correct IR fixed point even as our UV errors increase.

Sec. VI | Discussion
We conclude this manuscript with a discussion of the errors and some comments regarding generalization to higher dimensions.

VI.A. Understanding the sources of errors
We have shown that using tensor networks can turn the renormalization group from a qualitative idea to a powerful computational tool.Nonetheless, if our results in the continuum limit are precise, they do no yet reach machine precision.This brings the question of what, if anything, limits the accuracy now and perhaps in a future with larger computing resources.
First, for a fixed spacing a, tensor networks methods are extremely efficient.Although the evidence is not fully conclusive, it is believed that for the type of systems we consider, the Gilt-TNR algorithm makes errors that decrease exponentially in the maximum bond dimension χ [31].Other tensor network methods have a provable exponentially convergent error for local observables [45,46], making the previous observation believable.The cost of the algorithm scales in O(χ 7 ) in d = 2 [31], which is stiff yet still polynomial.Hence, away from the continuum limit, the error can likely be arbitrarily reduced at a moderate cost (see app.B for brief discussion of empirical error estimation).
Going to the continuum limit brings two additional sources of errors as compared to the setup with fixed lattice spacing: (i) The bond dimension of the initial tensor grows linearly in the inverse lattice spacing a for a fixed error, (ii) the RG has to be run for more steps as a → 0 to reach the same IR scale of the field theory.
One could think that (i) is a mere artifact of our discretization, which works well for small values of the field and strong coupling.But this is in fact independent of the basis chosen.Indeed, as we refine the lattice, the theory becomes free, and thus its spectrum becomes that of the free boson, which is continuous.Hence, as we get closer to the continuum limit, the singular values of the tensor decay slower and slower, and more and more of them have to be kept for a fixed error tolerance.
Difficulty (ii) is less of a problem for our method because the scale we consider grows exponentially in the number of RG steps, and thus reducing the lattice spacing only incurs a logarithmic increase in the computational cost.Even for a very fine grained lattice, finite size effects (or IR cut-off effects) can thus be made negligible.This would be much less favorable for methods like the corner transfer matrix renormalization group (CTMRG) [47][48][49], where the scale grows only linearly in the number of iterations.
For the smallest values of λ we probed in this paper, we have reached the point where UV complications, related to driving the system towards the free boson continuum, start to outweigh the IR difficulties, related to the running of the RG flow near the Ising CFT fixed point.In particular, for λ = 0.0025, the initial tensor is required to have a bond dimension that is significantly larger than the maximal one we allow in subsequent steps, in order to make sure that the field discretization error remains negligible.
While we do not believe there is any unacceptable bottleneck with current techniques, new methods to deal with the UV will likely need to be found to make the algorithm more efficient.One could try to use better discretizations of the gradient square term in the action so as to reduce lattice artifacts [50].But the locality of the tensor network structure makes the introduction of next-nearest-neighbor terms less straightforward than in Monte-Carlo.It could also be useful to have a better theoretical control of the expansions of physical quantities away from the continuum as a function of λ.If these were well understood enough, in particular regarding peculiar logarithmic corrections, safer extrapolations to the continuum could be made without the need for extreme fine-graining.
Ironically, the regime we struggle to capture is one that is trivial for perturbation theory, because it is essentially free.We speculate that there must be a way to use the best of the two approaches: Running the flow perturbatively for a while, until the theory is weakly coupled but not infinitely weakly coupled, and then run the hard stronglycoupled part of the flow with our method.We note that the recent bosonic-TRG of Campos et al. [51], which uses continuous variables and Gaussian integrations to renormalize free field theories, seems that it could deal precisely with the regime that is hard to capture with our method.However, how to hybridize the two approaches is not obvious to us at the moment.
Let us finally mention that as we increase the number of physical dimensions to 3+1 and go to QFTs that are no longer super-renormalizable but just renormalizable, the previous UV difficulties will become milder.Indeed, even for an asymptotically free theory like quantum chromodynamics, the coupling constant is believed to run to zero only logarithmically in the lattice spacing (instead of polynomially in super-renormalizable ones).This means that we can start from a lattice spacing dozens of orders of magnitude smaller than the physically relevant length scales while keeping a theory that is not too weakly coupled and thus does not require a tensor with a pathologically large number of singular values.However, as we will see, increasing the number of dimensions brings other algorithmic difficulties.

VI.B. Going beyond two dimensions
In the future, the challenge is not so much to improve the precision in (1+1)d than to go to higher dimensions.It seems natural to tackle φ 43 next, to increase difficulty progressively.Indeed, it is still a well defined superrenormalizable QFT, although its construction is slightly less trivial.The IR critical point is believed to also belong to the Ising universality class, and thus its critical exponents, which are useful for sanity checks, are known to very good precision thanks to impressive progress in the conformal boostrap program [52][53][54][55].
Increasing the number of dimensions does not present any major conceptual difficulties for the lattice Monte-Carlo approach, but all non-perturbative deterministic approaches suffer in one way or another.We compare their prospects with those of our own approach.
Renormalized Hamiltonian truncation methods are promising but encounter a number of difficulties.First, in (2+1)d, the λφ 4 coupling term is less relevant, which makes the method a priori converge slower.Second, it is no longer the case that the continuum Hamiltonian can be defined straightforwardly with normal ordering, and additional counter terms are needed in (2+1)d, requiring a non-trivial regularization (see however the recent advance [56]).Finally, the number of momenta to consider for a fixed energy cutoff is squared from (1+1)d to (2+1)d.
Matrix product state are easily generalized to projected entangled pair states (PEPS) in (2+1)d.The latter have been used to solve non-trivial problems to impressive precision [57][58][59][60] but on the lattice and sufficiently away from criticality.The method is variational and limited to low bond dimensions beyond which the convergence is too slow.Near the continuum limit, the correlation length becomes infinite and we expect that larger physical and bond dimensions will be required, making PEPS less accurate.
Generalizing our method to higher dimensions is also not straightforward.In (2+1)d, this would amount to contracting efficiently a cubic tensor network.To do so, we would proceed as in (1+1)d, by combining local truncations that remove short-range correlations with a coarse-graining procedure.Although the Gilt algorithm is also applicable for three-dimensional networks, 4 the corresponding numerical cost grows stiffly from O(χ 7 ) to O(χ 12 ).Indeed, the Gilt matrices must now be computed with respect to a cubic environment so as to remove threedimensional short-range correlations.Due to this rather challenging cost, the performance of such an algorithm has not been conclusively demonstrated yet.
In the end, for strongly coupled QFT in (2+1)d and a fortiori (3+1)d, no deterministic method has a demonstrable edge so far.However, we believe that tensor network renormalization approaches, which are accurate for universal and non universal quantities in (1+1)d and suffer from no crippling difficulty in (2+1)d, are a promising option that deserves to be explored further.

CD would like to thank Markus Hauru for countless discussions regarding the renormalization of tensor networks, collaboration on closely related topics, as well as useful comments about this manuscript. This project has received funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme through the ERC Starting Grant WASCOSYS (No. 636201), as well as the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) under Germany's Excellence Strategy -EXC-2111 -390814868.
App.A | More on the field discretization Let us bound the coefficients of T abcd .First, because log(n!) is convex, we have that where n = (a + b + c + d)/2, so that T is dominated by its diagonal.Next, so long as µ 2 0 > −4 (which is always verified near the continuum limit where µ 0 → 0), h is easily bounded: which, e.g. from Stirling's formula, goes to zero faster than exponentially.This decay is only asymptotic and for λ small enough, the coefficients of T first grow exponentially, before they are killed off by the factorial.The smaller λ, that is the closer to the continuum limit we are, the later this tipping point occurs.This crossover can be estimated crudely with Stirling's formula: which is negative when n ≥ 8/λ.Thus, we see that K max has to grow only linearly in λ −1 to insure that we are in the right asymptotic regime where the truncation error decreases super-exponentially.
As we go to the continuum limit and λ → 0, we can get an estimate of the relevant range K min , K max of indices to consider.It is much narrower than the naive 0, K max , because K min also grows linearly with λ −1 for a fixed error threshold.Indeed, up to a normalization and for even integers, the right-hand side of (A2) is the Poisson distribution of parameter 4/λ.The latter converges to a Gaussian with standard deviation σ ∝ λ −1/2 when λ → 0. Hence, as we go to the continuum limit and for a fixed error threshold, K = K max − K min grows at most as λ −1/2 .Equivalently, the bond dimension of the initial tensor grows linearly with the inverse lattice spacing a −1 .

App. B | Convergence in bond dimension χ and choice of Gilt
As mentioned in the main text, the error bars quoted in fig. 3 are obtained from the fluctuations of our results as a function of the bond dimension.As we increase the bond dimension, we numerically observe that our results converge with apparent exponentially decaying oscillations around a limiting value.Since we have no theoretical means to estimate the error generated throughout the simulation, we use the amplitude of these oscillations as a ballpark estimate, which is most likely an overestimation of the error induced by the algorithm when using the highest bond dimension, namely χ = 110.More precisely, we compute the critical coupling for χ = 90, 95, 100, 105, 110 and use the standard deviation of these five points as error bar, which we center around the value at χ = 110.It turns out that the explanation provided above is slightly more subtle in practice, due to the fact that the bond dimension is not the only parameter of our algorithm.Indeed, the Gilt component of our renormalization scheme is tuned via the threshold parameter Gilt , which we recall is responsible for the trade-off between global accuracy and local truncation.Although there are no rigorous theoretical arguments to decide which value this parameter should be set to, we have several concordant heuristics.
From the definition (23) of the Gilt matrix, we understand that if Gilt is too small, then Gilt will not truncate sufficiently and UV information will creep into the IR, and if it is too large, then the local error will be too important.By simply monitoring the amount of truncation performed by Gilt, we can already determine that choosing Gilt in the range 10 −6 to 10 −8 for bond dimensions in the range [90, 110] is suitable.Furthermore, Gilt should be chosen such that the error caused by the algorithm remains dominated by the coarse-graining step.Otherwise, increasing the bond dimension would have a very limited effect.In practice, this signifies that there is no optimal choice of Gilt given a wide range of χ, and that as χ is increased, Gilt should be decreased accordingly.
Our last heuristic is the expectation that, given a fixed bond dimension, small variations around the optimal value of Gilt should have a negligible effect on the renormalization flow.In other words, we expect to find a range of values for Gilt that yield the same final result.Naturally, the span of this range depends on the bond dimension we use.We do observe this phenomenon in practice and further notice that, close to the optimal value, a smaller Gilt tends to favour the symmetry-broken phase, whereas a larger Gilt tends to favour disordered phase.As such, the optimal Gilt is typically found at a stationary point.
Concretely, we applied this program to determine the optimal Gilt for bond dimensions χ = 90, 95, 100, 105, 110.In order to simplify the presentation, we made the choice of using the same Gilt (χ) for every λ.Although this might not be optimal for every λ, we ensured that these values perform well for the whole range we consider.We found that Gilt (90) = 10 • 10 Finally, let us point out that this whole procedure is very specific to the algorithm we use.Indeed, by choosing Gilt , we only decide indirectly the bond dimension to which the bonds are truncated to during the Gilt step.Although this requires the preliminary study presented in this appendix, this approach has the advantage of being particularly versatile.Indeed, the bond dimension adapts automatically according to whether truncation can be done with more or less error, which depends on the amount of short-range information available.This turns out to be particularly relevant for our study as the initial tensors for the smallest values of λ typically have a bond dimension that is larger than the maximal one that can be numerically dealt with during the coarse-graining step.We observe that the first step of Gilt significantly reduces the bond dimension and allows the algorithm to proceed.

FIG. 3 .
FIG. 3. Critical coupling fc = λ/µ 2 c as a function of λ, computed for χ = 110 and Gilt = 5 • 10 −8 .The error bars are estimated as described in app.B. The first fit up to linear order does not include the point λ = 0.1 as it is expected to be valid only for smaller values of the coupling.The right-hand panel shows a narrower view of the same figure close to λ = 0.

TABLE I .
Comparison of several estimates of the critical coupling constant f cont.