Entanglement Renormalization for Quantum Field Theories with Discrete Wavelet Transforms

We propose an adaptation of Entanglement Renormalization for quantum field theories that, through the use of discrete wavelet transforms, strongly parallels the tensor network architecture of the \emph{Multiscale Entanglement Renormalization Ansatz} (a.k.a. MERA). Our approach, called wMERA, has several advantages of over previous attempts to adapt MERA to continuum systems. In particular, (i) wMERA is formulated directly in position space, hence preserving the quasi-locality and sparsity of entanglers; and (ii) it enables a built-in RG flow in the implementation of real-time evolution and in computations of correlation functions, which is key for efficient numerical implementations. As examples, we describe in detail two concrete implementations of our wMERA algorithm for free scalar and fermionic theories in (1+1) spacetime dimensions. Possible avenues for constructing wMERAs for interacting field theories are also discussed.

1 Background: Entanglement Renormalization Ongoing advances in quantum technologies have opened up the possibility that quantum field theories (QFTs) might become simulatable in quantum computers in the medium-tolong term future, offering new methods to address long-standing and intractable problems that we have no hope to solve with classical computers.With such prospects, much thought and research efforts have been recently directed towards developing tools that are deliberately bypassed in the traditional, Euclidean space treatments of QFTs.Such tools include the construction of explicit Hilbert space representations of ground and excited states of QFTs-including their entanglement patterns-and real-time Hamiltonian evolution.
In building and understanding these tools, it pays off to borrow some of the wisdom developed in condensed matter physics, which has traditionally used the framework of wavefunctions, entanglement, and coarse graining to understand the properties and dynamics of condensed matter systems.In particular, a significant insight that emerged in the early 90s was that entanglement plays a significantly role in how fundamental degrees of freedom reorganize themselves as effective degrees of freedom under Renormalization Group (RG) evolution.
Until the early 90s, there were notorious difficulties in numerically implementing realspace block-spin renormalization of quantum spin chains.In block-spin renormalization, a coarse-graining step consists of mapping a block of sites of the fine-grained lattice into a single site of the coarse-grained lattice via a decimation procedure that reduces the dimensionality of the Hilbert space.In Wilson's numerical RG [1], this is implemented via diagonalization of the Hamiltonian restricted to the block's Hilbert space and its subsequent truncation, keeping only the lowest energy eigenstates of the block.This procedure, however, generically fails to retain the whole support of the ground state, and provides rather unreliable results for a variety of lattice systems [2][3][4], with the exception of Hamiltonians which can be put into a special form.In 1992, White and Noack [5] understood the fundamental difficulty with this approach: by neglecting interactions of the block with the rest of the lattice, Wilson's numerical RG fails to capture entanglement between neighboring blocks, and therefore cannot preserve ground state properties, regardless of any reasonable increase in the number of states kept.
The shortcomings of the standard numerical RG approach were overcome by White's Density Matrix Renormalization Group (DMRG) [6].Instead of the block's Hamiltonian spectrum, this approach considered the spectrum of the block's reduced density matrix, and truncated its Hilbert space by discarding eigenstates of the reduced density matrix with eigenvalues smaller than a given threshold, determined by the desired level of accuracy.In this way, the entanglement properties of the ground state could be retained to any desired precision and a proper RG flow could be generated.
Eventually, the DMRG was understood in terms of the quantum information theoretic concept of tensor networks (TNs), and recast as a variational optimization method over the class of TNs known as Matrix Product States (MPS) [7,8].Essentially, TNs represent quantum many-body states in terms of networks of interconnected tensors, which in turn capture the relevant entanglement properties of the wave function.TNs provide an unbiased approach to numerical quantum many-body physics; they have been shown to be highly efficient and accurate tools to compute physical observables, while also being free from the sign problem that plagues Monte Carlo sampling methods in Euclidean space (such as in, e.g., Lattice QFT), and also suitable to simulate real-time dynamics.
The success of TNs in describing properties of quantum many-body systems can be attributed to the efficiency of optimized TNs in identifying and keeping track of the relevant degrees of freedom.This, in turn, is directly related to the fact that the low energy states of realistic, gapped local Hamiltonians are not generic states in the Hilbert space: they are heavily constrained by locality, and their entanglement is limited by area laws.In particular, a key insight of TNs is that, for systems described by local Hamiltonians, coarse graining must be preceded by the removal of short-range entanglement by local unitary transformations, so that short-distance degrees of freedom can be safely traced out while preserving ground state entanglement properties [5,9].A particular class of TNs in which this is self evident and intuitive is called MERA, an acronym for Multiscale Entanglement Renormalization Ansatz [10].
In MERA, short range entanglement across the boundaries of adjacent blocks is first removed by local unitary transformations, or disentanglers, and then a coarse-graining transformation can be applied to trace out short distance degrees of freedom while safely preserving ground state entanglement properties.An important consequence of the use of disentanglers is that the dimension of the effective sites no longer needs to grow with each Figure 1.A (non-technical) depiction of the MERA architecture.From bottom to top, it can be viewed as an RG flow, whereas from top to bottom, it can be viewed as the construction of the ground state ansatz described in the text.
coarse graining step, and therefore an efficient description of critical fixed points becomes feasible [11,12].MERA has been successfully applied to study strongly interacting systems in low dimensions, including the emergence of symmetry breaking order and topological order [13,14], quantum critical points [12,15], real time evolution, fermionic and anyonic systems (which are intractable by quantum Monte Carlo) [16,17], as well as a lattice realization of the holographic principle [18,19].
The construction of a MERA to represent the ground state wavefunction of a system essentially amounts to a reverse RG flow.Consider, as a matter of example, a lattice system governed by a local Hamiltonian.The starting point of the MERA is an unentangled state of n lattice sites, e.g., the product state |P ⟩ = |0⟩ ⊗n .The reverse RG flow starts by acting on |P ⟩ with quasi-local unitary transformations 1 to create entanglement between neighboring sites.These quasi-local unitaries are called entanglers, and their collective action over the entire system is represented by an operator U 0 such that (1.1) A "fine-graining" mapping then follows by the introduction of n new lattice sites intercalated with the original sites, effectively reducing the lattice spacing by half.These new sites are also initialized as a product state, so that the augmented, finer-lattice state is given by The next step is to entangle the newly introduced sites with the rest of the lattice via new quasi-local entanglers represented by U 1 , yielding New layers of fine-graining and entangling operations can then be successively applied to flow deeper into the UV, until a ground state ansatz at the desired UV resolution (or until a fixed point) is reached, namely, The parameters of the entanglers U j are variationally optimized (by minimizing the expectation value of the Hamiltonian) to recover the ground state of the lattice system.A (non-technical) depiction of the MERA described above is shown in figure 1.
In the example described above, the entanglement introduced between adjacent sites at scale j gives rise to entanglement between 2 (j UV −j) sites in the final ansatz |Ψ UV ⟩.Hence, this iterative approach of introducing quasi-local entanglement at each scale allows one to capture long distance correlations with a series of manageable, quasi-local unitary operations.In essence, MERA works because (i) entanglement is quasi-local in position space (this stems from interactions being local) and (ii) the entanglement between effective degrees of freedom at a given resolution is insensitive to shorter-range entanglement at much higher resolutions (this stems from the decoupling principle [20], i.e., the insensitivity of a low energy effective theory to its UV completion).
The general properties aforementioned are not exclusive to condensed matter systems described by lattices; indeed, locality of interactions and decoupling are also properties of continuum relativistic systems such as local QFTs.Hence, one would expect that the MERA principles could be generalized to build variational classes for QFTs.And indeed, in 2011, Haegeman et al. [21] proposed such a generalization, which they coined "continuum MERA," or cMERA for short.The elements of the original, discrete MERA were conceptually modified in their adaptation to the continuum.In particular, Haegeman et al. did not propose a discretization of the field theory to more closely parallel the lattice structure and countable layers of the original MERA; instead, cMERA was formulated directly in the continuum, with the system's resolution determined by a Wilsonian UV cut-off Λ UV .In lieu of a discrete layered architecture, continuous scaling transformations were used to flow the system from one scale to another.In particular, the role of the fine-graining mapping was formally implemented by infinitesimal dilation operations to increase the UV cut-off Λ UV under RG flow.Finally, the action of the unitary entangler operator was implemented in a basis of momentum-space degrees of freedom, effectively generating correlations at all scales up to the shortest distance cut-off, 1/Λ UV .While formally attractive and initially promising, the cMERA formulation turned out to be unwieldy and did not lead to any notable results, except in examples of free field theories [22,23].Specifically, the path-ordered dilation transformations required to implement changes in scale proved to be impractical for any non-trivial (i.e., interacting) field theory.More problematic, however, was the fact that the implementation of entangling operations in momentum space obscured the locality of entanglement and failed to explicitly capture a key feature of MERA, namely, the hierarchical pattern of entanglement between localized, neighboring degrees of freedom at each scale.
The shortcomings of the cMERA formulation mentioned above were never addressed in the many attempts to implement and/or generalize it to interacting field theories [24][25][26][27][28].
This hindered any meaningful progress in the study of QFTs with cMERA since it was first proposed over a decade ago.In particular, no useful/practical entangler ansatzes for interacting QFTs have been fully demonstrated to date.
In this article, we re-examine the generalization of the MERA principles for field theories, and propose an alternative to cMERA which preserves, in explicit form, key features of Entanglement Renormalization, namely, (i) locality of entanglement and (ii) the layered architecture of the original MERA.We coin this alternative MERA adaptation to QFTs "wMERA."A key underlying principle of wMERA is the use of Discrete Wavelet Transforms (DWT) to discretize continuum quantum fields into countable degrees of freedom localized in position space2 .Intuitively, DWTs encode low-and high-pass filters that, upon a choice of scale, separate the degrees of freedom of a quantum field into short-range fluctuations and long-range fluctuations.This simple basis choice provides a natural framework for RG evolution, while keeping locality and separation of scales in explicit form, which makes the MERA generalization quite straightforward.In particular, the fine-graining transformations are implemented via Inverse Wavelet Transforms, which add shorter-range (wavelet) degrees of freedom at each inverse RG step to flow the system to the UV.
The outline of this article is as follows.In section 2, we provide a brief overview of cMERA and highlight its main shortcomings.In section 3, we introduce the concept and basic machinery of discrete wavelet transforms, and the specific class of Daubechies wavelets.In section 4, we provide the generic formulation of wMERA using DWTs, the concrete implementation of which is then illustrated in section 5 for two concrete examples, namely, a free scalar field and a free Dirac fermion in (1+1)d.Our concluding remarks and outlook are presented in section 6.

Brief overview of cMERA
In 2011, Haegeman et al. [21] introduced cMERA as a formal adaptation of MERA to the continuum.Briefly, in cMERA the ground state wavefunction of a system at its infrared (IR) fixed point, |Ω IR ⟩, is related to its microscopic counterpart at an ultraviolet (UV) scale s UV , |Ψ(s UV )⟩, via the following unitary transformation: Above, P denotes path-ordering in scale s, L is the generator of scaling transformations, and K(s) is the generator of entanglement at scale s.Imposing renormalization group (RG) invariance of physical observables, they obtained the RG equation for a generic hermitian operator Ô(s): Haegeman et al. [21] illustrated the cMERA construction for a few examples of free field theories 3 .For instance, for a free scalar field theory in (d + 1) dimensions, the entangling exponent in (2.1), K(s) + L, was given by the following Gaussian generators: and where φ( ⃗ k) and π( ⃗ k) are the momentum-space modes of the scalar field and its canonical conjugate momentum, respectively; m ϕ is the mass of the scalar field; and Λ is an ad hoc UV cut-off.While formally elegant and concise, the cMERA formulation turned out to be cumbersome even for the simplest case of free field theories.As we shall see in section 5, the ad hoc cut-off Λ and the path-ordered architecture of scaling transformations in (2.1) is completely unnecessary for free field theories.More importantly, however, they are impractical for any nontrivial (i.e., interacting) field theory, for which the entangler generator K(s) must be generalized beyond quadratic order in the canonical fields to be able to capture nontrivial n-point correlation functions beyond the two-point function.
Compounded by these obstacles is the fact that all incarnations of cMERA attempted so far in the literature were, to the best of this author's knowledge, implemented in momentum space.Since interactions are not local in momentum space, neither is entanglement.More problematic is the fact that interactions entangle momentum-space degrees of freedom across vastly different scales.Hence, at the conceptual level, cMERA constructions in momentum space obscure the hierarchical pattern of entanglement between localized degrees of freedom in position space at each scale.At the practical level, momentum-space cMERA is difficult to implement numerically, since it does not share the original MERA's layered architecture, nor the quasi-locality of its entanglers.
To address these obstacles, we propose an alternative adaptation of MERA to the continuum that is directly formulated in position space, and that does away with the use of continuous path-ordered scaling transformations.In our proposed formulation, which we coin "wMERA," the layers of the MERA are obtained via a Multiresolution Analysis, whereby the quantum fields are expanded in a wavelet basis.There is several practical advantages of wMERA over cMERA: (i) the wavelet discretization of quantum fields automatically regulates the theory in the UV; (ii) the locality of degrees of freedom preserves the quasi-locality of entangler ansatzes; (iii) separation of scales is explicitly built-in, enabling the layered architecture of MERA to capture the hierarchical pattern of entanglement across different scales; and (iv) with the degrees of freedom and scales automatically discretized, numerical methods developed for the original MERA can be adapted to wMERA.
Before describing wMERA in detail, we first we provide a brief introduction to Discrete Wavelet Transforms and, in particular, to Multiresolution Analysis with Daubechies wavelets, in section 3.

Brief overview of DWTs with Daubechies wavelets
Wavelet transforms offer a flexible framework to deconstruct the features of a function, image, or sequence in a hierarchy of resolutions.They are extensively used in signal processing and data compression.
In quantum field theories, on the other hand, Fourier Transforms (FTs) have been traditionally used to decompose the degrees of freedom of a quantum field into momentum modes.While such decompositions are optimal for problems dealing with scattering of plane waves, the delocalized nature of Fourier modes makes FTs highly inefficient in representing spatially localized features, and, by extension, the pattern of ground state entanglement of local, interacting field theories.A suitable alternative to FTs, which we will use in our wMERA formulation, are discrete wavelet transforms (DWTs), which will allow us to decompose the degrees of freedom of a quantum field into localized modes in a hierarchy of resolutions 4 .
In order to intuitively understand some of the key features of DWTs, we begin by describing Haar wavelets, which are arguably the simplest wavelet bases there are.

Haar wavelets
The Haar wavelet bases are generated from discrete translations and rescalings of the Haar scaling function, which is essentially the unit box function centered around x = 1/2: Suppose one wishes to represent square integrable functions on the real line, L 2 (R), with a resolution of ∆x = 1.An orthonormal basis for such decomposition can be constructed from discrete translations of the Haar scaling function s(x): Above, T is the unit translation operator.
x f(x) Now, suppose one wishes to improve the resolution of the decomposition basis by a factor of 2, i.e., to ∆x = 1/2.In this case, the appropriate Haar basis would consist of rescaled scaling functions which are narrower by a factor of 2 (to improve resolution) and taller by a factor of √ 2 (to preserve unit normalization): Above, D is the dyadic rescaling operator.An example in which a generic function is decomposed in this basis in a sequence of improved resolutions is shown in figure 2.
It is easy to see that the space S 0 spanned by the original basis {s n (x)} is a subset of the space S 1 spanned by the improved-resolution basis {s 1 n (x)}: The orthogonal complement of S 0 in S 1 is a nonempty set called the wavelet subspace W 0 : W 0 is spanned by the Haar wavelet functions w n (x), which are generated from discrete translations of the mother wavelet w(x): In other words, the bases span the same subspace S 1 in L 2 (R).Their elements are related by the following orthogonal transformation: whose inverse is: The transformation in (3.8) defines the Haar wavelet transform (WT), which can be recast in the following form: (3.12) The projection operators L and H in (3.11) and (3.12) are, respectively, the Haar low-pass filter and the Haar high-pass filter.They are defined in terms of the finite array of weight coefficients {h n } and {g n }.With this notation, the inverse wavelet transform (IWT) in (3.9) can be rewritten as: The Haar wavelet transform and its inverse are depicted in figure 3.
If one wishes to represent functions with resolution ∆x = 1/2 r (where the resolution index r is always an integer), a Haar basis can be constructed from r dyadic rescalings of s n (x), yielding basis function elements with width equal to 1/2 r : In particular, the subspace S r+1 spanned by the basis {s r+1 n (x)} n∈Z contains all lower resolution subspaces S r , S r−1 , ..., S 0 in a nested form: In full analogy with (3.5), the orthogonal complement of S r in S r+1 is the wavelet subspace W r : which is spanned by a basis of rescaled wavelets: The generalization of the (inverse) wavelet transform relations between the two equivalent bases {s r+1 n (x)} n∈Z and {s r n (x)} n∈Z ∪ {w r n (x)} n∈Z (3.18) spanning S r+1 is, expectedly, similar to the r = 0 expressions in (3.10) and (3.13): and Combining (3.15) and (3.16), we have a multiresolution decomposition of L 2 (R): i.e., for a base resolution choice r, the set of functions form an orthonormal basis for L 2 (R) with arbitrarily high resolution 5 .The scaling functions s r n (x) encode structure up to scale 1/2 r , and the wavelets w r+j n (x) capture structure at finer resolutions 1/2 r+j+1 .
It is easy to see that the multiresolution decomposition in (3.21) and its associated basis in (3.22) offer the desired architecture upon which to build the layers of a MERA for a QFT ground state.The discretized degrees of freedom of a quantum field ϕ r i , obtained from smearing ϕ(x) over the scaling functions s r i (x/a 0 ), heuristically describe collective field fluctuations over distances of order O(a 0 /2 r ).Shorter-distance fluctuations of the field at higher resolutions r ′ > r are captured by ϕ(x) projections over the wavelet functions w r ′ n (x/a 0 ).Inverse wavelet transforms are a means to recombine scaling-and waveletdegrees of freedom into shorter-distance collective modes, whereas wavelet transforms are a means to tease out UV fluctuations that need to be disentangled to enable coarse-graining.
Unfortunately, the Haar multiresolution decomposition has shortcomings that make it not completely suitable for this intended MERA application.In particular, the Haar scaling and wavelet functions are not differentiable, rendering the kinetic term in the Hamiltonian for discretized field degrees of freedom ill-defined.The naive lattice QFT discretization, which strongly parallels the Haar discretization, goes around this problem by adopting the ad hoc dictionary at resolution ∆x = a 0 /2 r , where a 0 is the reference lattice spacing for r = 0, and the discretized kinetic term coefficients are given by The problem is that the "Haar smearing" interpretation of ϕ r i as is inconsistent with the scaling dependence of the discretized kinetic term in (3.24).Specifically, consistency requires that coarse graining K r (through low-pass filter projections) should recover the kinetic term at lower resolutions.However, that is not the case.One can verify by inspection that which is the incorrect scaling.This problem is specific to wavelet classes whose basis elements are not differentiablenot a fatal flaw of Multiresolution Analysis in general.It can be bypassed by adopting a class of wavelets with continuous first derivatives.Indeed, Haar wavelets belong to a broader family of wavelets called Daubechies wavelets, named after their inventor Ingrid Daubechies [38].Other types of wavelets within the Daubechies family retain many of the simple and attractive properties of Haar wavelets while satisfying our differentiability requirement, as we shall discuss shortly.In what follows, briefly review of the definition and properties of Daubechies wavelets.

Daubechies DN wavelets
Daubechies wavelets provide a family of orthonormal bases for all square integrable functions on the real line, L 2 (R).Each type of Daubechies wavelets is characterized their number N of non-vanishing weights {h j }-hence the nomenclature Daubechies "DN " wavelets.The Haar wavelets, which have N = 2 non-vanishing weights h 0 = h 1 = 1/ √ 2 (see (3.11)), are also called Daubechies D2 wavelets.
The Daubechies DN basis elements are functions with compact support at different spatial resolutions, including basis elements which vanish outside an arbitrarily small open set.They are fractal functions, since any element of the basis can be generated from a single scaling function s(x) by applying a finite number of discrete translations and rescalings to s(x), and forming linear superpositions of these rescaled and translated functions.In particular, the DN scaling function satisfies a linear renormalization group equation: where the weights {h j } are a sequence of N non-vanishing numbers (to be discussed shortly), and the unit translation (T ) and dyadic rescaling (D) operators have been previously introduced in (3.2) and (3.3).Essentially, s(x) is the fixed point of an operation which takes the weighted average of a finite number of translated copies of itself and rescales it to half of its original support.The scaling function s(x) has compact support on the finite interval [0, N − 1].Further scaling functions with finer resolution are obtained by rescalings and translations of s(x): (3.28) Indeed, the subspace spanned by {s r n (x)}, S r , has a finest resolution scale of 1/2 r and contains all lower resolution subspaces S r−1 , S r−2 , ..., S 0 in the same nested form as in (3.15).
The wavelet subspaces are generated from the mother wavelet w(x), which also descends from the scaling function s(x) via: Above, the weights {g j } are related to the weights {h j } in (3.27) by: In full analogy with the case of Haar wavelets, a basis of translated and rescaled mother wavelets {w r n (x)} will span the DN wavelet subspace W r .A multiresolution decomposition of L 2 (R) with Daubechies DN wavelets then follows with the same structure of (3.21) and (3.22).
Additionally, the action of the low-and high-pass filters in (inverse) wavelet transforms for Daubechies DN wavelets is formally the same as for Haar wavelets: IWT with the difference that, above, the low-pass and high-pass filters L and H are determined by the DN weights that enter in (3.27) and (3.30): The weights {h j } for Daubechies DN wavelets are determined by three conditions: (1) Unit area for s(x): (2) Orthonormality of the basis {s n (x)}: (3) Ability of the basis {s n (x)} to locally represent k-degree polynomials with 0 ≤ k ≤ (N/2 − 1).In other words, there must exist a sequence {c n } such that: Given that wavelets and scaling functions are orthogonal to each other, the property above is equivalent to the requirement that the wavelets possess N/2 −1 vanishing moments: It then follows from (3.29) and (3.30) that this last condition on the weights {h j } can be written as: There are two solutions to conditions (3.34), (3.35), and (3.38) above, namely, {h j } and {h ′ j }.They are related by reverse ordering, i.e., h ′ j = h N −1−j , and their corresponding s(x) and s ′ (x) are mirror images of each other.
For our wMERA formulation, we would like to choose a Daubechies DN basis with a minimum requirement of a continuous first derivative.The smallest N value for which this condition is satisfied is N = 6 [39,40].Bases of smoother functions with higher differentiability could be chosen at the expense of additional computational cost, since such bases would have a larger number of weights {h j }.For our present purposes, we will consider Hamiltonians exhibiting at most first derivatives, and therefore, for the remainder of this paper, we will adopt Daubechies D6 wavelets for our wMERA formulation.
The weights for Daubechies D6 wavelets are given by: Figure 4 shows the D6 scaling function s(x) and mother wavelet w(x).
In anticipation of our wMERA study of free scalar and fermion theories in section 5, we end this section with a compilation of overlap integrals of D6 scaling and wavelet functions and their first derivatives, which will determine the kinetic-term coefficients of the discretized (1+1)-d Hamiltonians we will consider.For a derivation of the numerical values of the integrals listed below, see, e.g., [41].For scalar fields, there will be three types of kinetic-term coefficients, namely: Discrete translational invariance implies that these kinetic coefficients depend only on the separation between lattice sites, i − j.In addition, the compact support of s r n (x) and w r n (x) on the finite interval n/2 r ≤ x ≤ (n + 5)/2 r implies that the kinetic coefficients vanish for separations |i − j| ≥ 5. Defining the nonzero K j coefficients are given by: and illustrated in figure 5 (left).
The other kinetic coefficients can be obtained by applications of low-and high-pass filters (defined through (3.33), (3.39), and (3.30)): For a fermion field, the four types of kinetic coefficients will be: Unlike the scalar field case with symmetric K ss kinetic coefficients, in the fermionic case the D ss kinetic coefficients are antisymmetric.Defining the non-zero D j coefficients are given by and illustrated in figure 5 (right).The other fermionic kinetic coefficients can be obtained by applications of low-and high-pass filters: 4 Generic wMERA formulation The decomposition of the Hilbert space of a QFT into localized degrees of freedom organized by a hierarchy of scales makes the connection between wMERA and its discrete counterpart, MERA, more straightforward than the cMERA approach.In this section, we outline the generic architecture of wMERA using DWTs.To make the notation concrete, we assume a scalar field theory in (1+1) spacetime dimensions, but this assumption can be easily generalized to theories with more fields, either bosonic or fermionic.The generalization to higher (d+1) spacetime dimensions can also be made by considering wavelet bases in L 2 (R d ); this, however, is beyond the scope of this paper.
In the previous section, we introduced the resolution index r associated with a lattice spacing a 0 /2 r for the scaling and wavelet bases.In particular, under our unspecified choice for unit convention, the reference lattice spacing for r = 0 was set to a 0 = 1.Although the choice of a 0 is completely arbitrary, it is useful to identify it with some natural scale of the theory.For the sake of specificity, we will adopt a 0 = 1/m ϕ , where m ϕ is the mass of the scalar field in our toy example.Given this choice, we have that for an arbitrary resolution r, the "lattice site index" n of the scaling and wavelet functions represents a displacement from the origin by n integer steps of 1/(2 r m ϕ ).In addition, in order to preserve the dimensionless inner product of the scaling functions, we will also explicitly pull out a dimensionful normalization factor in their definition: The discretized degrees of freedom of a scalar field Φ(x) are then obtained from smearing the field in our chosen wavelet basis: The corresponding discretized degrees of freedom of the canonical conjugate momentum Π(x) are: They obey the usual canonical commutation relations, all other commutators vanishing.With the definitions above, the DWT representation of the scalar field Φ(x) at a fixed resolution (r+1) can be written as: and similarly for the canonical conjugate momentum Π(x): Alternatively, the discretized degrees of freedom can be recast in terms of creation and annihilation operators by defining: ) Above, σr and σr † are the creation and annihilation operators for "collective excitations" at scales 1/(2 r m ϕ ), whereas ωr and ωr † are the creation and annihilation operators for the corresponding "fluctuating" wavelet modes.They obey the commutation relations: with all other commutators vanishing.The width ∆ appearing in (4.7a), (4.7b) is arbitrary; it will be adjusted by the entangler when the wMERA parameters are variationally optimized.An obvious initial choice is ∆ = m ϕ .Note that the discretized degrees of freedom also obey (inverse) wavelet transform relations.Specifically, let {ŝ, ŵ} denote any one of the following pairs: { φ, φ}, {π, π}, {σ, ω}, or {σ † , ω † }.The (inverse) wavelet transform relations for each of these pairs are given by: IWT The D6-discretization of quantum fields also leads to discretized Hamiltonians at different resolutions.Consider, for the sake of illustration, a generic Hamiltonian for our toy example of a scalar field in (1+1)d: Plugging expressions (4.5a) and (4.6a) into (4.11) and integrating over x, we obtain the discretized system's Hamiltonian at resolution r in terms of scaling modes: where the kinetic term coefficients K ss have been defined in (3.43) and (3.44).
In (4.11), we left the interaction potential V Φ unspecified.Assuming V Φ = λ Φn , the corresponding discretized interaction potential in (4.12) will have the generic form: where The discretized Hamiltonian can also be written in terms of creation and annihilation operators of scaling modes by using (4.7a).Assuming ∆ = m ϕ for convenience, we have: Finally, through the use of IWTs defined in (4.10)-or, alternatively, by using the field decompositions in (4.5b) and (4.6b) combined with (4.7)-the discretized Hamiltonian can also be cast in terms of scaling and wavelet modes.Concretely, where the kinetic term coefficients K sw and K ww have been defined in (3.46) and (3.47), respectively.
The wavelet decomposition of the field degrees of freedom also implies a decomposition of the Hilbert space of this theory into orthogonal subspaces at different scales: Above, H r is the Hilbert subspace containing the "long-range" scaling modes describing collective excitations of the scalar field over scales 1/(2 r m ϕ ).It is identified with the Fock space generated by creation operators σr † acting on the product state defined by: Likewise, the finer resolution Hilbert subspaces h r ′ (r ′ ≥ r) contain modes describing shortrange quantum fluctuations of the scalar field at scales 1/(2 r ′ +1 m ϕ ).Each of these subspaces is identified with the Fock space generated by the creation operators ωr ′ † acting on the product state defined by: With the Hilbert space decomposition described in (4.17)-(4.19),we can sketch the main ingredients of wMERA: (1) At very low resolutions r IR ≪ 0, the discretized degrees of freedom describe fluctuations at lengths much longer than the correlation length of the system (that is, L ≫ 1/m ϕ in our toy example).At such scales, the nonlocal terms in the Hamiltonian-including the kinetic term-are suppressed, and the ground state exhibits exponentially negligible spatial entanglement.Hence, the far infrared ground state is well-approximated by a product state ansatz6 , At resolutions close to the correlation length of the system, (that is, r ∼ 0 and L ∼ 1/m ϕ in our toy example), entanglement is still relatively localized (i.e., it is only non-negligible between neighboring sites, and decays exponentially for site separations |n−m| ≫ 1).Hence, any resolution close to r = 0 is a good starting point for building the wMERA.Choosing, for concreteness, r = 0 as the initial resolution, the ground state can be described by the following ansatz state in the Hilbert space H 0 : where Û0 is a quasi-local unitary entangler ansatz.The (a priori unkown) parameters of Û0 are obtained variationally by minimizing the expectation value of the discretized Hamiltonian at r = 0: This constitutes the first wMERA layer.
(2) To build the second layer, the first step is to enlarge the Hilbert space to include the finer resolution wavelet subspace h 0 via the fine-graining mapping W 1 : Accordingly, the ansatz is augmented by wavelet states initialized in a product state with the rest of the system.By abuse of notation, we write: (3) The second step is to entangle the newly introduced wavelet states with the original, lower-resolution scaling states.This is parameterized by a quasi-local entangler V1 : Similarly to the first layer, the parameters of the entangler ansatz V1 are variationally optimized by minimizing the expectation value of the discretized Hamiltonian at r = 1: where, above, the r = 1 Hamiltonian Ĥ1 σω is expressed in terms of scaling and wavelet modes as in (4.16), and (4) The state in (4.25) is the ansatz for the ground state at resolution r = 1 defined in the Fock space of r = 0 scaling and wavelet states.As we will justify shortly, at this point it is convenient to perform an inverse wavelet transform T1 −1 to re-express this ansatz as a state in the Fock space of r = 1 scaling states: Note that the variational optimization of the entangler V1 can be alternatively performed after the IWT step is applied.Noting that the r = 1 Hamiltonians Ĥ1 σ and Ĥ1 σω (generically defined in (4.15) and (4.16)) are related by and defining we can alternatively perform the variational optimization of the r = 1 entangler parameters by minimizing the following expectation value: (5) Finally, these operations can be performed iteratively to construct the layers of the wMERA until one reaches either a fixed point, or the desired UV resolution r UV : This state describes the vacuum of the theory at scale r UV provided that the parameters of the transformations Vr are, at each layer, chosen so as to minimize the expectation value of the system's Hamiltonian Ĥr .If these conditions are not sufficient to uniquely determine the ansatz parameters, one can also demand that expectation values of other physical observables are preserved under RG flow, such as n-point correlation functions with n > 2.
We expect that the variational optimization of the wMERA parameters in (4.32) can be efficiently implemented numerically with unitary ansatzes Vr that are quasi-local and sparse.
If this is indeed the case, then computational techniques for the variational optimization of MERA tensor networks developed for discrete systems should be adaptable to wMERA for studying QFTs.
In the next section, we will show that quasi-local and sparse wMERA unitaries are indeed sufficient to describe ground state entanglement in two simple examples: a free scalar field and a free Dirac fermion in (1+1) spacetime dimensions.
We end this section by justifying the need for inverse wavelet transforms at every resolution layer in the wMERA formulation.Note that if these transformations were not performed, the description of the system would incur two problematic features: (1) a nonuniform lattice structure, and (2) discretized Hamiltonians Ĥr becoming increasingly nonlocal as the number of wMERA layers is increased.To illustrate this situation, consider a resolution layer r = r 0 + n, where n ≫ 1 and r 0 is the initial wMERA layer.In the absence of IWTs, the degrees of freedom at this layer would consist of the scaling modes, ϕ (r 0 ) , and n + 1 types of wavelet modes: φ (r 0 ) , φ (r 0 +1) , . . ., φ (r 0 +n) .These n + 2 different types of modes would have different lattice densities, with ϕ (r 0 ) and φ (r 0 ) having the sparsest lattices (with lattice spacing a 0 ), and φ (r 0 +n) having the densest lattice (with lattice spacing a 0 /2 n ). Figure 6.Depiction of the lattice describing the system at a deep UV layer if inverse wavelet transforms were not incorporated in the wMERA architecture.The degrees of freedom of the first (bottom layer) consist of scaling modes (depicted in yellow) and wavelet modes (depicted in dark blue).The next layers would iteratively introduce shorter-and shorter-distance wavelet modes (depicted in dark blue) with successively denser lattices.In this lattice architecture, a scaling mode in the first layer would couple to O(n 2 ) wavelet modes, where n is the number of wMERA layers that support this particular resolution.
In addition, any given scaling mode ϕ (r 0 ) j would couple to a total of 5n 2 + 23n + 36 /2 other modes, namely, In other words, the nonlocality of the Hamiltonians Ĥr 0 +n would grow quadratically with the number of wMERA layers.This situation is depicted in figure 6.It is easy to see that not only this lattice structure becomes quickly impractical, but the ansatz would require unitary entanglers that more and more nonlocal at each layer.The IWTs address these difficulties by redefining the degrees of freedom when moving from one layer to the next so that the uniformity of the lattice structure is preserved, and the quasi-locality of the discretized Hamiltonians Ĥr is restored at each layer.

Concrete wMERA Examples
For the remainder of this article, we work through the concrete implementation of wMERA for a couple examples of free field theories in (1+1) spacetime dimensions: a scalar field and a Dirac fermion.For each example, we begin by deriving the exact entangler directly in the continuum without invoking cMERA complications such as ad hoc UV cut-offs or path-ordered integrals over scaling transformations.We then solve these theories using our proposed variational wMERA approach: we obtain the optimized numerical results for the wMERA entanglers, which we then use to compute the two-point correlation functions for these theories for both spacelike and timelike separations.

Free theory of a scalar field in (1+1)d
Consider a free scalar field theory in (1+1)d.Its Hamiltonian is given by where the scalar field Φ(x) its canonical conjugate momentum Π(x) obey the equal-time commutation relations Φ(x), Π(y) = i δ(x − y); all other equal-time commutators vanishing.
Since free theories only have Gaussian correlations, it should come as no surprise that their entanglers are Gaussian7 (i.e., the entanglers' exponents are quadratic functions of Φ and Π).Specifically, writing the field and its canonical conjugate momentum in terms of position-space creation and annihilation operators: we would like to find the entangler Û that relates the product state P to the ground state 0 : where the product state P has, by definition, no entanglement in position space and can be thought of as the generalization of the state in (4.18) to the continuum: â(x) P = 0 ∀ x. (5.5) In this case, the entangler ansatz is a generator of Bogoliubov transformations (a.k.a.multimode squeezing), and has the generic form: or, alternatively in terms of position-space creation and annihilation operators, The entangler exponent ζ(x − y) defines the entangling parameter that needs to be variationally optimized.The easiest way to proceed is to go to Fourier space by defining: (5.9) Note that â † k and âk in (5.8b) are not the creation and annihilation operators of momentummode Hamiltonian eigenstates.In particular, they do not annihilate the ground state, but, because of (5.5), they do annihilate the initial product state: (5.10) In Fourier space, the entangler and the Hamiltonian are given by and where ω 2 k has the usual definition ω 2 k ≡ k 2 + m 2 ϕ .Recalling the commutation relations âk , â † q = 2π δ(k − q), one can easily derive the action of the entangler on the creation and annihilation operators, and the entangler-transformed Hamiltonian, Obtaining the vacuum expectation value of the Hamiltonian is now straightforward.Using (5.14) and (5.10), we have: The final step is to minimize the vacuum expectation value of the Hamiltonian with respect to entangler exponent ζ k :

.16)
We can now use (5.9) and (5.16) to obtain the exact entangler exponent in position space, (5.17) It is also informative to look at the entangler-transformed Hamiltonian using the exact entangler in (5.16) and (5.17): where, above, (5.19) While the use of FTs in this exercise was convenient to obtain the exact entangler, this trick unfortunately only works well in the case of free theories for which the entangler is Gaussian (as we argued extensively in sections 1 and 2).For interacting theories, we expect that working directly in position space through the use of DWTs will make calculations much more tractable.For the remainder of this subsection, we show how this actually works by solving this theory using wMERA.
We work with the D6-discretized degrees of freedom defined in (4.2), (4.3), and (4.7), whose dynamics are governed by the discretized Hamiltonians (4.15) and (4.16) with vanishing interaction potential (i.e., V = 0).To recall our notation, at any given resolution r the product state P σ r is related to the ground state by a unitary entangler Ûr : 0 σ r = Ûr P σ r . (5.20) In our free theory example, the ansatz for Ûr is a generator of Bogoliubov transformations between scaling modes: (5.21) The entangling exponent ζ (r) is a Hermitian matrix whose entries satisfy the following relation due to discrete translational invariance.The entangler Ûr transforms scaling modes in the following way: (5.23b) The other type of entangler discussed in section 4 was the unitary operator Vr that entangled scaling and wavelet modes at the wMERA layer with resolution r.While Vr is still a Gaussian entangler, it has a more general parametrization than Ûr .Concretely, , (5.25) or one can perform the optimization after the IWT reorganizes the scaling and wavelet modes σ (r−1)( †) and ω (r−1)( †) into σ r( †) scaling modes.While this choice is completely up to the developer of the optimization algorithm, this author found the second option easier to implement numerically, for two main reasons.First, there is a significant amount of redundancy in the Vr ansatz parameters ν , one can easily anticipate that these transformations will have a more complicated form compared to (5.23).
Hence, we choose to proceed with the wMERA optimization by instead minimizing the parameters of each layer after the IWT is performed.That is, we will optimize Ĥr σ Ûr P σ r . (5.26) Note, however, that this does not mean we intend minimize all ζ (r) ij parameters of the ansatz in (5.21) at once, since this would not take any input from previous wMERA layers, and would lead to challenges with numerical efficiency and accuracy.Instead, the crucial point to note is that the "long range" entanglement encoded in ζ (r) ij has already been optimized in the previous wMERA layers.The way to take advantage of that is through the following: using WTs, we can recast Ûr in terms of scaling and wavelet degrees of freedom8 ,

Tr Ûr Tr
Above, the ξ (r−1) exponents are the WTs of the entangling exponent ζ (r) , (5.28c) Note that they too have a restricted form due to discrete translational invariance: (5.29) The IWT relating ζ (r) to ξ (r−1) is given by This expression is key for our algorithm to minimize Ĥr σ .Specifically, at a fixed wMERA layer with resolution r, we treat ξ ww(r) i−j and ξ sw(r) i−j as free parameters, and fix the condition where, above, we have used the short-hand notation Using (5.32) and (4.18), we then have 2ζ . (5.34) As it turns out, the minimization of this vacuum expectation value can be performed analytically.Differentiating (5.34) with respect ζ (r) and setting it to zero, we obtain whose solution is given by (5.36) However, as previously stated, our actual goal in this demonstration of the wMERA algorithm is to minimize (5.34) numerically using the parametrization in (5.30) and the input in (5.31); we will only use the exact solution for ζ (r) ij in (5.36) to validate our numerical results.We can further manipulate (5.34) using (5.33), (3.43), and (3.44) to obtain: (5.37) The expression above has two remarkable aspects.First, even though it is a trace over a square matrix of infinite dimension, the infinite summation over one of the indices factors out (this is due to discrete translational invariance).Second, the r.h.s. of (5.37) depends only on a small, finite number of elements of the matrix exponential of ζ (r) (this is due to the Daubechies discretization, whose scaling and wavelet functions have compact support in position space).In particular, because of the form of (5.37), there is no need to perform a finite volume truncation of the system in order to minimize Ĥr σ with respect to the variational parameters.However, a truncation is still needed in the number of variational parameters, of which there is in principle an infinite number.One may question whether any such truncation would preserve the properties of the ground state.To argue that this is the case, note that the matrix ζ (r) i−j is the discrete counterpart of the entangler in the continuum, ζ(x − y), whose exact solution is given in (5.17).In particular, ζ(x − y) becomes exponentially suppressed for separations larger than the correlation length of the system, i.e., for |x − y| ≳ m −1 ϕ .Therefore, we expect that the entries ζ (r) i−j will likewise be exponentially suppressed for site separations |i − j| ≳ 2 r .
Recalling that the free parameters of ζ (r+1) are given by ξ ww(r) and ξ sw(r) as defined in (5.30), in practice our truncation of free parameters will amount to setting where i (r) max and j (r) max are finite integers which can be increased to improve numerical accuracy.Remarkably, we will see shortly that i (r) max and j (r) max do not need to grow with resolution r, and can be kept small (i.e., in the range of ∼ 3-5) for all wMERA layers while still retaining excellent accuracy 9 .This is in contrast with the truncation of entangling parameters ζ (r) i , which must grow with r to allow the same desired degree of accuracy to be retained, as we just argued in the previous paragraph., defined in (5.21) and (5.28), for a few illustrative choices of resolution r.We recall that the lattice spacing is related to the resolution index r via a r = 1/(2 r m ϕ ).
In appendix A, we discuss a method to obtain the entries of Exp 2ζ j as a function of a truncated array ζ i .This, combined with (5.30), can plugged into (5.37)for the numerical variational optimization.
The exact solutions for ζ (r) i−j , ξ ww(r) i−j , and ξ sw(r) i−j are shown in figure 7 for a few illustrative choices of resolution r.The plots clearly show, as previously stated, that the range of the entangling parameters ξ ww(r) and ξ sw(r) does not grow exponentially with increasing resolution.Indeed, their optimal values not only quickly saturate at resolutions close to r = 0, but their range decays exponentially with increasing lattice site separations, which justifies the truncations in (5.38).This is in stark contrast to the scaling entangling parameters ζ (r) , whose range grows exponentially with resolution.
The results of our numerical wMERA optimization show excellent agreement with the exact solutions within the truncated range of parameters in our algorithm.Specifically, in obtaining our results, we truncated the ξ ww i , ξ sw j , and ζ k arrays at i max = 3, j max = 4, and k max = 21, respectively, in all the layers of our calculation, which ranged from r = −3 to 10.
For a fixed resolution layer r, our numerical results for ζ  ij also allow us to reconstruct the scalar field two-point function in the continuum.For spacelike separations, it suffices to compute the equal-time twopoint function; the answer for any other spacelike interval can be obtained by a Lorentz transformation.Recalling the DWT representation of the scalar field Φr (x) at a fixed resolution r in (4.5a), we have: where, in the last equality, we used (4.1), (4.7a), (4.18), and (5.)), calculated using the numerical output of our wMERA optimization.The solid black curve shows the exact continuum solution for ζ(x − y) (see (5.17)).Bottom: the reconstructed continuum scalar twopoint function at equal times, Φ(x) Φ(y) , calculated using the numerical output of our wMERA optimization.The exact continuum solution for Φ(x) Φ(y) is given by the solid black curve (see (5.47)).In both plots, the functions were calculated for different ranges in m ϕ |x − y| using different wMERA layers, indicated by the resolution index r.We recall that the lattice spacing is related to the resolution index r via a r = 1/(2 r m ϕ ).
Our final task in solving this free scalar theory with wMERA is to obtain the scalar twopoint function at finite-time separations.This requires implementing real-time evolution, which is also straightforward in the wMERA framework.The first step is to obtain the entangler-transformed Hamiltonians at different resolutions using the variational solution for the entangler parameters.Plugging (5.36) into (5.32),we obtain: where (5.42) We can also define the entangler-transformed time-evolution operator10 : Its action on the scaling modes at resolution r is easily derived: We now have all the ingredients to obtain the scalar field two-point function at resolution r for finite-time separations : Φr (x, t) Φr (y, 0) = 0 σ e i t Ĥr σ Φr (x) e −i t Ĥr σ Φr (y) 0 σ r (5.45) The expression above, which generalizes (5.40), is restricted to the range |t Piecing together Φr (x, ∆t) Φr (x, 0) at different resolutions, we recover the continuum result for Φ(x, ∆t) Φ(x, 0) at timelike separations over several orders of magnitude in ∆t, as shown in figure 9.
It is important to note that our wMERA results for Φr (x, t) Φr (y, 0) should track the continuum solution over scales longer than the resolution r, but otherwise show some degree of unphysical structure at scales of order |x − y| ∼ (2 r m ϕ ) −1 .Indeed, in figures 8 and 9, we chose to plot our results for a discrete array of points {x i , y i } satisfying With these definitions, we can rewrite (5.45) as: k .
The unphysical structure of this reconstructed two-point function at scales |x i − y i | ∼ (2 r m ϕ ) −1 can be seen by varying the parameter ϵ above over its range of 0 ≤ ϵ < 1 (in other words, varying ϵ will probe the lack of smoothness of the D6 scaling function).Hence, Φ(x, t) Φ(x, 0) , calculated using the numerical output of our wMERA optimization.The real and imaginary parts are shown in the top and bottom plots, respectively.The apparent derivative discontinuity at m ϕ t = 1 in the bottom plot is an artifact of the hybrid scale of the horizontal axis: note that the scale is logarithmic for m ϕ t < 1, and linear for m ϕ t > 1.This choice was made so that the comparison between the numerical results and the exact continuum solution (solid black curves) is displayed with better clarity.In both plots, the functions were calculated for different ranges in m ϕ t using different wMERA layers, indicated by the resolution index r.We recall that the lattice spacing is related to the resolution index r via a r = 1/(2 r m ϕ ).
for every plot in this paper 11 showing continuum functions reconstructed from numerical wMERA results, we select the value of ϵ that minimizes these unphysical fluctuations.We can compare our numerical wMERA results for the scalar two-point function with the exact analytical expression for the Feynman propagator in position space, given by: where t > 0, Θ is the Heaviside step function, H 0 is a Hankel function of the second kind, K 0 is a modified Bessel function of the second kind, and s is the spacetime interval, i.e., s 2 ≡ t 2 − x 2 .
As shown in figures 8 and 9, the agreement between our numerical results and the exact continuum solution is quite satisfactory considering our aggressive truncation choices.That said, one can see that the accuracy and precision of our results slowly degrade with increasing r.This is not a shortcoming of wMERA itself, but an expected artifact of the approximations made in our numerical code.Specifically: (i) we kept the ζ |k| coefficients to 3 (see appendix A for details); and (iii) we truncated the Exp 2ζ (r)  l array to l (r) max = 36 for all layers.These choices, which were not optimized, were made so we could run our code on a personal laptop with a dual Intel i5 core processor in a couple of hours (the code was written in Wolfram Mathematica [43]).However, with more computational resources, one can relax the truncation choices to overcome the degradation in accuracy and precision with increasing resolution.

Free theory of a Dirac fermion in (1+1)d
Our second and final wMERA example is a free theory of a Dirac fermion in (1+1) spacetime dimensions.Many steps of this exercise (and their justifications) have a strong parallel with the scalar field theory case just considered.Therefore, in this section, we will skip over details that have already been dissected in the previous section.
We start by laying out our notation.We choose conventions such that the Dirac spinor and γ-matrices are given by: (5.48) The system's dynamics is determined by the Hamiltonian and equal-time anticommutation relations all other equal-time anticommutators vanishing.We define the initial product state upon which the wMERA will be built as the unentangled state in position space that is annihilated by χ and ψ: χ(x) P = ψ(x) P = 0 ∀ x. (5.51) Our goal is to find the entangler Û that relates the product state P to the Hamiltonian ground state 0 : 0 = Û P . (5.52) The ansatz for Û is a Gaussian operator (serving the analogous role for fermions that the generator of Bogoliubov transformations does for scalars [44]).Its parametric form in generic (d+1)-spacetime dimensions (with d odd) is given by: (5.54) A key difference compared to the scalar case is that, for fermions, the entangling exponent in (5.54) is anti-symmetric, that is, To proceed with the variational optimization of the entangler, we go to Fourier space.Defining and we can write the Hamiltonian as: (5.57) The entangler in Fourier space is given by: where, above, we used the fact that The action of the entangler on the spinor components is easily derived 12 : (5.59b) 12 Recalling that, in Fourier space, the anticommutation relations are given by χk , χ † q = 2π δ(k − q) and ψk , ψ † q = 2π δ(k − q).
With (5.57) and (5.59), we obtain the entangler-transformed Hamiltonian, Using (5.60) and (5.51), we then obtain the vacuum expectation value of the Hamiltonian, (5.61) Minimizing ⟨ Ĥ⟩ with respect to entangler exponent ζ k , we get: (5.62) Using (5.56), we finally obtain the fermionic entangler exponent in position space, Finally, we can plug in the exact entangler expression into (5.60) to write the entanglertransformed Hamiltonian as = dx dy Ω(x − y) χ † (x) χ(y) − ψ(x) ψ † (y) (5.64b) = dx dy Ψ(x)Ω(x − y) Ψ(y). (5.64c) where Ω(x − y) has already been defined in (5.19).We now move on to the discretized theory to build the wMERA.Using the Daubechies D6 wavelet basis 13 to smear the fermion fields, we define: so that the fermionic field representation at resolution r is given by14 : (5.66) Plugging (5.66) into (5.49),we obtain the discretized Hamiltonian at resolution r: where the kinetic term coefficients D ss ij have been defined in (3.52) and (3.53).The ansatz for the unitary entangler that transforms the initial product state into the ground state of the discretized Hamiltonian in (5.67) can be written as: where entangling exponent ζ (r) is an anti-Hermitian matrix, and the operator Rr (θ), with a single additional variational parameter θ, is given by:

.69)
This operator mixes the two components of the Dirac spinor locally.It turns out that, since we are working with the Dirac representation of γ-matrices, θ = 0. Had we chosen to work with the Weyl representation, we would find that θ = π/4.For simplicity, in our wMERA exercise we will set θ = 0 from the start, although including θ as a ansatz parameter in our wMERA optimization would have been straightforward and would not have affected the final results.Setting Rr (θ) = 1 in (5.68), the action of the entangler Ûr on the fermionic scaling modes is given by: (5.70b)As with the scalar example in subsection 5.1, given a wMERA layer r, our free variational parameters ξ ww(r−1) and ξ sw(r−1) will be embedded in ζ (r) via an IWT: while ξ ss(r−1) will be fixed by the result of the optimization of the previous layer through (5.72) Similarly to the scalar example, the entangling matrices ζ (r) , ξ ww(r) and ξ sw(r) have their forms restricted due to discrete translational invariance.However, unlike the scalar example, ζ (r) and ξ ww(r) are anti-Hermitian: (5.73) Plugging (5.76) into (5.74),we obtain the entangler-transformed Hamiltonian in terms of the optimized entangler parameters: with (5.81) The entangler-transformed time-evolution operator is then given by: (5.82) It transforms the fermionic scaling modes as : (5.83b) Using (5.70) and (5.83), we can obtain the fermionic two-point correlation functions at any resolution r.For instance, χr (x, t) χr † (y, 0) = 0 e i t Ĥr σ χr (x) e −i t Ĥr σ χr † (y) 0 r (5.84) Analogous steps can be followed to obtain As usual, the expressions above are restricted to their regimes of validity, |t Piecing together the reconstructed two-point functions (5.84)-(5.86) at different resolutions, we can compare them with the exact continuum solution for Ψ(∆x, ∆t) Ψ(0, 0) in (5.79) at both spacelike and timelike separations over several orders of magnitude in ∆x .The reconstructed continuum fermionic two-point function at equal times calculated using the numerical output of our wMERA optimization.The top and bottom plots show, respectively, Tr Ψ(x) Ψ(y) = χ(x) χ † (y) − ψ † (x) ψ(y) and Tr i γ 1 Ψ(x) Ψ(y) = i χ(x) ψ(y) + i ψ † (x) χ † (y) .The exact continuum solution is given by the solid black curves (see (5.79)).In both plots, the functions were calculated for different ranges in m ψ (x − y) using different wMERA layers, indicated by the resolution index r.We recall that the lattice spacing is related to the resolution index r via a r = 1/(2 r m ψ ). and ∆t, respectively.Our results are shown in figures 10, 11, and 12.We find very good agreement between the wMERA and exact solutions for spacelike separations (see figure 10).For timelike separations, there is good agreement at a resolutions comparable to or lower than the system's correlation length, m −1 ψ .However, at shorter time scales, there are two competing effects worth mentioning: the first is a loss of accuracy and precision stemming from the truncation choices in our numerical algorithm.We saw this effect in the scalar case as well, and it can be easily mitigated with modified truncations choices.The second effect is an oscillatory behavior of the wMERA two-point function around the The reconstructed continuum fermionic two-point function χ(x, t) χ † (x, 0) for timelike separations calculated using the numerical output of our wMERA optimization.The real and imaginary parts are shown in the top and bottom plots, respectively.The oscillatory behavior of the wMERA results around the exact continuum solution, given by solid black curves, stems from presence of fermion doublers in the spectrum of the discretized theory-this effect becomes significant at short time scales m ψ t < 1.In both plots, the functions were calculated for different ranges in m ψ t using different wMERA layers, indicated by the resolution index r.We recall that the lattice spacing is related to the resolution index r via a r = 1/(2 r m ψ ).
exact solution (see figures 11 and 12).This effect is actually physical, and it reflects the fact that the continuum and discretized theories have a different spectrum of states-an effect commonly know as the "fermion doubler problem" in discretized field theories with fermions.Indeed, our Daubechies D6 discretization is subject to the Nielsen-Ninomiya theorem [45,46].This means that, in the massless limit, our D6-discretized theory of Dirac fermions has indeed twice the number of zero-modes compared with the continuum theory [47].This explains the presence of oscillatory features in our wMERA results: since the Hamiltonian mass term becomes less important as the system flows to the UV, the system's correlation function should approach the massless solution at very high resolutions; therefore, it should not be surprising that the effects of the doubler states become more and more pronounced at shorter and shorter time scales.
It is notable, however, that modulo these UV oscillations due to doublers, the two-point function for D6-discretized theory actually tracks the shape of the continuum solution, a feature that is not shared by the naive lattice field theory discretization based on Haar wavelets.Indeed, for the sake of comparison, we also obtained the wMERA solution for the Haar-discretized fermionic theory, and found that (i) the resulting two-point function had a completely different shape than the continuum solution-this also follows from solving the Haar-discretized theory using other methods, such as discrete Fourier transforms; and (ii) the real-time Hamiltonian evolution of the Haar-discretized theory suffered from much higher noise than the real-time evolution of the D6-discretized theory.This noise was present at all scales and was not mitigated by increasing the truncations orders in the numerical algorithm.This seems to indicate that the D6 discretization might be a better choice than the traditional, Haar discretization of fermionic theories.It is not clear to this author, however, whether this preliminary conclusion will hold once these theories are properly regularized to remove the doubler states.
Given the ubiquity of the fermion doubler problem in lattice field theory, many methods have been developed, in various contexts, to rid the spectrum of fermionic doubler states (see, e.g., [48]).At least some of these methods should be adaptable to the Daubechiesbased discretization.Exploration of these possibilities, however, is deferred to future work.

Conclusions and Outlook
Traditionally, the concept of entanglement has been completely obscured in the analytical and numerical treatments of QFTs.More recently, entanglement has emerged as an important element in the context of the holographic principle, black holes, and gravity [49][50][51][52][53], and as a useful measure to quantify complexity and decoherence in quantum field theories [54][55][56][57][58][59].In this article, we hope to have conveyed another important role that entanglement can play, namely, as a computational tool 15 to develop new numerical methods for quantum field theories, which could potentially lead to novel insights.
In this study we introduced wMERA, an adaptation of entanglement renormalization for quantum field theories that preserves several key features of the original MERA, namely, (i) discretized degrees of freedom localized in position space; (ii) a discrete, layered architecture of the ansatz paralleling an inverse RG flow; and (iii) the quasi-locality of entanglers at each layer.We considered two simple free field theory examples for which we demonstrated the concrete implementation of the wMERA algorithm in detail.In particular, we showed that wMERA allowed us to fully extract the two-point correlators for spacetime separations spanning several orders of magnitude.Notably, we were able to efficiently perform real-time evolution over many decades in ∆t by working in the appropriate resolution for the ansatz and for the system's Hamiltonian.In other words, the use of the RG flow played a key role in the efficiency of these computations.
In these specific wMERA examples, however, we did depart from the original MERA implementation in one respect: instead of having the entanglers acting on the Hilbert space (à la Schrödinger picture), we chose to have the entanglers transform the field operators (à la Heisenberg picture).This, combined with translational invariance, allowed us to preserve the systems' infinite length (that is, we were able to bypass finite "volume" artifacts).In addition, we were able to preserve the infinite dimension of the local Hilbert spaces for bosonic degrees of freedom.Whether this "Heisenberg picture" approach will be easily generalizable to wMERA implementations for interacting field theories remains to be seen; either way, efficient numerical MERA algorithms in the "Schrödinger picture" have been developed for condensed matter systems, and there is no reason to believe that they could not be successfully adapted for (suitably regularized and discretized) relativistic field theories.
Our motivation for proposing wMERA is obviously not to have an alternative method for solving free field theories.The real goal is to use this approach to explore new variational methods for solving challenging problems in quantum field theory (such as strongly coupled dynamics) as well as the incorporation of RG flow in quantum algorithms (to, e.g., improve their efficiency and accuracy, including in the preparation of initial states, in the simulation of real-time evolution, and in the readout and interpretation of the simulation's output.) The next important step towards this goal is explore non-Gaussian entanglers for interacting QFTs, and to understand their properties in different dynamical regimes.Possible avenues would be to explore hybrid entangling layers composed of the intercalation of ultralocal non-Guassian entanglers with quasi-local Bogoliubov transformations, which parallels the perturbative calculation approach; and, for systems that exhibit dualities or confinement, to explore the introduction of "ancillary" degrees of freedom to represent the dual or confined states, and devise entanglers whose action is to map the original states to the dual or confined states.Another important next step would be to generalize the wMERA formulation to (3+1)d QFTs by exploring optimal wavelet classes in 3 spatial dimensions.In this exploratory phase, the wMERA results would have to be validated against standard perturbative calculations in the weakly coupled regime, or against exact solutions, when available.Importantly, the matrices M (k) satisfy the following relations: , M (j) = 0 ∀ i, j and M (i) M (j) = M (i+j) .(A.3) These relations allow us to trade the basis of matrices M (k) for powers of a commuting variable (a c-number) denoted by m.That is, in (A.2), we make the replacement

Figure 3 .
Figure 3. Left panel: a pictorial illustration of the Haar wavelet transform, which decomposes a signal into low-resolution "collective" features (green scaling functions) and high-resolution features (orange wavelet functions).Right panel: the inverse Haar wavelet transform, which recombines this decomposition into a basis of finer-resolution collective features.
depend on a smaller number of parameters, as it should became clear shortly).Second, while we have not shown explicitly how the Vr ansatz in (5.24) transforms the scaling and wavelet modes σ

. 31 )
that is, we use as input the numerical values for ζ (r−1) obtained from the optimization of the previous layer.While in the exact solution the relation(5.31)  is only approximately true, it is quite accurate for lattice site separations |i − j| ≳ 5.For shorter lattice site separations |i − j| ≤ 4, the optimization of the free parameters ξ ww(r) and ξ sw(r) can compensate for any inaccuracy introduced by(5.31).With the parametrization of the entanglers now fully specified, we can proceed to write the vacuum expectation value to the Hamiltonian in(5.26)  in terms of ζ (r) ij .First, we obtain the entangler-transformed Hamiltonian in (4.15) (with V = 0) using (5.23):

Figure 7 .
Figure 7.The optimized values of the entangling parameters ζ (r) ij , ξ ww(r) ij , and ξ sw(r) ij ij allow us to recover the continuum result for the entangling exponent ζ(x − y) restricted to the range |x − y| ∼ O(1 − 10) × (2 r m ϕ ) −1 : −j .Piecing together the functions ζ (r) (x − y) at different resolutions, we can reconstruct the continuum entangling exponent ζ(x − y) over several orders of magnitude in |x − y|, as shown in figure 8. Our results for ζ (r)

Figure 8 .
Figure 8. Top: the reconstructed continuum entangling exponent ζ(x − y) (see (5.7)), calculated using the numerical output of our wMERA optimization.The solid black curve shows the exact continuum solution for ζ(x − y) (see(5.17)).Bottom: the reconstructed continuum scalar twopoint function at equal times, Φ(x) Φ(y) , calculated using the numerical output of our wMERA optimization.The exact continuum solution for Φ(x) Φ(y) is given by the solid black curve (see(5.47)).In both plots, the functions were calculated for different ranges in m ϕ |x − y| using different wMERA layers, indicated by the resolution index r.We recall that the lattice spacing is related to the resolution index r via a r = 1/(2 r m ϕ ).

Figure 10
Figure 10.The reconstructed continuum fermionic two-point function at equal times calculated using the numerical output of our wMERA optimization.The top and bottom plots show, respectively, Tr Ψ(x) Ψ(y) = χ(x) χ † (y) − ψ † (x) ψ(y) and Tr i γ 1 Ψ(x) Ψ(y) = i χ(x) ψ(y) + i ψ † (x) χ † (y) .The exact continuum solution is given by the solid black curves (see(5.79)).In both plots, the functions were calculated for different ranges in m ψ (x − y) using different wMERA layers, indicated by the resolution index r.We recall that the lattice spacing is related to the resolution index r via a r = 1/(2 r m ψ ).

Figure 11 .
Figure 11.The reconstructed continuum fermionic two-point function χ(x, t) χ † (x, 0) for timelike separations calculated using the numerical output of our wMERA optimization.The real and imaginary parts are shown in the top and bottom plots, respectively.The oscillatory behavior of the wMERA results around the exact continuum solution, given by solid black curves, stems from presence of fermion doublers in the spectrum of the discretized theory-this effect becomes significant at short time scales m ψ t < 1.In both plots, the functions were calculated for different ranges in m ψ t using different wMERA layers, indicated by the resolution index r.We recall that the lattice spacing is related to the resolution index r via a r = 1/(2 r m ψ ).
This property enables the matrix ζ to be represented by an infinite array of parameters ζ k through the following decomposition: