{'fraction_non_alphanumeric': 0.039265842907724544, 'fraction_numerical': 0.009668328636462171, 'mean_word_length': 4.3844836977667105, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 4, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 9, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Quantum Criticality with the Multi-scale Entanglement Renormalization Ansatz

The goal of this manuscript is to provide an introduction to the multi-scale entanglement renormalization ansatz (MERA) and its application to the study of quantum critical systems. Only systems in one spatial dimension are considered. The MERA, in its scale-invariant form, is seen to offer direct numerical access to the scale-invariant operators of a critical theory. As a result, given a critical Hamiltonian on the lattice, the scale-invariant MERA can be used to characterize the underlying conformal field theory. The performance of the MERA is benchmarked for several critical quantum spin chains, namely Ising, Potts, XX and (modified) Heisenberg models, and an insightful comparison with results obtained using a matrix product state is made. The extraction of accurate conformal data, such as scaling dimensions and operator product expansion coefficients of both local and non-local primary fields, is also illustrated.

Tensor network states are powerful variational ansätze that can be used to characterize the low-energy properties of quantum many-body systems on a lattice. The premise of tensor network approaches is to parameterize a many-body wave-function by using a collection of tensors connected into a network. The number of parameters required to specify these tensors is much smaller than the exponentially large dimension of the system's Hilbert space, in such a way that very large (and even infinite) lattices can be considered.
Tensor network states can be broadly classified into two sets according to the geometry of the underlying networks [1]. In the first set, the network reproduces the physical geometry of the system, as specified by the pattern of interactions in the Hamiltonian. For instance, the matrix product state (MPS) [2][3][4], an ansatz for D = 1 dimensional systems, consists of a collection of tensors connected into a chain; similarly, its generalization for lattices in D > 1 dimensions, known as projected entangled pair states (PEPS) [5][6][7], consists of a collection of tensors connected according to a D dimensional lattice. In contrast, a second class of tensor network states aim at reproducing the holographic geometry of the system. The latter spans an additional dimension used to parameterize the different length scales (or, equivalently, energy scales) relevant to the description of the many-body wave-function. Thus the multi-scale entanglement renormalization ansatz (MERA) [8][9][10][11][12][13][14][15][16][17][18][19][20][21][22][23] for a lattice system in D dimensions consists of a network in D + 1 dimensions.
The simplest and most widely employed tensor network state is the MPS. The MPS underlies the remarkable success of the density matrix renormalization group (DMRG) algorithm [24][25][26][27], which for almost two decades has dominated the numerical study of D = 1 dimensional quantum systems, providing very accurate results for low energy properties. DMRG has not only become ubiquitous in condensed matter physics but has also found application in other fields involving quantum many-body systems, such as quantum chemistry [28]. Further algorithms based upon the MPS have also been developed, such as the time evolving block decimation (TEBD) [29,30] algorithm, later reformulated as time-dependent DMRG [31,32], which allows the simulation of certain low-energy dynamics for relatively long times.
In this manuscript we discuss the application of the MERA to study critical systems in D = 1 dimensions (although most of the present formalism is directly applicable to D > 1 dimensions). Given the success of MPS-based methods such as DMRG and TEBD for quantum systems in D = 1 dimensions, it is natural to ask whether an alternative approach is actually needed. A clear answer to this question is obtained by discussing the short-comings of the MPS representation for critical systems as well as by exploring the benefits of including scale invariance directly into a tensor network state, something that is possible with the MERA but not the MPS.
Critical systems typically lack a characteristic length scale and are thus invariant under changes of scale. One manifestation of scale invariance in critical systems is that correlators decay polynomially, in sharp contrast with gapped systems, where they decay exponentially with a characteristic correlation length. It is well-known, however, that an MPS with a finite bond dimension χ (where χ indicates the size of the tensors) can never properly capture the scale invariance of a critical state. Indeed, a duly optimized MPS possesses an intrinsic finite correlation length ζ ≈ χ κ [33,34], where κ is a constant that depends on the universality class of the phase transition under consideration, such that correlators decay exponentially at length scales larger than ζ. Thus, while the MPS can accurately approximate short-range properties of a critical ground state, it necessarily fails to capture its properties at asymptotically large distances. [In practice, however, the cost of MPS-based approaches scales only as O(χ 3 ) with the bond dimension χ. This means that one can use a very large value of χ, which often allows the critical behavior of a system to be accurately captured up to some very large length scale ζ.] On the other hand, the MERA can explicitly capture the scale invariance of critical systems [8, 10, 12, 15-17, 19, 22], a feature that has significant advantages, both conceptual and practical. Tensors in the MERA are organized in layers, where each layer corresponds to a different length (or energy) scale. In an infinite system, scale invariance is then easily enforced by choosing all layers of tensors to be identical. The resulting ansatz is referred to as the scale-invariant MERA. Certain structural properties of the scale-invariant MERA, such as the polynomial decay of correlators and the logarithmic growth of block entanglement entropy [1,9], already hint at its suitability to represent critical systems.
In addition, this ansatz offers direct access to the scaling operators of a critical theory, namely those operators that transform into themselves under scale transformations. As a scale-invariant/covariant object, a scaling operator must act non-trivially on a region of the system that has no characteristic length scale. In a (1+1)-dimensional conformal field theory (CFT) [35][36][37] (corresponding to the continuum limit of a critical quantum system in D = 1 spatial dimensions), the support of a scaling operator can therefore only be one of three possibilities: (i) an infinite line, (ii) a single point, or (iii) a semi-infinite line. The first type of support corresponds to a global internal symmetry of the CFT's Hamiltonian. The second type of support is seen to correspond to local scaling operators, associated to local excitations. Finally, the third type corresponds to non-local (or semi-local) scaling operators, associated e.g. to domain wall excitations. Going back to the lattice, scaling operators are distorted by the presence of a finite lattice spacing, but they can still be directly extracted from the scale-invariant MERA. Thus, on the lattice, (i) a global internal symmetry is implemented by an infinite string of identical single-site operators; (ii) local scaling operators are now supported on a small number of sites (the specific number depends on the MERA scheme); and (iii) non-local operators mix elements of the two previous objects: they consist of a semi-infinite string of identical single-site operators (the same ones that implement an internal symmetry) completed with some local operator at the end of the string. Importantly, the scaling dimensions and fusion rules of the scaling operators on the lattice are the same as in the continuum. As a result, a relatively simple and inexpensive calculation with the scale-invariant MERA can be used to obtain remarkably accurate estimates of the conformal data characterizing the underlying CFT.
The rest of the manuscript is organized in sections as follows. Sect. II introduces the key aspects of entanglement renormalization, which is the renormalization group transformation for lattice models on which the MERA is based. We describe how local operators transform under the coarse-graining transformation and discuss basic aspects of the MERA, including the evaluation of the expectation value of local observables, and briefly compare different MERA implementations in D = 1 dimensions. Sect. III addresses the role of spatial and internal symmetries in tensor network states, and compares how different symmetries can be enforced in MPS/PEPS and MERA. While translation invariance is naturally implemented in an MPS and PEPS, it can only be approximately enforced on the MERA. Sect. IV specializes on the implementation of scale invariance in the MERA and discusses how the scaling operators of a critical theory can be extracted from it. In Sect. V we demonstrate the performance of the scale-invariant MERA for a number of critical quantum spin chains, namely Ising, Potts, quantum XX and Heisenberg models. Specifically, Sect. V A compares ground state energy and two-point correlators obtained with MERA and MPS. Interestingly, MPS and MERA approaches seem to complement each other. For a given computational cost, the MPS is seen to provide more accurate estimates of the expectation value of a local observable, such as the ground state energy. However, the MERA is seen to provide a better characterization of long-range properties, such as two-point correlators at long distances. The advantages of the scale-invariant MERA are then further illustrated in Sect. V B by extracting, in the concrete context of the quantum Ising model, the conformal data (including scaling dimensions, fusion rules for scaling operators and central charge) of the underlying CFT.  1. (i) The coarse-graining transformation U , a specific implementation of entanglement renormalization, is comprised of isometries w and disentanglers u and maps blocks of three sites from the initial lattice L into a site of the coarser lattice L . (ii) The tensors w and u are constrained to be isometric, see also Eq. 3.

II. ENTANGLEMENT RENORMALIZATION AND THE MERA
In this section we first recall a few basic aspects of entanglement renormalization (ER), the coarse-graining transformation upon which the MERA is based. We then introduce the MERA and review a few of its features.

A. Foundations of Entanglement Renormalization
For concreteness, we mostly consider a specific implementation of ER that produces the so-called ternary MERA (where three sites are coarse-grained into one effective site). In Sect. II C we also discuss other MERA schemes.
Let L denote a D = 1 dimensional lattice made of N sites, where each site is described by a Hilbert space V of finite dimension d, so that the vector space of the lattice is V L ∼ = V ⊗N . We consider a coarse-graining transformation U that maps lattice L to a coarser lattice L where L is made of N/3 sites, each with a vector space V of dimension χ, so that V L ∼ = V ⊗N , and where transformation U decomposes into local transformations, known as disentanglers u and isometries w, according to Fig. 1(i). More specifically, if we partition the initial lattice L into blocks of three sites, then the disentanglers u are first applied across the boundaries of neighboring blocks, followed by the isometries w, which map each block of three sites into a single effective site of the coarser lattice L . Disentanglers and isometries are required to satisfy isometric constraints, namely where I and I are the identity operator in V and V , respectively, see Fig. 1(ii). Note that, by construction, the disentangler u is also unitary, that is uu † = I ⊗2 . The dimension χ of the Hilbert space V can be chosen to be different than d, provided that χ ≤ d 3 (as demanded by the above isometric constraint on w). In general, choosing a larger dimension χ, i.e. retaining a larger effective Hilbert space V for each coarsegrained site, yields a more accurate RG transformation, one that better preserves the low energy properties of the system. An important property of the coarse-graining transformation U is that it preserves locality. Let o(r 1 , r 2 ) be a local operator defined on two contiguous sites (r 1 , r 2 ) of lattice L. Under coarse-graining, the operator o(r 1 , r 2 ) becomes where the only non-trivial part of U † o(r 1 , r 2 )U is a new operator o (r 1 , r 2 ) supported on two contiguous sites (r 1 , r 2 ) of lattice L . Notice that operator o (r 1 , r 2 ) remains local (i.e., it is supported on two sites) thanks to both the specific decomposition of U into disentanglers and isometries, and the isometric constraints of these tensors, Eq. 3, which ensure that most of the tensors of U in U † o(r 1 , r 2 )U annihilate to identity with their conjugates in U † , as shown Fig. 2(i). In view of this fact, it is most convenient to introduce left, center and right ascending superoperators {A L , A C , A R }, as shown Fig.  2(ii), which directly produce the two-site coarse-grained operator o from the two-site operator o, as given by one of the following, where the specific choice of ascending superoperator to be used depends on the location of the operator o on the lattice L.
We may now concatenate the coarse-graining transformation a number T of times to obtain a sequence of coarser lattices, where we use superscripts in square brackets to denote the level of coarse-graining, with the initial lattice L [0] ≡ L.
through application of the appropriate ascending superoperator A, as per Eq. 5.

B. Foundations of the MERA
We have just seen that the ER transformation U can be used to coarse-grain local operators, producing a renormalization group (RG) flow for local operators, Eq. 7. As a linear (isometric) map from V L to V L , U can of course also be used to coarsegrain quantum states. More important for us, however, is to consider an inverse RG flow of states. Let us assume that we have the sequence of ER transformations U Let us assume that the number of levels T is chosen T ≈ log 3 (N ), such that the maximally coarse-grained lattice L [T ] contains a small number of sites and hence the state ψ [T ] can also be described with a small number of parameters. Then the multiscale entanglement renormalization ansatz (MERA) is the class of states ψ [0] that can be represented as Eq. 9, for some choice of U [0] , U [1] , . . . , U [T −1] and |ψ [T ] . For instance, Fig. 3 depicts the MERA, organized into T = 3 layers, for a state ψ [0] on a lattice of N = 54 sites, Let us now count variational parameters. We will assume for simplicity that for any value of τ = 0, 1, · · · , T − 1, the dimension of the vector space V [τ ] is χ. Recall that the transformations {U [τ ] } themselves are comprised of local tensors, the disentanglers u and isometries w, each specified by χ 4 parameters. Since in an N -site lattice there are O(N ) disentanglers u and isometries w (distributed in O(log(N )) layers), the MERA depends on O(N χ 4 ) parameters. In Sect. III we will see that, through incorporation of spatial symmetries into the ansatz, this number can be reduced to O(χ 4 ), which is independent of N and allows for the study of infinite systems.
We have therefore established that the MERA can be specified with a number of parameters that is much smaller than the dimension of the Hilbert space V L , which grows exponentially in the number N of sites. However, for this ansatz to be useful, we also need to be able to efficiently extract information about |ψ [0] in Eq. 9. For any local operator o [0] (r 1 , r 2 ), and due to the very peculiar causal structure of its tensor network, it is actually possible to efficiently compute the expectation value from the MERA. Let us define the past causal cone of a site in lattice L [τ ] as the set of tensors and indices that can affect the state on that site. By construction, in a MERA the causal cone of any site of L [τ ] is seen to involve just a constant (that is, independent of N ) number of sites of any other lattice L [τ ] for τ > τ , a property that we refer to by saying that the past causal cone has bounded 'width'. Fig. 3 displays the past causal cone of two sites (r 1 , r 2 ) in a ternary MERA, which only involves two sites of every lattice L [τ ] . This property allows for the efficient computation of local reduced density matrix ρ [0] (r 1 , r 2 ), from which the expectation value can be obtained. The computation of local reduced density matrices is simplified through the introduction of left, center and right descending superoperators, {D L , D C , D R }, which are the adjoints of the ascending superoperators of Eq. 5. Let us assume that we have the density matrix ρ [τ ] (r 1 , r 2 ) describing the state on two contiguous sites (r 1 , r 2 ) of lattice L [τ ] . Then, as shown Fig. 4, the descending superoperators may be used to compute the two-site reduced density matrix ρ [τ −1] on certain sites of the lattice Thus, through repeated use of the appropriate descending superoperator in Eq. 13, we can compute the reduced density matrix ρ [0] (r 1 , r 2 ) of any two contiguous sites (r 1 , r 2 ) of the original lattice L [0] by 'lowering' the density matrix through the appropriate causal cone. For instance, the reduced density matrix ρ [0] (r 1 , r 2 ) on the two sites (r 1 , r 2 ) of lattice L [0] shown in Fig. 3 can be computed as The cost of calculating a local reduced density matrix ρ [0] , as in the example of Eq. 14, is proportional to the number T ≈ log(N ) of layers in the MERA, hence scales with system size N as O(log(N )). Two point correlators can also be evaluated using similar manipulations, see Ref. 15.

C. Choice of MERA Scheme
There are many possible ways of implementing the MERA in D = 1 dimensions, of which the ternary MERA described in Sect. II B is just one example. Fig. 5 displays a ternary MERA together with two other possible implementations: the binary MERA scheme (in terms of which the first ER proposals [8,9] were formulated) and a modified binary MERA scheme (with half the amount of disentanglers as the previous binary scheme). While all MERA schemes function similarly on a conceptual level, the computational efficiency may differ between schemes. For instance, the cost of evaluating the expectation value of local observables as described Sect. II B scales, in terms of the bond dimension χ, as O(χ 9 ) for the binary MERA, as O(χ 8 ) for the ternary MERA and as O(χ 7 ) for the modified binary MERA. On the other hand, the binary MERA scheme has more disentangling power than either the ternary or modified binary schemes and, for any given χ, produces a more accurate representation of ground states. It is therefore not obvious which MERA implementation will give the best numeric results for a fixed computational budget. However, a direct comparison of performance, see Sect. B, shows that the modified binary scheme of Fig. 5(iii) is the most efficient scheme. Consequently, this scheme is used for the obtaining the numeric results presented in Sect. V. However, for the sake of simplicity, we shall continue to discuss theoretical aspects of MERA in terms of the ternary MERA.

III. SYMMETRIES IN TENSOR NETWORK STATES
Consider a many-body state that is invariant under some symmetry transformation. In approximating this state with a tensor network state, we would like to preserve the original symmetry as much as possible. In this section we examine the types of symmetries that can be enforced upon the MPS and PEPS, and upon the MERA. We also examine whether the presence of the symmetry can be exploited for computational gain. We begin by discussing spatial symmetries, followed by global internal symmetries. The results are summarized in Table I.

A. Spatial Symmetries
For simplicity we discuss only two typical spatial symmetries: invariance under translations in homogeneous systems and invariance under changes of scale in e.g. critical systems. Let us first consider them in the context of the MPS for D = 1 dimensions and PEPS for D > 1 dimensions.
In an inhomogeneous MPS/PEPS, one associates a different tensor A r to each site r of the lattice L. Hence, in a lattice made of N sites, the total number of tensors in the tensor network is also N . In the presence of translation invariance (either in a finite system with periodic boundary conditions or in an infinite system), this symmetry can be incorporated into the MPS/PEPS by choosing all the tensors to be a copy of the same tensor A, i.e. A r = A, see Fig. 6(i) for an MPS. Translation invariance is in this way exactly preserved. It can also be exploited to reduce computational cost from On the other hand, it is not clear how scale invariance could be enforced in these tensor networks. For an MPS this is unlikely to be possible at all because, as we mentioned earlier, a finite bond dimension χ already implies the presence of an effective finite correlation length ζ ≈ χ κ [33,34].
Let us now consider spatial symmetries in the MERA. Recall that a generic MERA on an N site lattice is arranged into T ≈ log(N ) layers of tensors, and contains O(N ) different tensors, as depicted in Fig. 6(iia). Suppose now that the state to be approximated by the MERA is translation invariant. Then we can choose all the tensors in each layer to be the same, so that layer U [τ ] is characterized by a single pair of tensors u [τ ] and w [τ ] , see Fig. 6(iib). In this way translation invariance can be exploited to reduce computational costs from O(N ) to O(log(N )). Notice, however, that this choice of tensors does not enforce translation invariance, because the structure of the coarse-graining is not homogeneous (different sites are positioned in inequivalent positions with respect to the disentanglers and isometries). The final effect is examined in Fig. 7. A MERA characterized by a single pair of tensors u [τ ] and w [τ ] for each layer, where these tensors are filled with random coefficients (compatible with the isometric constraints of Eq. 3), is highly non-translation invariant, with e.g. oscillations in the expectation value of the energy of the order of 0.1 for the Hamiltonian H Ising of Eq.  36. Still, these violations of translation invariance decrease significantly once the tensors are optimized so as to minimize the expectation value of the translation invariant Hamiltonian H Ising . Indeed, they become of order 10 −5 for χ = 4 and decrease with increasing χ. [In practice one can efficiently average the expectation value of a local observable over all possible lattice positions in order to further reduce the effect of these small violations of translation invariance]. We conclude that translation invariance can be exploited to reduce computational costs, but it can only be reproduced approximately. It is not known whether it can be enforced exactly.
Instead, enforcing scale invariance in the MERA is straightforward. This is accomplished by choosing all disentanglers and isometries to be copies of a single pair u and w, see Fig. 6(iic), which further reduces the number of parameters and the computational cost of MERA algorithms from O(log(N )) to O(1), allowing infinite systems to be considered. The scale-invariant MERA will be discussed in more detail Sect. IV.
To summarize, in the MPS/PEPS we can enforce and exploit translation invariance but not scale invariance, whereas in the MERA we can enforce and exploit scale invariance but only exploit (i.e., we cannot enforce) translation invariance. Thus both MPS/PEPS and MERA have potential advantages over each other, depending on whether exact translation invariance or exact scale invariance is more important for the problem under consideration. Once the MERA has been optimized for the ground state of HIsing, it more closely approximates translation invariance; bond energies now differ from the average by less than 0.1%. As the bond dimension χ of the MERA is increased, the optimized wavefunction better approximates translation invariance; for χ = 16 the bond energies differ from the average energy by less than 0.001%.

B. Global Internal Symmetries
A second important class of symmetries are those involving internal degrees of freedom, such as Z 2 spin flips and U (1) or SU (2) spin rotations simultaneously applied on all the sites of a spin model. Such symmetries can be enforced and exploited in all tensor networks.
Let us assume that the Hamiltonian H of our lattice model is invariant under a symmetry group G, where Γ g ≡ · · · V g ⊗ V g ⊗ V g · · · is an infinite string of copies of a matrix V g , with V g a unitary representation of G, and let |ψ be the ground state of H, which we will assume to have the same symmetry, i.e. Γ g |ψ = |ψ (or, more generally, Γ g |ψ = e iφ |ψ ). We can then ensure that the symmetry is also exactly preserved in a tensor network approximation to |ψ by using symmetry preserving tensors [38,39]. For instance, for the MERA, we choose the disentanglers u and isometries w such that, where V g acting on different indices may actually denote different (in general, reducible) representations of G, see also Fig. 8. The use of symmetry preserving tensors implies that the tensors are block diagonal when viewed in a certain basis and thus contain less free parameters than generic tensors. This reduction in the number of parameters can be exploited to significantly decrease computational costs. Symmetries, and in particular a truncated version of the operator Γ g , also play an important role in the description of non-local scaling operators, as discussed in Sect. IV D.

IV. SCALE-INVARIANT MERA
We have already introduced the scale-invariant MERA: in a lattice L with an infinite number of sites, N → ∞, it consists of infinitely many layers of tensors, where all the disentanglers and isometries are copies of a unique pair u and w. In this section we enumerate two significant structural properties of the scale-invariant MERA and review how one can compute a local reduced density matrix, from which the expectation value of a local operator can be evaluated. Then we discuss the three types of scale-invariant (or covariant) objects one can extract from it.

A. Basic Properties
Two basic features of the scale-invariant MERA in D = 1 dimensions match wellknown properties of the ground state of a critical system. Firstly, the entanglement entropy S L of a block of L contiguous sites can be seen to scale as the logarithm of L [9], which is compatible with the critical scaling [40,41], where c is the central charge of the CFT. Secondly, correlation functions can be seen to decay polynomially [9], as it is also expected of critical correlators. Interestingly, these two properties of the scale-invariant MERA follow from simple geometric considerations, namely by studying minimally connected regions and geodesic paths in the (discrete) holographic geometry generated by the tensor network [1].

B. Transitional Layers
In a practical computation (see Sect. A) it is customary to consider a scale-invariant MERA with some small number M of translational layers Secondly they also allow to diminish the effect of RG irrelevant terms in the critical Hamiltonian H of the system. Such terms violate scale invariance but become less and less important at larger length scales. The number M of transitional layers required depends on the amplitude and scaling dimensions of the irrelevant operators present in H, and is often determined by trial and error. For the sake of simplicity, in the rest of this section we shall focus on the case of a purely scale-invariant MERA with no transitional layers.

C. Local Density Matrix
The computation of the (average) local density matrix for two contiguous sites of lattice L is of central importance in the present formalism. The density matrixρ is required both in order to extract the expectation value of a local operator o from the scale-invariant MERA and to optimize its tensors so as to approximate the ground state of a critical Hamiltonian H.
Here we consider the evaluation of the expectation value o(r, r + 1) of a local observable o(r, r + 1). As discussed in Sect. III A, the MERA in not manifestly translation invariant. Thus the expectation value o(r, r + 1) can artificially vary with the position of site r in the lattice. To mitigate this effect, rather than evaluating the expectation value at a particular lattice position r we will instead evaluate the average expectation value over all lattice sites, Notice that this average expectation value can be expressed in terms of the average two-site reduced density matrixρ introduced in Eq. 19, In Sect. II B we described the use of the left, center and right descending superoperators, {D L , D C , D R }, to compute the reduced density matrix ρ [0] (r, r + 1) from a finite MERA. In particular, it was argued that obtaining the density matrix ρ [0] (r, r + 1) required application of a specific sequence of left, center and right descending superoperators that depended on the causal cone associated to sites (r, r + 1), see Eq. 14 for an example. The average density matrixρ can be seen to follow from using the average descending superoperatorD, defined as in order to descend the density matrix through the 'average' causal cone. That is, given the average density matrixρ [τ ] at level τ , the average density matrixρ [τ −1] at lower level τ − 1 is obtained asρ In an infinite system, n → ∞, the MERA has T → ∞ layers, and the average density matrixρ is obtained fromρ whereρ is simply the dominant eigenoperator of the descending superoperatorD, which is independent ofρ [T ] . As a manifestation of scale invariance, this is also the two-site density matrix of any coarse-grained lattice L [τ ] , that isρ [τ ] =ρ for any τ ≥ 0. More details on the computation ofρ can be found in Sect. A.

D. Scale-invariant Objects
The scale-invariant MERA offers direct access to objects of a critical theory that are invariant (more generally, covariant) under a change of scale. In D = 1 dimensions we can identify three classes of such objects, as depicted Fig. 10(i-iii). Next we discuss them in some detail.

Symmetry Transformations
Let us assume that the critical ground state is invariant under some symmetry group G, as implemented by the symmetry transformations Γ g ≡ · · · V g ⊗V g ⊗V g · · · introduced in Eq. 15, where g ∈ G. The first type of objects that transform into themselves under changes of scale correspond precisely to the infinite strings Γ g , for which we have Indeed, if we use symmetry preserving tensors in the MERA, as per Eq. 16, then the string Γ g commutes with each layer U of the MERA, Γ g U = U Γ g , or equivalently the string Γ g remains invariant under coarse-graining, U † Γ g U = Γ g , as shown in Fig.  10(i).

Local Scaling Operators
The second class of objects with a simple transformation rule under coarse-graining, which can be easily extracted from the scale-invariant MERA, are local scaling operators φ α , fulfilling where λ α is some constant.
For simplicity, here we focus on one-site scaling operators [One can also compute two-site scaling operators, but they lead to the same scaling dimensions and fusion rules]. As depicted in Fig. 10(ii), a one-site operator o located at certain points on the lattice is coarse-grained into another one-site operator o . This coarse-graining can directly be implemented with the one-site ascending superoperator, which we call the (one-site) scaling superoperator S in the scale-invariant setting, see also Fig. 10(iv). Iteration produces an RG flow for one-site operators: The scaling operators φ α and their corresponding scaling dimensions ∆ α , can then be obtained by simply diagonalizing the scaling superoperator S [14,16].

Non-Local Scaling Operators
Let us assume again that the critical ground state represented by the scale-invariant MERA is invariant under a symmetry group G, as implemented by the symmetry transformations Γ g ≡ · · · V g ⊗ V g ⊗ V g · · · introduced in Eq. 15, where g ∈ G, and that the tensors of the MERA have been chosen to preserve this symmetry, as per Eq. 16. We can then identify a third class of objects with a simple transformation rule under changes of scale, namely non-local scaling operators φ g,α , to be defined below, which fulfill where λ g,α is some constant. To see how these scaling operators come about [22], let us first introduce non-local operators o g of the form, where Γ g is a semi-infinite string made of copies of V g and o is a one-site operator attached to the open end of Γ g . Notice that, under coarse-graining, o g can be mapped into another non-local operator o g of the same type, since the semi-infinite string Γ g commutes with the coarse-graining everywhere except at its open end, as illustrated in Fig. 10(iii). Thus we can study the sequence of coarse-grained non-local operators o g −→ o g −→ o g · · · by just coarse-graining the local operator o with the modified one-site scaling superoperator S g of Fig. 10(v), In particular we can diagonalize the modified scaling superoperator S g , to obtain non-local scaling operators φ g,α of the form Notice that for g = I we recover the local scaling operators φ α of Eq. 29. Importantly, the scaling dimensions ∆ α and ∆ g,α of the both local and non-local scaling operators φ α and φ g,α (as well as their operator product expansion coefficients, see Ref. 16) are the same in the lattice than in the continuum. Therefore by extracting properties of the scaling operators on the lattice, we can characterize the CFT that describes the critical theory in the continuum. As demonstrated by the benchmark results of Sect. V B, a relatively simple and inexpensive MERA simulation can actually be used to obtain remarkably accurate conformal data of the underlying CFT.

V. BENCHMARK RESULTS
In this section we benchmark the performance of the scale-invariant MERA by applying it to study of the ground state of several well-known quantum spin chains at criticality. The models we analyze are the critical Ising model [42,43], the critical three-state Potts models [44], the XX model [45] and a Heisenberg zig-zag chain [46] (the Heisenberg model with a next-nearest neighbor coupling), corresponding to the following Hamiltonians: H Potts = r Z (r) −X(r)X † (r + 1) −X † (r)X(r + 1) H XX = r X(r)X(r + 1) + Y (r)Y (r + 1) H Heis.Zig−Zag = r S(r) · S(r + 1) + J 2 S(r) · S(r + 2) where X, Y , Z are Pauli matrices, S = [X, Y, Z], and whereZ,X are three-state Potts spin matrices given byZ The next-nearest neighbor coupling in the Heisenberg zig-zag chain is set at the critical value J 2 = 0.24116; at this value the model is scale-invariant [46]. Note that, although the standard Heisenberg model (with J 2 = 0) is quantum critical, the Hamiltonian contains a marginally irrelevant contribution that breaks scale invariance, which here we remove by adding the next-nearest neighbor coupling.
In the present calculation, we have used the modified binary MERA scheme, depicted in Fig. 5(iii), with either M = 2 or 3 transition layers. This ansatz is optimized by minimizing the expectation value of the energy density for each of the above Hamiltonians, by using the optimization algorithm described in Sect. A. In terms of the bond dimension χ MERA , the cost of optimizing the modified binary MERA scales as O(χ 7 MERA ), but can be reduced to O(χ 6 MERA ) through use of an approximation in the tensor network contractions as described Sect. C. We have computed the ground states of the four models over a range of values of χ MERA up to χ MERA = 150. Each simulation took under a week on a 3GHz dual-core workstation with 32Gb of RAM.
For purposes of comparison, we have also computed the ground state of the four critical spin chains using an infinite, translation invariant (with a 4-site unit cell) MPS. The MPS tensors are optimized with a variational approach similar to the iDMRG algorithm [47]. The computational cost scales with the bond dimension χ MPS of the MPS as O(χ 3 MPS ). We have computed the ground states of the four models over a range of values of χ MPS up to χ MPS = 1536.
In both the MERA and MPS calculations we have employed symmetry preserving tensors, as described in Sect. III B, to enforce (some of) the global internal symmetries of these models. Specifically, Z 2 symmetric tensors have been used for the Ising model; Z 3 symmetric tensors have been used for the Potts model (Z 3 is a subgroup of the full S 3 symmetry of this model); and U (1) symmetric tensors have been used for both the quantum XX and Heisenberg zig-zag chains (again, U (1) is a subgroup of the full SU (2) symmetry of the Heisenberg zig-zag chain).
In the first part of the benchmark, Sect. V A, we compare ground energy and twopoint correlators obtained from MERA and MPS, and discuss the relative merits of each approach. Then in Sect. V B we demonstrate the extraction of conformal data from the scale-invariant MERA for the critical Ising model, following the approaches of Refs. [16,22].

A. Comparison with MPS
Here we compare the performances of MPS and scale-invariant MERA for the computation of ground state energy and two-point correlators.

Ground Energy
For both Ising and quantum XX models, the exact ground energy per site is E = −4/π, while for the three-state Potts and Heisenberg zig-zag chains we use an MPS with χ MPS = 1536 to estimate the ground energy per site at E Potts = −2.4359911239(1) and E Heis.Z.Z. ≈ −1.607784749(1).
Then a first observation is that for both MERA and MPS, for a given bond dimension χ the larger the central charge c the larger the error in the energy is. A second similarity is that for both MERA and MPS the error ∆E in the energy scales polynomially in χ, i.e. to a good approximation, A linear fit in Fig. 11 produced the estimates for the coefficients a and b displayed in Table II. In the large χ regime, the error ∆E is dominated by the coefficients b. Interestingly, the ratio b MERA /b MPS for the four models produces very similar results, namely 1.80, 1.74, 1.72, and 1.56 for the Ising, Potts, XX and Heisenberg zig-zag models respectively. Given that the average ratio is b MERA /b MPS ≈ 1.72, we conclude that in the large χ limit a similar error in the energy for MERA and MPS is obtained if that is, if χ MPS = (χ MERA ) 1.72 . Taking into account that the number of variational parameters in the MERA and MPS scale as (χ MERA ) 4 and (χ MPS ) 2 , this comparison shows that in the large χ limit the MPS requires less variational parameters than the scale-invariant MERA in order to obtain a similar accuracy in the ground state energy. It is tempting to extend this comparison to computational costs. A first step in this direction is to note that each iteration in the optimization of MERA and MPS scales (naively) as (χ MERA ) 6 and (χ MPS ) 3 , from which it would be tempting to conclude that MPS algorithms require a lower computational budget than MERA algorithms to obtain the same accuracy in the ground state energy. However, there are important multiplicative prefactors k MERA and k MPS modifying the naive scaling of costs. In both cases, one is required to find the dominant eigenvector of a transfer matrix. But while in the case of the MERA this transfer matrix (or scaling superoperator) has a well-defined gap, implying that k MERA is essentially independent of χ MERA , in the case of the MPS the gap in the transfer matrix closes to zero with increasing χ MPS , and the prefactor k MPS actually also grows with growing χ MPS . Therefore a proper comparison of computational costs requires first a careful characterization of the dependence of k MPS in χ MPS , which is beyond the scope of the present manuscript.

Two-Point Correlators
Let us now compare the accuracy of two-point correlators produced by the scaleinvariant MERA and MPS. For both approaches, we consider the ground state of the quantum XX model and compute the correlator, whereâ is a fermionic operator defined in terms of spin operators aŝ The correlation function of Eq. 44 has the exact expression which is obtained from mapping the quantum XX model to a free-fermion model [48].
The correlation function C(d) decays polynomially, as expected from the ground state of a quantum critical model.
The figure illustrates quite clearly that while the MPS can only approximate polynomial correlations up to a finite length scale ζ = χ κ [33], with κ ≈ 1.38 for the quantum XX model, correlations in the MERA decay polynomially at all length scales [9,16]. Fig.  12(ii) displays the relative error in correlators, Here it is seen that a MERA and an MPS that produce the same accuracy in the ground state energy produce also similarly accurate correlators at short distances, but the relative error in the correlator grows much slower with distance in the case of the MERA. For instance, the χ = 96 MERA reproduces correlators up to d = 10 6 sites with relative error ∆C < 2 × 10 −4 , whereas the χ = 512 MPS, although possessing a similar ground energy, reproduces similarly accurate correlators only up to d ≈ 100 sites.
To summarize, we have seen that the MPS is more efficient than the scale-invariant MERA when it comes to ground state energies of critical Hamiltonians, in that it requires less variational parameters to achieve the same accuracy. However, the MERA produces better correlators at large distances, and it is therefore better suited to characterize asymptotic behaviors, such as the polynomial decay of correlations, from which one could in principle extract the critical exponents of the theory. However, as exemplified in Sect. IV D with the computation of scaling dimensions for local and non-local operators, critical exponents and other conformal data can actually be extracted more directly by analyzing the scaling superoperator. We illustrate this next.

B. Evaluation of Conformal Data: The Ising Model
As an example of extraction of conformal data from the scale-invariant MERA, here we identify the whole operator content (local and non-local primary fields) of the CFT corresponding to the quantum critical Ising model, reproducing the analysis of Ref. 22. Similar results have also been obtained for the three-state Potts and quantum XX models in Refs. [16,22].
The Hamiltonian H Ising of Eq. 36 has a global internal Z 2 corresponding to flipping all the spins. That is, G = Z 2 and g ∈ {+1, −1}, with V +1 = I and V −1 = Z, and The tensors u and w that comprise the scale-invariant MERA are chosen to be parity preserving; each index i of tensors u and w decomposes as i = (p, α p ), where p labels the parity (p = +1 for even parity and p = −1 for odd parity) and α p labels the distinct values of i with parity p. For tensors u, w to be parity preserving it is ensured that, e.g. u j1,j2 i1,i2 = 0 if p(i 1 )p(i 2 )p(j 1 )p(j 2 ) = −1, in accordance with Eq. 16. An operator O acting on the spin chain has parity p if The local scaling superoperator S g=1 and the non-local scaling superoperator S g=−1 , see Fig. 10(iv-v), are obtained from the optimized scale-invariant MERA with bond dimension χ = 32. The scaling superoperator S g=1 is diagonalized to find the local scaling operators φ +1,α together with their scaling dimensions ∆ +1,α , while the scaling superoperator S g=−1 is diagonalized to find the non-local scaling operators of the form together with their scaling dimensions ∆ −1,α . Table III compares the exact scaling dimensions of the primary fields of the Ising CFT and their numerical estimates obtained from the MERA, which reproduce the former with 4 to 6 digits of accuracy in all cases. In Fig. 13 we plot the scaling dimensions of magnitude ∆ ≤ 2.5 obtained from the scale-invariant MERA, which correspond to the primary fields and their descendants, organized both according to the locality g and the parity p of the corresponding scaling operators. Local scaling operators (g = +1) with even parity (p = +1) form the two conformal towers of the primary fields identity I and energy of the Ising CFT, whereas those with odd parity (p = −1) form the conformal tower of the primary field spin σ. Non-local scaling operators (g = −1) with even parity (p = +1) form the conformal tower of the disorder operator µ, and those with odd parity (p = −1) are organized according to two towers corresponding to the fermion operators ψ andψ. The numerical results from the scale-invariant MERA are seen to accurately reproduce the smallest scaling dimensions, those with ∆ ≤ 2.5, from all six conformal towers of the Ising CFT [36,37]. We have also computed the OPE coefficients C αβγ for all primary fields, obtained by analyzing three-point correlators as described in Refs. [16]. Table IV shows the numerical estimate of all non-vanishing OPE coefficients. Once again the results match their exact values to within several digits of accuracy. Thus, not only have we been able to identify the entire field content {I, , σ, ψ,ψ, µ} of the Ising CFT from a simple and rather inexpensive analysis of a quantum spin chain, but we can now also identify  0.068% all possible subsets of primary fields that close a subalgebra by inspecting Table IV. Indeed, it follows that we have the fusion rules (as well as other, such as σ × µ = ψ +ψ, etc) from where we see that {I, } and {I, , σ} close subalgebras of local primary fields, whereas {I, , µ} and {I, , ψ,ψ} close subalgebras that contain both local and non-local primary fields, where locality is relative to the spin variables.

VI. CONCLUSIONS
In this manuscript we have presented an introduction to the scale-invariant MERA and its application to the study of quantum critical systems. The main strength of MERA, when applied to quantum critical systems, is that it can explicitly incorporate scale invariance. This facilitates enormously the computation of scaling dimensions (equivalently, critical exponents) and of other properties that characterize a quantum phase transition.
Direct comparison with an MPS shows that, while the later is more efficient at computing local observables such as the ground state energy, the MERA produces significantly more accurate correlators at long distances. In addition, from the MERA it is straightforward to identify the scaling operators of the theory, as well as their scaling dimensions and operator product expansion, producing accurate conformal data that can be used to unambiguously identify the underlying CFT.
Here we have considered homogeneous systems. However, the scale-invariant MERA has been successfully generalized to critical systems where translation invariance is explicitly broken, as it is the case of a critical system with a boundary, with an impurity, or the interface between two critical systems [49][50][51]. In all these scenarios translation invariance is no longer present, but exploitation of scale invariance still produces a MERA algorithm with a cost O(1) (that is, independent of the lattice size N ), and therefore infinite systems can be addressed. Again, scaling operators associated to boundaries, defects and interfaces can be easily extracted from the simulations.
Finally, much of the MERA formalism for critical systems in D = 1 dimensions is also directly applicable to D = 2 dimensions, including the characterization of scaling dimensions. However, due to significantly larger computational costs, so far only small values of χ have been used in actual computations. Thus, further progress needs to be made in reducing computational costs before the scale-invariant MERA becomes a viable approach to quantum criticality also in D = 2 dimensions.

Appendix A: Optimization algorithm for the scale-invariant MERA
In this section we describe an algorithm to optimize the scale-invariant MERA to approximate the ground state of a critical system. It is based on modifying the optimization algorithm for a regular MERA of Ref. [15] and has been previously sketched in Refs. [15,16], although the present implementation differs in some details.
A scale-invariant MERA is composed of a small number M of transitional layers U [0] , U [1] , . . . , U [M −1] followed by an infinite sequence of (identical) scale-invariant layers, here denoted with an asterisk, U * . We shall henceforth use similar asterisk notation for all operators and tensors associated to the scale-invariant layers. Each of the transitional layers U [τ ] are characterized by a single isometry w [τ ] and disentangler u [τ ] ; similarly, the scale-invariant layers U * are characterized by a single isometry w * and disentangler u * . The goal of this section is to describe how these tensors can be optimized in such a way that, given a critical Hamiltonian, the scale-invariant best approximates its ground state.
We begin, in Sect. A 1, by describing the key building blocks of the optimization algorithm, which are (i) coarse-graining of the Hamiltonian; (ii) fine-graining of the density matrix; and (iii) optimization of one tensor of the MERA. Then, in Sect. A 2, we describe how these algorithmic building blocks can be put together to form an iterative optimization scheme. Finally, Sect. A 3 describes computational tricks to improve convergence and accuracy.

Building blocks a. Coarse-graining of the Hamiltonian
A key part of the MERA optimization algorithm is the coarse-graining of the Hamiltonian. Given an initial Hamiltonian H [0] = r h [0] (r, r + 1) that decomposes as a sum of identical local terms h [0] (r, r + 1) = h 0 , we wish to construct, for any level τ , the coarse-grained Hamiltonian H [τ ] = r h [τ ] (r, r + 1), defined as The coarse-graining of an operator that decomposes as a sum of local operators, such as a local Hamiltonian, can be achieved with the ascending superoperator formalism introduced in Sect. II A for the coarse-graining of a single local operator. The coarse-grained Hamiltonian coupling h [τ +1] is obtained by enacting the (left,center,right) ascending superoperators A L , A C and A R on the Hamiltonian term h [τ ] , whereĀ is the average of the three ascending superoperators. The diagrammatic representation of the tensor network described by Eq. A2 is shown in Fig. 15(i).

b. Fine-graining of the density matrix
Another key part of the MERA optimization algorithm is fine-graining the density matrix. Given the average two-site density matrixρ [τ ] defined on lattice L [τ ] , we wish to construct the average two-site density matrixρ [τ −1] on lattice L [τ −1] . This is accomplished through the average descending superoperatorD introduced in Sect. IV C, which acts asρ From iteration of Eq. A3 one can construct the average two-site density matricesρ [τ ] for all τ < τ . The diagrammatic representation of the tensor network described by Eq. A3 is shown in Fig. 15(ii).
c. Optimization of one tensor of the MERA In order to approximate the ground state of a Hamiltonian H, the disentanglers u and isometries w that define a scale-invariant MERA |Ψ should be chosen to minimize the energy E ≡ Ψ|H|Ψ . We shall proceed by updating one tensor of the MERA at a time, while holding the rest of the tensors fixed. Given a tensor to be updated, we note that the energy E depends quadratically on that tensor and its conjugate. Here we employ the linearized optimization scheme of Ref. [15], which we sketch next (a justification for linearizing the cost function and further details can be found in Ref. [15]).
The update of an isometry w relies on computing its environment Υ w , which represents a factorization of E = Ψ|H|Ψ to isolate its dependence on w, with k 1 an irrelevant constant. [Υ w is the derivative of Ψ|H|Ψ with respect to w while keeping w † fixed]. The updated isometry w that minimizes the linearized energy -that is, such that tr( where V 1 and V 2 are obtained from the singular value decomposition (SVD) of the environment, We then proceed by replacing that particular isometry w with w in the MERA, which now represents a new state |Ψ . We emphasize that the energy E = Ψ |H|Ψ is computed by replacing both w and w † with w and w † , and it is therefore not given by tr(w Υ w ) + k 1 . In other words, w does not minimize the expectation value of the Hamiltonian (in fact, the energy could even rise!). However, the approach works well in practice. The update of a disentangler u follows an analogous procedure, involving the SVD of the corresponding environment Υ u .
The environments of an isometry w [τ ] and disentangler u [τ ] at layer τ depend only on a small number of other tensors, as depicted in Figs. 15(iii)-(iv).

Optimization algorithm
For concreteness, we describe the optimization of a scale-invariant MERA with M = 2 transitional layers as depicted Fig. 14(ii). The ansatz is then completely characterized by three isometries w [0] , w [1] , w * and three disentanglers u [0] , u [1] , u * , where the isometry w * and the disentangler u * define the scale-invariant layers. For purposes of the optimization algorithm, one should also store in computer memory the two-site Hamiltonian couplings h [0] , h [1] , h [2] , h ♦ . Here h ♦ is a scale-averaged Hamiltonian, to be introduced later in Eq. A11. We also store the average two-site density matrices ρ [1] ,ρ * , where we recall thatρ * is the two-site density matrix any of the scaleinvariant layers. We now proceed to describe the optimization algorithm for scaleinvariant MERA, as outlined Fig. 14(iii). In an iteration of the algorithm, comprising Sect. A 2 b, Sect. A 2 c and Sect. A 2 d, each of the stored tensors is updated to a new version of itself denoted with a prime, e.g.

a. Initialization
The preliminary step of the algorithm is to initialize the tensors which define the MERA. It is most often sufficient to initialize the tensors to be random, though in certain situations, for instance if one has prior knowledge about the ground state of the Hamiltonian, it is useful to use that information. Random isometries w [0] , w [1] , w * can be obtained through singular value decomposition of an appropriately sized random rectangular matrix. The disentanglers u [0] , u [1] , u * can also be initialized randomly, or as the identity operator. Notice that the local Hamiltonian term h [0] is given as the input of the algorithm (it describes the critical Hamiltonian whose ground state is to be found), whereas the Hamiltonian terms h [1] , h [2] , h ♦ and density matrices ρ [1] ,ρ * will be generated during the first iteration of the algorithm.
b. Updating the density matrices The first step of each iteration is to update the two-site reduced density matrices. As explained Sect. IV C, the scaling density matrixρ * is defined as the dominant eigenoperator of the average descending superoperatorD * associated to the scaleinvariant layers,ρ * =D * (ρ * ) . (A7) One can solve Eq. A7 forρ * with a sparse diagonalization technique such as the Lanczos method, where the scaling density matrixρ * from the previous iteration (if any) can be reused as the starting point for the Lanczos method to accelerate convergence. The updated density matrixρ [1] is then obtained by fine-grainingρ * with the descending superoperatorD [1] ,ρ and obtain the updated Hamiltonian coupling h [1] , as described Sect. A 1 a. The second transitional layer U [1] is then updated, in the same manner as described here for the first transitional layer, to obtain updated tensors w [1] , u [1] , h [2] .

d. Updating the scale-invariant layers
The final step of the iteration is to update the isometry w * and disentangler u * associated to the scale-invariant layers of the MERA. Updating the scale-invariant layers is more complicated than updating the transitional layers since the scaling tensors w * and u * do not appear in just one layer but in all layers U [τ ] for τ ≥ M , with M the number of transitional layers). To update these tensors, we construct linearized environments that average the contributions coming from all length scales τ ≥ M . Computation of these scale-averaged environments is simplified [15] by the introduction of the scaleaveraged Hamiltonian term h ♦ , which for the case of M transitional layers is defined as where the factor of (1/3) τ −M arises from the fact that layer U [τ ] contains three times as many tensors as layer U [τ +1] [15]. Note that for a critical Hamiltonian (with negligible or sufficiently suppressed RGirrelevant terms) and a scale invariant MERA that has already been optimized, scaleinvariance implicitly assumes that h [τ +1] ∝ h [τ ] for all τ ≥ M , and therefore the scale-averaged Hamiltonian h ♦ is simply proportional to h [M ] . However this property will only hold for a MERA which has been properly optimized; during the optimization procedure it is necessary to explicitly average over all scale-invariant layers when constructing environments. Computing the scale-averaged Hamiltonian h ♦ directly through Eq. A11 is not feasible due to the infinite summation. One strategy to obtain an approximation to h ♦ is to use a partial summation of Eq. A11 which stops at some τ = T . In practice, however, a useful estimate of the updated scale-averaged Hamiltonian h ♦ is already obtained from the operator h ♦ from the previous iteration through from which, as usual, one obtains the updated tensors w * and u * by SVD. This concludes a single iteration of the optimization algorithm for scale-invariant MERA; the next iteration can begin again at Sect. A 2 b.
where I is the identity operator on two sites, such that the shifted Hamiltonian H α is negative defined. This can be achieved by choosing α in Eq. A15 as the largest eigenvalue of h [0] . Having a negative defined Hamiltonian ensures that the optimization targets the low-energy, i.e. ground state, subspace (the optimization algorithm, based upon extremizing the energy of the state, could otherwise target the high-energy subspace).
It is also useful to similarly shift the spectrum of the coarse-grained Hamiltonian that arise during the optimization, by replacing α where α as the largest eigenvalue of h [τ ] , since this is seen to speed-up convergence during the optimization of the MERA.

c. Converging the number of transitional layers
In order to obtain an accurate approximation to the ground state of critical Hamiltonian H with a scale invariant MERA it is important to use enough transitional layers to sufficiently suppress the effect of any RG irrelevant terms present in H, which can break scale-invariance at short length scales. For a given critical H the appropriate number transitional layers is not known a priori, and must be found by trial and error. We proceed in the following way. Suppose that we have already optimized a MERA with M = 2 transitional layers, that is, characterized by U [0] , U [1] , U * . We can then add an additional transitional layer, M = 3, so that now the MERA is given by U [0] , U [1] , U [2] , U * , where we initially set U [2] = U * , and re-optimize the tensors as to best approximate the ground state of Hamiltonian H. We then keep adding more transitional layers, and re-optimizing the ansatz, until addition of an extra transitional layer does not produce significantly different results under re-optimization in terms e.g. of the ground energy, scaling dimensions, etc.
In Sect. II C three different implementations of the MERA in D = 1 dimensions were described: the binary MERA, the ternary MERA and the modified binary MERA. The three schemes differ in how the computational cost (incurred e.g. in optimizing the ansatz and in computing the expectation value of local observables) scales with the bond dimension χ, namely as O(χ 9 ), O(χ 8 ) and O(χ 7 ), respectively. In addition, the three schemes also differ in the strength of disentangling, which means that for the same value of χ the numerical accuracy of the results will vary across the three schemes. In this section we investigate which of the schemes provides the most accurate ground state energy for a fixed computation cost.
For this purpose, we have optimized each of the three MERA schemes so as to approximate the ground states of the critical Hamiltonians defined Eqs. 36-39 in Sect. V. Fig. 17 shows the relative energy error ∆E, defined as as a function of bond dimension χ. For all three schemes the relative energy error ∆E is seen to scale as a power of the bond dimension χ, for some coefficients 'a' and 'b' that depend both on the MERA scheme in use and the spin model under consideration. The sets of these coefficients, obtained from linear fits in Fig. 17, are displayed in Table V. For each of the four spin models separately, the data shows that while the coefficient 'a' depends considerably on the MERA scheme (with variations by one or two orders of magnitude), the coefficient 'b' is remarkably constant across the schemes (with variations of about 15%). These results indicate that, while using a scheme with more powerful disentangling (such as the binary MERA scheme) can reduce the energy error by a considerable, fixed multiplicative factor (as expressed by a significantly smaller constant 'a'), the power 'b' of χ for a given critical model is not significantly affected by the disentangling power of the scheme, even though the latter does affect significantly how the computational cost scales with χ. This suggests that the modified binary MERA scheme, which has the smallest disentangling power, will be the most cost effective at large χ, since its cost scales as a smaller power of χ and the accuracy scales roughly as the same power of χ as the other schemes. Fig. 18 shows ∆E directly as a function of the computational cost for the Ising model. Similar results are obtained for the other models. The costs have been approximated to be C = χ 9 for binary, C = χ 8 for ternary and C = χ 7 for modified binary MERA, where C is measured in floating-point operations. The exact expression for the cost would include a multiplicative constant of order one that does not change significantly across the schemes and has been ignored for simplicity. The plot confirms that, for large χ, the modified binary MERA scheme, with cost C = χ 7 , is the scheme that provides the most accurate results for a fixed computational cost, whereas for smaller values of χ the binary scheme is the most cost effective. The crossover occurs at a cost of order C ≈ 10 7 floating-point operations per iteration, which is well within the capabilities of a small workstation. Therefore, in most practical numerical applications the modified binary MERA is a better scheme than the others.  In Appendix B, a numerical study determined the modified binary MERA scheme, which can be optimized with a leading computational cost which scales as O(χ 7 ), to be the optimal 1D MERA implementation of those considered. In this Appendix we describe how the cost of implementing the modified binary MERA scheme can be reduced from O(χ 7 ) to O(χ 6 ) through the use of approximations in the tensor network contractions required for its optimization. The steps required to optimize the modified binary MERA are analogous to those described for the ternary MERA in Sect. A. Fig. 19 shows four closed tensor networks (tensor networks without any open index). The tensor networks required to optimize the modified binary MERA (which include the ascending and descending superoperators, and the environments for single tensors) can be generated from these closed tensor networks through removal of a single tensor. The cost of contracting the two networks in Figs. 19(i-ii) scales as O(χ 7 ), whereas cost of contracting the two networks in Figs. 19(iii-iv) scales just as O(χ 6 ) contraction cost. If we could reduce the cost of the two first networks also down to O(χ 6 ), then that would be the overall leading cost of the modified binary MERA. Let us then see how this can be accomplished.

Insertion of Projectors
Let us consider a rankχ projector P , decomposed as the product of an isometric tensor v and its conjugate v † ,  where V χ is a χ-dimensional vector space and Vχ is aχ-dimensional vector space for someχ ≤ χ 2 . In place of using the original tensor networks of Fig. 19(i-ii) in the optimization algorithm, with contraction cost O(χ 7 ), we shall use the tensor networks of Fig. 19(v-vi), which have been modified through inclusion of the projector P and whose cost is now of order O(χ 5χ ). In the reminder of this Appendix we first explain how to optimally choose the isometric v so that the modified tensors networks yield environments that best approximate the environments obtained from the original networks; and then we argue, based on numerical evidence, thatχ may indeed be chosen as O(χ) without significant loss of accuracy, thus reducing the overall cost of optimization down to O(χ 6 ).
FIG. 20. The spectrums of one-site and two-site density matrices, obtained from the scaleinvariant fixed point of a χ = 16 MERA optimized for the ground state of the quantum XX model. The two site density matrix can be truncated down to itsχ = 50 most significant eigenvalues, out of the total χ 2 = 256 eigenvalues, without significant loss of accuracy (as gauged by the smallest eigenvalue of the one-site density matrix).

Optimizing the Projector
Given the two-body density matrix ρ [τ ] at level L [τ ] of the MERA, the proper choice for the isometric v [τ ] is that which maximizes the trace of the density matrix,  Fig. 19(v-vi) will give the same result as the original tensor networks of Fig. 19(i-ii). Likewise the other tensor environments (such as environments of disentanglers and isometries) generated from the modified tensor networks will be the same as those generated from the original networks.
In principle the isometric v [τ ] that projects onto the most significant subspace of ρ [τ ] , as described Eq. C2, could be obtained directly from the spectral decomposition of the density matrix ρ [τ ] for each τ ; however, computing the exact reduced density matrix ρ [τ ] from ρ [τ +1] is an O(χ 7 ) operation and should be avoided. By expressing the density matrix ρ [τ ] in terms of ρ [τ +1] , u [τ ] and w [τ ] , we can write the trace in Eq. C2 as the tensor network in the left hand side of Fig. 19(vii). From that tensor network, we can extract the environment for v with cost O(χ 5χ ). Therefore, by using the linearized single tensor optimization described in Sect. A 1 c we can iteratively optimize v at a cost O(χ 5χ ).

Rank of Projector
The introduction of a properly optimized rankχ projector P was argued to be equivalent to truncation of the two-site density matrix to retain itsχ most significant eigenvalues; this leads to truncation error ε = 1 − tr P [τ ] ρ [τ ] that should be kept sufficiently small. A good indication of exactly how small ε should be kept is to ensure that it is the same size, or smaller than, the smallest eigenvalue of the one-site density matrix, which is indicative of the degree of accuracy of the MERA under consideration. Thus the relevant question becomes, for a MERA with bond dimension χ, what is the necessary rankχ of the projector P required to maintain this sufficient degree of accuracy?
Since the MERA represents states with at most a logarithmic scaling of the entropy with block size, see Sect. IV A, the entropy of the two-site density matrix is much less than twice the entropy of the one-site density matrix. Therefore, it is to be expected thatχ χ 2 . [Notice thatχ = χ 2 would be consistent with an entropy that scales linearly with the block size]. Indeed, the numerical evidence suggests thatχ may be chosen as O(χ); for the simulations with the modified binary MERA scheme of Sect. V, it was found that keepingχ ≈ 3χ gave a sufficient level of accuracy for all four critical models under consideration, and over the large range of bond dimensions χ analyzed. Given that this relation held over a range of bond dimensions between χ = 4 and χ = 150, it seems likely that the relation between χ andχ is indeed linear, or at least very close to linear for ground states of critical systems. An example is shown in Fig. 20, which plots the spectrum of one-site and two-site density matrices for χ = 16 quantum XX model, where it can seen that choosingχ ≈ 3χ is sufficient to ensure that the truncation error ε is of the same order as the smallest eigenvalue of the one-site density matrix.
Thus, though not rigorously justified, the available evidence indicates thatχ may be chosen as O(χ), hence the overall cost of optimizing the modified binary MERA has been reduced to O(χ 6 ) cost. The reduction in cost comes at the price of both introducing new tensors v [τ ] that must be updated with each iteration and introducing a controlled amount of approximation into the tensor network contractions. Although we have only described how the cost of the modified binary MERA scheme can be reduced through introduction of a projector P into the tensor network diagrams, the same approach could be employed to potentially reduce the cost of any MERA scheme.