Rigorous free fermion entanglement renormalization from wavelet theory

We construct entanglement renormalization schemes which provably approximate the ground states of non-interacting fermion nearest-neighbor hopping Hamiltonians on the one-dimensional discrete line and the two-dimensional square lattice. These schemes give hierarchical quantum circuits which build up the states from unentangled degrees of freedom. The circuits are based on pairs of discrete wavelet transforms which are approximately related by a"half-shift": translation by half a unit cell. The presence of the Fermi surface in the two-dimensional model requires a special kind of circuit architecture to properly capture the entanglement in the ground state. We show how the error in the approximation can be controlled without ever performing a variational optimization.


I. INTRODUCTION
Characterizing quantum phases of matter and computing their low-temperature physical properties are two of the central challenges of quantum many-body physics. Part of the challenge is that many phases and phase transitions are best characterized in terms of their pattern of entanglement, as opposed to, say, their pattern of symmetry breaking [1]. As an extreme case, topological phases of matter (see Ref. [2] for a review) are solely characterized by entanglement, and since a topological state is gapped, it is insensitive to any local perturbation [3,4] and need not break any symmetry. By contrast, metallic states with Fermi surfaces have many low-energy excitations, but they can also be usefully characterized in terms of their anomalously high entanglement. In this work, we are concerned with such metallic states.
A useful idealization is to focus on ground-state physics, where the entanglement entropy of spatial subsystems provides powerful nonlocal diagnostics of phases. For example, most known ground states obey an "area law": The entanglement entropy of a subsystem scales as the boundary size of the subsystem [5]. Topological states have a subleading (in subsystem size) shape-independent "topological entanglement entropy" term that partially characterizes the state [6,7]. Another diagnostic is the scaling of entropy with subsystem size itself: Metals logarithmically violate the area law when they have a Fermi surface [8][9][10].
However, discussions based on entanglement entropy are only the beginning. We can more fully characterize the pattern of entanglement in a quantum state by giving a recipe for physically constructing it from unentangled degrees of freedom (d.o.f.). Such a recipe is called a quantum circuit and constitutes a sequence of local unitary transformations that produce the desired state from an unentangled state. While symmetry-breaking states, e.g., a ferromagnet with all spins aligned, can be caricatured by an unentangled state (or, more realistically, by a state in which only nearby d.o.f. are entangled), it is known that topological states must have a high degree of nonlocal entanglement as measured by the minimal number of circuit layers needed to prepare them from an unentangled state [11]. The anomalously high entanglement in metallic states similarly implies that any circuit that prepares a metallic state starting from an unentangled state must have many layers.
In terms of calculations, quantum circuits often yield efficient classical algorithms for computing physical properties such as correlation functions. Moreover, they give rise to a local description of the multipartite entanglement structure in terms of multilinear maps. As such, they form an important subclass of tensor network states, a very successful variational class of quantum states that have been shown to be applicable in situations when other methods fail, e.g., because of a fermion sign problem in Monte Carlo methods (for a review, see Refs. [12,13]).
In terms of experiments, quantum circuits provide a precise operational procedure, implementable on a sufficiently versatile quantum simulator or quantum computer, to prepare interesting states. For example, while it might be difficult to directly cool a system to its ground state, a quantum circuit provides an alternative way to directly engineer the ground state.
In this work, we provide quantum circuits that, when acting on a suitable unentangled state, prepare approximations to the metallic ground states of certain one-and twodimensional noninteracting fermion Hamiltonians. This is a nontrivial task in part because these states of matter are highly entangled and violate the area law. We work primarily in the thermodynamic limit of an infinite lattice, but our constructions can also be applied to finite-size systems. The circuits themselves are composed of layers having a hierarchical renormalization group structure, in which the size of the system is doubled (or halved) after each layer. In one dimension, the scheme is a version of the multiscale entanglement renormalization ansatz (MERA) [14], while in two dimensions, it yields an interesting branching network structure [15]. Such renormalizationgroup-inspired quantum circuits have been useful in describing a variety of gapless states in one dimension [16][17][18][19], and inroads have been made in two-dimensional systems [20][21][22][23][24]. They have also been instrumental for recent progress in our understanding of the holographic duality [25][26][27][28]. As in the pioneering work [29], our construction is based on discrete wavelet transforms, although we take a somewhat different approach here.
The key features of our work are the following. Our circuits are hand-crafted-no variational optimization is used-and come with provable error bounds. The essential physics is a certain "half-shift" condition discussed in detail below [30]. Our two-dimensional circuit provides a representation of a state with a Fermi surface, albeit on a special lattice and at a special filling where the Fermi surface has no curvature. The error scaling with circuit size is consistent with the hypothesis that an exact renormalizationgroup circuit can be implemented using a quasilocal circuit with rapidly decaying tails [31]. Together, our results address a long-standing challenge: to rigorously construct tensor networks for gapless states of matter with Fermi surfaces. This paper is organized as follows. We first briefly set the stage for our work and review the basics of noninteracting fermion systems and some wavelet terminology (Sec. II). Then, we discuss the one-dimensional model and construct MERA quantum circuits that approximate its ground state (Sec. III). Next, we generalize to the two-dimensional model, which has a square Fermi surface (Sec. IV). We then explain how to obtain a priori error bounds for our constructions (Sec. V), and we illustrate the effectiveness of our results through numerical experiments (Sec. VI). We conclude in Sec. VII.

II. SETUP
Throughout this paper, we work in the context of noninteracting fermion systems. At the single-particle level, the systems can be described by a Hilbert space spanned by a basis of single-particle states or modes ϕ α . These data depend on the physical setup, but in the cases below, ϕ α can be taken to be a single-particle state in which the particle is localized at site α in some lattice Z d . The many-body description is achieved by second quantization, i.e., passing to creation, a † α , and annihilation, a α , operators that create and destroy a fermion in the single-particle state ϕ α .
The Hamiltonians we study will all consist of fermion bilinears of the form H ¼ α;β can be regarded as a single-particle Hamiltonian acting on the mode space. The eigenstates of such a "quadratic" Hamiltonian are many-body quantum states of fermions that obey Wick's theorem. This, in turn, implies that these states are uniquely specified by the two-point correlation function G α;β ¼ ha † α a β i. In particular, the ground state of H is obtained from the state with no fermions by diagonalizing the matrix h ð1Þ α;β and placing one fermion in each mode with negative single-particle energy. It is therefore possible to carry out much of the analysis of the ground state at the level of single-particle states. In particular, the circuit approximation of the ground state is constructed from single-particle unitaries u ¼ e iz . The corresponding "quadratic" many-body unitary U ¼ e iZ with Z ¼ P α;β a † α z α;β a β then acts as U † a α U ¼ P β u α;β a β . Any state obeying Wick's theorem can always be prepared from an unentangled state (consisting of fermions localized to sites) by acting with such a "quadratic" unitary.
The models we study are translation invariant, so we immediately know how to diagonalize the single-particle Hamiltonian h ð1Þ using the Fourier transform (up to a small eigenvalue problem in the case of having several modes per site). However, the Fourier transform is not a local unitary transformation (it mixes modes arbitrarily far away), so it fails to fully capture the special physics of the ground state and its real-space entanglement renormalization structure. For example, the Fourier transform typically takes an unentangled state to a state with volume law entanglement. Importantly, however, the ground state is invariant under basis transformations within the filled and empty spaces (negative-and positive-energy eigenspaces, respectively), where the former is also known as the Fermi sea. We can therefore prepare the same state by filling modes that are suitable linear combinations of negative-energy eigenstates, chosen to be approximately local in real space. Vice versa, we can approximate the ground state by filling strictly local modes that are approximately supported within the Fermi sea.
Wavelets [32] provide a suitable basis to construct such local modes. As first discussed in Ref. [29], the hierarchical structure of a wavelet transform provides the single-particle version of an entanglement renormalization quantum circuit. In one dimension, wavelet transforms can be specified by a scaling filter h s and a wavelet filter h w . An input signal ψ ∈ l 2 ðZÞ (the space of square summable sequences) is then decomposed by convolution and downsampling into a scaling output ψ s ½n ¼ P m h Ã s ½m − 2nψ½m and a wavelet output Intuitively, the wavelet filter should project onto details of a certain scale, while the scaling filter should project on all features up to this scale. The (discrete) wavelet transform is obtained by iterating this scheme: The scaling output is taken as the input signal for the next iteration. It decomposes the Hilbert space into orthogonal subspaces, each describing details at a certain scale. Its inverse reassembles the input signal from the wavelet outputs at all scales.
If we design the wavelet transform such that it separates negative-energy from positive-energy modes, we obtain a renormalization scheme from the "quadratic" unitary U RG corresponding to one step of the wavelet transform. If the filters have finite length, then U RG is a finite-depth quantum circuit [33], meaning it is composed of a finite number of layers of two-site unitaries. The unitary U RG constitutes one renormalization step: Given the ground state jψi of the Hamiltonian, where, on the right-hand side, jψi is the ground state on the renormalized lattice and jχi is some unentangled state. Crucially, the disentangled sites are interleaved with the renormalized lattice, and each unitary layer is a local transformation. By composing many layers of U RG , we thus obtain a quantum circuit that approximately prepares the ground state. The layout of the circuit is illustrated in Fig. 1. The bottom of the figure corresponds to the state jψi, each layer of red and green blocks constitutes the quantum circuit implementing U RG , the product states j1ij0i on half of the sites make up jχi, and the lines that go up into the next layer correspond to jψi on the other half of the sites, which can be identified with the renormalized lattice. To realize this approach, we still need to design finite-length filters h s , h w such that the wavelet transform separates negative-from positive-energy modes. We now discuss in detail how this can be done systematically and to arbitrarily high fidelity for two fundamental model systems.

III. FERMIONS ON THE DISCRETE LINE
We first consider the fermion nearest-neighbor hopping Hamiltonian on the one-dimensional infinite discrete line, After blocking neighboring sites using the modes b 1;n ¼ ð−1Þ n a 2n and b 2;n ¼ ð−1Þ n a 2nþ1 , corresponding to the even and odd sublattices, respectively, we can write In terms of momentum modes b j ðkÞ ¼ P n b j;n e −ink , the Hamiltonian is FIG. 1. MERA quantum circuit preparing an approximate ground state of the one-dimensional hopping Hamiltonian (1), using the notation of Ref. [29]. As discussed in Sec. III, the tensor network consists of an initial "preprocessing" step followed by repeated layers of "renormalization-group" circuits. Each repeated layer consists of a finite depth circuit followed by projection of half the d.o.f. onto a product state. In more detail, we first apply phase gates (yellow circles) and swapping sublattices in order to let the following gates act either on the even or the odd sublattice. The bulk of the MERA corresponds to two independent wavelet transforms, one for h s (even sublattice) and one for g s (odd sublattice). Each step of the two wavelet transforms gives rise to a layer (red boxes) of two parallel quantum circuits of depth K þ L, consisting of 2 × 2 unitary gates (solid or hatched blue boxes). At the end of each layer, we apply second quantized Hadamard unitaries (green boxes) that couple both sublattices, and we occupy the negative-energy modes. The figure illustrates the case of K ¼ L ¼ 1 and L ¼ 2 layers. As we increase K, L, and the number of layers L, we systematically obtain better and better approximations to the ground state (see Sec. V).
To design a quantum circuit for the ground state, it is convenient to diagonalize the single-particle Hamiltonian into negative-and positive-energy eigenmodes by using the where, importantly, we are free to choose a k-dependent phase. Note that the matrix dðkÞ is discontinuous around k ¼ 0 because of the sign function, but not around k ¼ AEπ, where the discontinuity in the sign is canceled by the discontinuity in the half-shift phase factor (and the result is even smooth). The many-fermion ground state corresponding to the diagonalized single-particle Hamiltonian is disentangled and can be prepared in a completely local fashion by filling the even sublattice, corresponding to the first component, while leaving the odd sublattice empty. We now show that the "quadratic" unitary corresponding to uðkÞ can be well approximated by a finite-depth quantum circuit. The Hadamard h 2 is not k dependent, and thus, its second quantization simply corresponds to a local unitary between neighboring sites of the original nonblocked lattice. Hence, it suffices to focus on the unitary dðkÞ, which is block diagonal between the even and odd sublattices. In view of the quantum circuit/wavelet correspondence discussed in Sec. II, we thus need to design a pair of wavelet transforms, acting on the even and odd sublattices and specified by filters h s , h w and g s , g w , respectively, whose Fourier transforms are related by One can verify that Eq. (5) is fulfilled if the corresponding scaling filters satisfy [36] h s ðkÞ ≈ e iðk=2Þ g s ðkÞ: ð6Þ The phase difference e ik=2 in Fourier space implies that the two scaling filters are related by a half-shift or half-delay in real space. Its appearance is not surprising, given the translation invariance of the original (unblocked) lattice [37]. It is easily seen that the outputs of the inverse wavelet transforms are then, at all levels, related approximately as in Eq. (3), as illustrated in Fig. 2; thus, they can be used to implement dðkÞ [38]. In other words, the same filters can be used throughout, and a scale-invariant circuit will be obtained.
Because of the discontinuity of the half-shift at k ¼ AEπ, a pair of local filters cannot satisfy Eq. (6) exactly. Fortunately, approximate solutions were studied in great detail in the context of filter design in the signal-processing literature [30,37]. Selesnick devised a general algorithm to construct filter pairs, indexed by two integers K, L and having length 2ðK þ LÞ, whose Fourier transforms have exactly equal magnitude and differ by a phase e iδðkÞ [30]. The parameter K determines the usual moment condition used in the wavelet literature, which controls the smoothness of the wavelets and the localization properties of the filters in momentum space. The difference between e iδðkÞ and the ideal half-shift is controlled by the parameter L and decreases quickly in the region around k ¼ 0. While e iδðkÞ is continuous at k ¼ AEπ and therefore necessarily deviates from the half-shift in this region, the support of the scaling filter is, in this same region, suppressed with increased K. This allows us to control the error of the approximation (6) FIG. 2. Fourier spectrum of the outputs of the inverse wavelet transform for levels l ¼ 1, 2, 3 and wavelet parameters K ¼ 4, L ¼ 6. Both filters have equal magnitude in momentum space (shaded regions). The relative phase difference of their Fourier transforms (solid lines) approximates the exact phase difference of the two components of Eq. (3) or, equivalently, of the negativeenergy eigenstates of the single-particle Hamiltonian (black dotted line). by increasing the parameters K and L (see right panel of Fig. 3).
We thus obtain entanglement renormalization quantum circuits by combining the circuits for the "quadratic" unitaries corresponding to the wavelet transforms, constructed using the procedure described in Sec. II, with the Hadamard unitaries and the disentangled ground state of the diagonal Hamiltonian (4). These circuits, illustrated in Fig. 1, are composed of self-similar layers, each of which is a quantum circuit of finite depth K þ L that consists of nearestneighbor 2 × 2-unitary matrices. This corresponds to a bond dimension χ ¼ 2 KþL if the circuit is represented in the standard form of a binary MERA, written in terms of single layers of disentanglers and isometries [14,39]. These quantum circuits allow us to rigorously approximate correlation functions of the ground state of the Hamiltonian (1) as discussed in Sec. V and illustrated numerically in Sec. VI.
From the perspective of the renormalization group, it is natural to consider the coarse-grained or renormalized Hamiltonian. Recall that the original single-particle Hamiltonian is of the form h ð1Þ ðkÞ ¼ eðkÞ( cosðk=2Þσ y − sinðk=2Þσ x ), where eðkÞ ¼ 2 sinðk=2Þ. Because of the downsampling, both the wavelet and scaling outputs couple h ð1Þ ðkÞ and h ð1Þ ðk þ πÞ (i.e., a single layer of a binary MERA is invariant under shifts over two sites). The Hamiltonian can be naturally divided into three terms-corresponding to the scaling modes, the wavelet modes, and the mutual "interaction" between scaling and wavelet modes, respectively, each of which are a free-fermion Hamiltonian. The wavelet Hamiltonian takes the exact single-particle form −ϵ l ðkÞσ z (after the additional local Hadamard transforms). Here, l denotes the level of the wavelet transform, viz. the layer of the MERA, and ϵ l ðkÞ > 0, so that its ground state is a product state in real space, obtained by filling the first mode on every site, in agreement with Eq. (4). If Eq. (6) is satisfied exactly, then the scaling Hamiltonian or renormalized Hamiltonian has the structure e l ðkÞ( cosðk=2Þσ y − sinðk=2Þσ x ), where only the eigenvalues AEe l ðkÞ change with the level l, and not the eigenvectors; in general, this is still approximately true. This is the proclaimed scale invariance, and it provides an alternative way to see that the same pairs of scaling and wavelet filters should be used in every layer. The coarsegrained dispersion relation e l ðkÞ does eventually reach a fixed point (up to a scaling 2 l that accounts for the rescaled lattice spacing), as illustrated in the left panel of Fig. 3. Note that there is also a residual wavelet-scaling interaction term, originating from the overlap between the momentum-space support of the wavelet and scaling filters, so the Hamiltonian is not exactly block diagonal. In particular, the dispersion relation e lþ1 ðkÞ is not simply e l ðk=2Þ (the lower half of the dispersion relation of the preceding level), and lim l→∞ 2 l e l ðkÞ does not simply converge to jkj for all k because of deviations around k ¼ AEπ (see also the left panel in Fig. 3). An exact block diagonalization would require filters with nonoverlapping support, which are therefore nonlocal in real space. While this behavior is more closely approximated with increasing K, the magnitude of the interaction term decays, at most, polynomially in K. However, full block diagonalization of H is too strong of a requirement and would allow the creation of arbitrary eigenstates by replacing the product states with a plane-wave state within a single layer. For the ground state itself, convergence of correlation functions is still exponential in K and L, as discussed in Sec. V.

IV. SQUARE LATTICE AND FERMI SURFACE
We now extend our construction to fermions hopping on the two-dimensional infinite square lattice: We again start by focusing on the single-particle domain and then later transform everything into second-quantized form. The two-dimensional problem we study is special because of the Fermi surface structure: The twodimensional fermion Green function ha † x;y a 0;0 i factorizes into two one-dimensional Green functions, one that depends on x þ y and one that depends on x − y. Thus, as in the one-dimensional case, we decompose the lattice into an even and an odd sublattice, now defined by demanding that the sum of both coordinates is even or odd, respectively; we likewise shift the Brillouin zone by momentum ðπ; πÞ, resulting in new mode operators b 1;x;y ¼ ð−1Þ xþy a xþy;x−y and b 2;x;y ¼ ð−1Þ xþy a xþyþ1;x−y , with corresponding momentum modes b i ðk x ; k y Þ, i ¼ 1, 2. Note that these momenta are now defined with respect to the even or odd sublattice and hence are rotated by 45 degrees with respect to the original lattice. This transformation effectively decouples the x and y directions, as the corresponding one-particle Hamiltonian is now of the form Its eigenvalues are products of the eigenvalues in the one-dimensional case, AE4 sinðk x =2Þ sinðk y =2Þ, with eigenmodes AEe iðk x þk y =2Þ As in the one-dimensional case, the eigenmodes exhibit a phase difference between the two sublattices corresponding to a half-shift in real space, but now the half-shift is in both lattice directions. The positive-and negative-energy eigenmodes are given by ϕ AEsignðk x Þsignðk y Þ ðk x ; k y Þ, respectively, and are thus discontinuous around both k x ¼ 0 and k y ¼ 0, as illustrated in the left panel of Fig. 4.
It is now clear that we can diagonalize the single-particle Hamiltonian with the unitary uðkÞ ¼ dðk x Þdðk y Þh 2 , where d is the block-diagonal unitary (3) and h 2 the Hadamard gate defined previously. We can implement using the tensor products of two one-dimensional wavelet transforms as before-one acting in the x direction and the other in the y direction. More specifically, let us denote by Wψ ¼ ψ s ⊕ ψ w a single step of the one-dimensional wavelet transform with filters h w , h s . Then, which we identify as a single step of the two-dimensional separable wavelet transform. In particular, the waveletwavelet component ψ ww corresponds to the filter h ww ðk x ; k y Þ ¼ h w ðk x Þh w ðk y Þ, and similarly if we use the one-dimensional filters g s , g w instead. Thus, Eq. (5) implies that h ww ðk x ; k y Þ ≈ −signðk x Þsignðk y Þe iðk x þk y =2Þ g ww ðk x ; k y Þ; which is precisely the desired phase relation between the two components of dðk x Þdðk y Þ (see left panel in Fig. 4). To obtain the tensor product of the two wavelet transforms, we now iteratively apply W ⊗ W to the scaling-scaling component ψ ss , as well as W ⊗ I to ψ sw and I ⊗ W to ψ ws [40]. The resulting transform is thus labeled by two levels, l x and l y , corresponding to the number of renormalization steps in each direction. After the second quantization procedure and converting these transformations into a quantum circuit, we obtain an entanglement renormalization quantum circuit of the form sketched in Figs. 5 and 6. This is a particular example of a branching MERA, which generalizes the MERA to allow for logarithmic corrections to the area law [15,41,42] and which was explicitly built with Fermi surfaces in mind. Unlike in the original proposal, our network has three branches instead of two. Indeed, after each layer, we are left with four branches, of which only one can be projected into a product state while the remaining three have to be analyzed further. However, only one of the three branches keeps branching in the higher levels. The other two are further disentangled by ordinary one-dimensional MERAs as in Fig. 1. This ensures that the ground state produced by our network satisfies an appropriate area law of the form SðRÞ ≃ R log 2 R for the entropy of the reduced density matrix of an R × R box. Indeed, let us first recall the estimation of the entanglement entropy in a onedimensional MERA. Each layer is a finite-depth quantum circuit that increases the entanglement entropy of a region by, at most, a constant amount c 1 , so we obtain FIG. 4. Left panel: The relative phase difference between the Fourier transforms of h ww , g ww (smooth surface), which approximates the exact phase difference of the two components of Eq. (8) or, equivalently, of the negative-energy eigenstates of the single-particle Hamiltonian (wireframe mesh). The colored shading of the coordinate plane indicates the momentum-space support of h ww and g ww . It is concentrated around k x ¼ k y ¼ AEπ and vanishes for k x ¼ 0 or k y ¼ 0. S 1D ðRÞ ≤ c 1 þ S 1D ðR=2Þ ≤ … ≤ c 1 log 2 R. For a regular two-dimensional MERA, every layer can increase the entanglement entropy of an R × R box by c 2 R, leading to S 2D ðRÞ ≤ c 2 R þ S 2D ðR=2Þ ≤ … ≤ 2c 2 R. Thus, the entanglement entropy in a regular two-dimensional MERA obeys a strict area law. In contrast, our branching MERA adds, in every layer, the entanglement contribution of a collection of one-dimensional MERAs in the horizontal and vertical directions. The resulting entanglement entropy is bounded by While only an upper bound, this estimate illustrates how a logarithmic violation of the area law can be obtained from the one-dimensional MERAs in each layer.
From the perspective of the renormalization group, the scaling-scaling branch gives rise to a renormalized Hamiltonian whose eigenmode structure is exactly the same as that of the original Hamiltonian, so it can indeed be further processed in a self-similar fashion. The eigenvalues of the renormalized Hamiltonian converge to a fixed point upon successive coarse-graining (see right panel of Fig. 4). The other two branches, resulting from a scaling filter in one direction and a wavelet filter in the other direction, give rise to the coarse-grained Hamiltonians depicted in Fig. 7. The structure of their eigenmodes is FIG. 5. Branching MERA quantum circuit preparing an approximate ground state of the two-dimensional hopping Hamiltonian (7). The bottom layer consists of applying onedimensional wavelet transforms in both the x and y directions of the 45-degree rotated lattice. This gives rise to four outputs, as illustrated in detail in Fig. 6. The wavelet-wavelet output corresponds to approximate eigenstates and hence can be projected out after a subsequent Hadamard transform that couples the even and the odd sublattice, indicated by yellow triangles. The mixed wavelet-scaling outputs are disentangled in one direction and therefore have to be processed further by one-dimensional transformations in the scaling directions, giving rise to two branches that take the form of binary trees as in Fig. 1. The scaling-scaling output contains no disentangled d.o.f. and is thus connected to a copy of the circuit on the renormalized lattice. purely one dimensional. Indeed, for both outputs, one direction is already of wavelet type, so we only have to apply the one-dimensional MERA in the other direction to obtain wavelet outputs at each scale.

V. RIGOROUS ERROR ESTIMATES
In Secs. III and IV, we constructed entanglement renormalization quantum circuits to approximately prepare the ground state of free-fermion Hamiltonians in one and two dimensions, and we gave a heuristic account of the improved quality of our approximations with increased circuit parameters K and L. We now discuss how this intuition can be turned into rigorous a priori error estimates. For simplicity, we only formulate our result in the one-dimensional case, but its statement and proof are completely analogous for two dimensions.
Our theorem is stated in terms of correlation functions of fermion creation and annihilation operators. Given a sequence f ∈ l 2 ðZÞ, we define the corresponding annihilation and creation operators via aðfÞ ¼ P n∈Z f½na n and a † ðfÞ ¼ P n∈Z f½n Ã a † n . We are interested in computing correlation functions of 2N creation and annihilation operators in a many-body state Ψ, Other orderings of operators can be obtained by using the canonical anticommutation relations faðfÞ; aðgÞ † g ¼ hfjgi and faðfÞ; aðgÞg ¼ 0. The number of creation and annihilation operators must be equal to obtain a nonvanishing result since we are interested in states that are invariant (up to an overall phase) under a global Uð1Þ (particle number) transformation of the form a α ↦ e iθ a α . For a pure state of a finite-size system, this invariance would simply imply that the state has a fixed number of particles. Let Dðff i gÞ denote the maximal support of any linear combination of the observables f i (e.g., n for an n-point function). We find that correlation functions of sparse observables are easier to approximate.
Our result is independent of any specific filter construction and only depends on the following parameters. Let h s and g s be two scaling filters of finite length M such that the half-delay condition (6) is approximately satisfied: We also assume that the filters generate corresponding multiresolution analyses with scaling functions bounded in absolute value by some constant B ≥ 1. Then, we have the following a priori error estimate: Theorem 1 Let jΩi denote the exact ground state of the Hamiltonian (1) and jΩ MERA i the many-body state prepared by L layers of the MERA quantum circuit constructed from two scaling filters as above. Then, we have the following error bound for correlation functions: For all f 1 ; …; f 2N with jf i j 2 ≤ 1, where the constant C is given by 2 3=2 ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi Dðff i gÞ p BM. Theorem 1 shows that correlation functions can be approximated to arbitrarily high fidelity for a MERA constructed from suitable scaling filters. As discussed in Sec. III, Selesnick's algorithm gives rise to such filters, parametrized by two integers K and L [30]. Their length is M ¼ K þ L, and we numerically find that B remains bounded, while ε decreases exponentially as we increase K and L (see Fig. 8). Thus, the error bound in Theorem 1 is likewise exponentially small if the number of layers L is sufficiently large. We illustrate this in Sec. VI below, where we numerically approximate the energy density and more general two-point functions. FIG. 7. Coarse-grained Hamiltonians for the wavelet-scaling and scaling-wavelet branches. We show the positive-energy branch (surface) and the relative phase difference between the two components of the positive-energy eigenmodes (color coding of the coordinate plane). Because the eigenmodes are independent of one of the two directions and take the form of the onedimensional problem studied in Sec. III, the ground state can be disentangled by a tensor product of one-dimensional MERAs as in Fig. 1.   FIG. 8. Illustration of the error term ε in Theorem 1. The error ε decreases exponentially for Selesnick's construction as the parameters K and L are increased.
It is instructive to consider a few features of Theorem 1. Suppose f 1 ½n ¼ δ n;x and f 2 ½n ¼ δ n;y . Then, Gðff i gÞ is simply the two-point function Cðx; yÞ ¼ ha † x a y i, which, in the true ground state, decays with jy − xj as a power law, C ∼ jy − xj −1=2 . Yet, Theorem 1 gives a bound that is independent of the separation jy − xj. This might seem puzzling since for a finite depth L, all correlations between operators separated by more than 2 L M vanish. However, at a distance of jy − xj ¼ M2 L , the two-point function is of order C ∼ M −1=2 2 −L=2 , which is consistent with Theorem 1: Dropping the second term in the square root, we still have More generally, the two terms in the square root in Theorem 1 have different physical interpretations. The first is associated with the convergence of the renormalizationgroup transformation, while the second is associated with the goodness of approximation of the phase relation. Indeed, Eq. (10) requires that the phase relation (5) is approximately correct or, when this is not the case for some k, that both h s ðkÞ and g s ðkÞ are small in magnitude (cf. Sec. III).
Our proof of Theorem 1 makes this intuition precise. We show that Eq. (10) guarantees that the single-particle modes prepared by the MERA are approximate eigenmodes, and the boundedness of the scaling function ensures that the truncation error decreases exponentially with the number of layers of the tensor network. Together, this implies that the two-point correlation functions of the states jΩi and jΩ MERA i are approximately equal. We then use a robust version of Wick's theorem [43] to show that higher correlation functions can likewise be approximated up to small errors. We refer to the Appendix for a rigorous mathematical proof.
It is remarkable that the error converges as L → ∞: Even though correlation functions now depend on an infinite number of "nonideal" (finite ε) layers, the total error is bounded. This is a consequence of the hierarchical renormalization-group structure of the network combined with the boundedness of the scaling functions.
Note that Theorem 1 does not provide an error estimate on the fidelity between the true ground state and the MERA state for an infinite system. Indeed, these two states are expected to necessarily be orthogonal in the thermodynamic limit since any finite error per unit volume will result in zero overlap as the system size is taken to infinity. Nevertheless, Theorem 1 proves that our construction can yield correlation functions that approximate those of the true ground state to arbitrary accuracy. Therefore, all intensive (not scaling with system size) physical properties that can be inferred from these are likewise well reproduced. Our results can thus be seen as another instance where we can rigorously construct tensor network states for critical systems or for quantum field theories if we focus on correlation functions, a point first raised in Refs. [44,45].
On a finite ring of size V, the one-dimensional model has an energy gap Δ ∼ 1=V. In such a situation, the infinitesystem L → ∞ circuit must be modified to fit into the finite-size system. We expect that there exists an analogue of Theorem 1 that guarantees correlation functions are well approximated for sufficiently small ε. Moreover, in this finite-size setting and with sufficiently small ε, the state jΩ MERA i can have high overlap with jΩi. Indeed, if P Ω ¼ jΩihΩj is the ground-state projector and E Ω is the groundstate energy, then we have Δð1 − P Ω Þ ≤ ðH − E Ω Þ and hence Thus, if the energies of jΩi and jΩ MERA i are within 1=polyðVÞ of each other, then the overlap jhΩjΩ MERA ij 2 is within 1=polyðVÞ of one. If, as suggested by our numerics, the error ε achieved by Selesnick's wavelets, say, decreases exponentially with minðK; LÞ, then one would have high overlap between a MERA state and the true ground state using a bond dimension scaling only polynomially in V.

VI. NUMERICAL RESULTS
Our construction can be used to effectively calculate physical properties in real space [46]. For example, consider the energy density of the approximate ground state. Its value for the MERA quantum circuit for the infinite onedimensional discrete line, truncated at depth L, is given by P L l¼1 2 −ðlþ1Þ e ðlÞ , with e ðlÞ ¼ −2 Re P n ϕ ðlÞ ½nϕ Ã ðlÞ ½n þ 1 the single-particle energy of a mode ϕ ðlÞ obtained from the lth layer. The scaling factor comes from the fact that at the lth level of the MERA, the density of d.o.f. is 2 −l , half of which are filled. This can be easily evaluated numerically, and it displays convergence to the true value −2=π, as illustrated in the left panel of Fig. 9. The numerical results are consistent with a power-law convergence in the effective bond dimension χ ¼ 2 KþL , in agreement with our discussion below Theorem 1. We find that our analytical construction systematically improves over the one from Ref. [29], but its energy-density estimate is outperformed by the variationally optimized non-Gaussian MERA from Ref. [47].
A similar procedure works for the two-dimensional square lattice. The energy density is now given by P L x l x ¼1 P L y l y ¼1 2 −ðl x þl y þ1Þ e ðl x ;l y Þ , where e ðl x ;l y Þ denotes the single-particle energy of a mode obtained from levels l x , l y of the quantum circuit, which, we recall, denote the number of renormalization steps in the x and y directions, respectively. In other words, minðl x ; l y Þ is the level at which we switch to a one-dimensional branch (cf. Fig. 5). It is useful to note that, since the two wavelet transforms involved are separable, the modes obtained on each sublattice are tensor products of one-dimensional modes, coupled only by the final Hadamard transforms. This allows us to carry out all computations in the onedimensional setting. Our numerical results are shown in the right panel of Fig. 9, and they again show power-law convergence in the effective bond dimension to the true value −8=π 2 .
As a last example, we consider a general two-point function Cðx; yÞ ¼ ha † x a y i. While the true ground state is translation invariant, Cðx; yÞ ≠ Cðy − xÞ for the approximate ground state prepared by the quantum circuit since the latter is not perfectly invariant under arbitrary lattice translations. For simplicity, we only discuss the one-dimensional case. As above, let ϕ ðlÞ denote a single-particle mode obtained from the lth level of the MERA quantum circuit. Then, we have Cðx; yÞ ¼ P L l¼1 P z∈Z ϕ ðlÞ ½x − 2 lþ1 zϕ Ã ðlÞ ½y − 2 lþ1 z. The inner sum corresponds to the different modes obtained from the lth level, obtained as translations of ϕ ðlÞ ; we note that only finitely many translations yield a nonzero summand since the ϕ ðlÞ are finitely supported. The result is shown in Fig. 10. Again, we find convergence to the exact solution Cðy − xÞ ¼ sinðπðy − xÞ=2Þ=ðπðy − xÞÞ. In particular, the two-point function becomes more and more translation invariant with increased K, L.

VII. CONCLUSIONS
In this work, we showed how wavelet theory can be used to rigorously construct quantum circuits that approximate metallic states of noninteracting fermions. Working directly in the thermodynamic limit, we showed that arbitrary correlation functions of fermion creation and annihilation operators can be approximated to high accuracy for an appropriate choice of circuit parameters. In a finite-size system, we argued, based on our numerics, that a tensor network with bond dimension scaling only polynomially with system size can achieve unit overlap with the true ground state in the large-system-size limit. Although such a bond dimension is high from the point of view of numerical calculations using a classical computer, from an information-theoretic point of view, it represents an astounding compression of the quantum state. At no point did we use a variational optimization to determine the circuit parameters, and the circuits we construct have a plain physical meaning. The essential physics arose from the structure of negative-and positive-energy eigenspaces and was encapsulated in a half-shift delay between pairs of wavelet filters. The design of such pairs of wavelets had already been carried out in the signal-processing community.
The constructions reported here are closely related to a forthcoming work by three of the authors; it uses wavelet technology to approximate correlation functions in a continuum quantum field theory, namely, the free Dirac field in 1 þ 1 dimensions. As in the case of the thermodynamic limit of the lattice system, the correct notion of approximation turns out to be the approximation of correlation functions instead of the approximation of states. Using the free Dirac field construction, it is also possible to construct MERA circuits that approximate the correlation functions of interacting Wess-Zumino-Witten field theories in 1 þ 1 dimensions.
There are many immediate directions for further development. It is of considerable interest to adapt existing wavelets or design new wavelets that can capture curved Fermi surfaces; then, we would truly be able to describe a general class of metallic states in two or more dimensions. This would, for example, enable us to address different chemical potentials in the square lattice model. It is also interesting to adapt our wavelet approach to describe Dirac . For comparison, we also display the relative error for the analytical construction from Ref. [29] and for the numerically optimized non-Gaussian MERA from Ref. [47]. points in two or more dimensions; the basic approach used here is clearly sound, but there is an interesting wavelet design problem to capture the physics of the filled Dirac sea. Another very interesting direction is interacting fermions. For example, similar in spirit to Slater-Jastrow wave functions, our noninteracting wavelet MERA construction might be used as the starting point for a variational class of wave functions for interacting metals. has independent unpublished data, arriving at a similar 2D network structure, and assisted in preparing the manuscript.

APPENDIX: PROOF OF THEOREM 1
We start by describing the setup in precise mathematical terms. Any pure gauge-invariant generalized free state Ψ can be described by an operator ψ ≥ 0, known as the symbol, such that the correlation functions are given by For pure states, ψ is a projection that we can think of as the projection onto the Fermi sea. To connect with physics notation, note that "gauge-invariant" effectively means that the density matrix Ψ is invariant under a global Uð1Þ (particle number) transformation of the form a α ↦ e iθ a α . Both the true ground state and the state prepared by the MERA are pure gauge-invariant generalized free states. Following the discussion in Sec. III, their symbols can be described as follows. We denote by v∶l 2 ðZÞ ⊗ C 2 → l 2 ðZÞ the unitary corresponding to the transformation ðb 1 ; b 2 Þ ↦ a and by mðθ w Þ∶l 2 ðZÞ → l 2 ðZÞ the Fourier multiplier by θ w ðkÞ ¼ −isignðkÞe iðk=2Þ , so the operator (3) can be written as d¼½ I 0 0 mðθ w Þ . Recall that u ¼ dðI ⊗ h 2 Þ, with h 2 the Hadamard gate. Then, the symbol of the true ground state jΩi is given by where jþi ¼ ðj0i þ j1iÞ= ffiffi ffi 2 p . Next, recall that we are given two pairs of filters, h s , h w and g s , g w . We denote the corresponding wavelet transforms, truncated at level L, by where the first L coordinates of the output correspond to the wavelet outputs and the last one to the remaining scaling output (see, e.g., Ref. [32] for an introduction to wavelet theory). Let us denote by p ðLÞ w , p ðLÞ s the projection onto the wavelet outputs and the scaling output, respectively. Then, the many-body ground state jΩi MERA prepared by the MERA quantum circuit has the symbol Let F⊆l 2 ðZÞ denote a subspace (which we will later take to be the span of the observables f i ). Let DðFÞ ≔ sup f∈F jfx ∈ Z∶fðxÞ ≠ 0gj denote the maximal support of any sequence in F. We denote by p F the orthogonal projector onto F and abbreviate as ωj F ≔ p F ωp F . As usual, we write ∥ − ∥ p for p-norms, ∥ − ∥ ∞ for supremum norms, and ∥ − ∥ for operator norms. We use mðθÞ∶l 2 ðZÞ → l 2 ðZÞ more generally for the Fourier multiplier by some periodic function θðkÞ.
We first prove that the truncation of the MERA only incurs an error that is exponentially small in L.
Lemma 1. Let h s ∈ l 2 ðZÞ be a scaling filter of length M such that the associated scaling function ϕ h ∈ L 2 ðRÞ is bounded. Let f ∈ l 1 ðZÞ. Then, For this, let us denote by Ωj F and Ω MERA j F the mixed gauge-invariant generalized free states with symbols ωj F and ω MERA j F , respectively, which capture all n-point functions with observables from F. It is clear from Eq. (A1) that where ∥ − ∥ 1 denotes the trace norm. We now use a result by Powers and Størmer to bound the trace-norm distance between generalized free states in terms of the trace-norm distance of their symbol. Specifically, we use Lemmas 4.1 and 4.6 of Ref. [43] to obtain the first inequality in Nδ p (as long as the right-hand side is smaller than 1=6); for the second inequality, we use dim F ≤ 2N, and the last one is Eq. (A8). We have thus established Eq. (A9) and, thereby, the theorem. □