Bosonic entanglement renormalization circuits from wavelet theory

Entanglement renormalization is a unitary real-space renormalization scheme. The corresponding quantum circuits or tensor networks are known as MERA, and they are particularly well-suited to describing quantum systems at criticality. In this work we show how to construct Gaussian bosonic quantum circuits that implement entanglement renormalization for ground states of arbitrary free bosonic chains. The construction is based on wavelet theory, and the dispersion relation of the Hamiltonian is translated into a filter design problem. We give a general algorithm that approximately solves this design problem and provide an approximation theory that relates the properties of the filters to the accuracy of the corresponding quantum circuits. Finally, we explain how the continuum limit (a free bosonic quantum field) emerges naturally from the wavelet construction.


Introduction
An important task in the study of quantum many-body systems is finding useful parameterizations of physically relevant quantum states. One successful approach is to consider so-called tensor network states, which are defined by contractions of local tensors according to a network or graph structure. This gives a natural way to prescribe the entanglement structure of the state, while retaining the ability to describe interesting states such as low energy states of local Hamiltonians. See [1][2][3][4] for reviews of tensor network states. Tensor networks are particularly useful to implement real-space renormalization methods for strongly interacting quantum many-body systems. In one spatial dimension, prominent examples are the density matrix renormalization group [5], with the associated tensor network class of matrix product states (MPS) [6] and entanglement renormalization [7], with the corresponding multiscale entanglement renormalization ansatz (MERA) states [7,8]. Entanglement renormalization implements a real-space renormalization by a local unitary transformation, decomposing a state into a product state and the renormalized state. By applying many such layers one can build a highly entangled state from product states. Scale-invariant MERA states are a good variational class for approximating ground states of critical quantum chains and one can extract the conformal data of the continuum limit conformal field theory of the system from the entanglement renormalization superoperator [9]. If the entanglement renormalization unitaries are implemented by low-depth local quantum circuits we will call this an entanglement renormalization circuit -see Fig. 1 for an illustration. This class of states can be prepared efficiently on a quantum computer, which makes them a promising ansatz class for variational optimization on a quantum computer. This latter perspective was introduced in [10], where the corresponding class was called DMERA. The contraction cost using known classical contraction algorithms of such DMERA states increases exponentially with the depth of the quantum circuit, compared to which the contraction of these states is exponentially faster on a quantum computer. Another appealing property of entanglement renormalization circuits is that they are robust to small errors, which makes them interesting candidates for noisy intermediate-scale quantum (NISQ) devices [10,11]. Entanglement renormalization circuits are appealing as their depth is logarithmic in the system size, and the circuit depth of a single layer typically scales polylogarithmically in the desired error, and they apply to gapless systems. See [12,13] for some other applications of tensor networks for quantum computing.
Unfortunately our analytic understanding of MERA in general and DMERA in particular is still limited (as compared to for instance MPS). One direction in which progress to analytic understanding has been made is in connection to wavelets. Wavelet Figure 1: The structure of an entanglement renormalization circuit. Each layer is a constant depth quantum circuit that is supposed to implement a real-space renormalization. Every layer takes as input the output of the previous layer and a product state, resulting in an entangled quantum state at the bottom. Layers further up in the figure correspond to structure at larger scales.
a signal as a linear combination of localized wave packets or 'wavelets' at different scales (as compared to the Fourier transform, which uses plane waves). This can be implemented iteratively: In each step the signal is decomposed into a high-frequency component (the 'details' of the signal) and a low-frequency component (the 'large scale structure' of the signal). The wavelet transform then proceeds iteratively on the low-frequency component of the signal. Wavelet theory has many and wide-ranging applications, from practical signal processing applications such as image compression [14] to mathematical analysis [15]. The procedure of the wavelet transform is very similar to real-space renormalization, and its original development was partially motivated by applications in real-space renormalization.
Recently it has been observed [16][17][18] that any finite wavelet transform can be written as a classical linear circuit whose fermionic second quantization gives rise to a free fermionic (D)MERA. Moreover, the continuum limit can be precisely related to the corresponding wavelet functions [19]. In [16] it was suggested that a similar result could also be true for free bosonic systems. In order to formulate what this entails, we will work with bosonic quantum circuits. This means that we have a set of bosonic modes, and a bosonic quantum circuit will be a sequence of operations acting locally on these bosonic modes. We restrict to the subclass of Gaussian or linear optics circuits, meaning that each local operation is implemented by time evolution with a quadratic Hamiltonian. This is an efficiently simulable subclass of all bosonic quantum circuits (upon adding non-Gaussian bosonic quantum gates, however, bosonic quantum circuits are able of universal quantum computation [20]). In contrast to more usual notions of quantum circuits and tensor networks, the Hilbert spaces are infinite dimensional. In particular, the usual definition of a tensor network with a finite bond dimension has no immediate analogue. However, finite-depth quantum circuits such as entanglement renormalization circuits of the form of Fig. 1 are still meaningful even in this infinite-dimensional bosonic setup. The notion of Gaussian bosonic entanglement renormalization has been introduced and studied in [21], in which an extensive explanation of the formalism can be found.

Main results
In this work we show that one can indeed construct a Gaussian bosonic entanglement renormalization scheme for bosonic quadratic one-dimensional Hamiltonians, using the second quantization of biorthogonal wavelet filters or perfect reconstruction filters. This extends the wavelet-MERA correspondence substantially. The resulting entanglement renormalization takes the form of a short-depth Gaussian bosonic circuit, providing evidence for the relevance of entan-glement renormalization circuits for preparing ground states of (near) critical quantum systems. Moreover we can relate, similar as in the fermionic case [17], properties of the biorthogonal wavelet transform to the resulting MERA state, and we prove a rigorous approximation theorem for the correlation functions of the MERA state. Interestingly, our formalism is not restricted to the scale-invariant case, but can be used to construct entanglement renormalization circuits for arbitrary translation invariant quadratic bosonic Hamiltonians. Given such a Hamiltonian, we explain how a corresponding (approximate) entanglement renormalization circuit can be found by solving a filter design problem. We also give a general method for constructing such filters, similar to the construction of the Daubechies wavelets. This is in contrast to the fermionic case, where the only known constructions are for massless (critical) fermions [16,17]. Finally, the continuum limit of the discrete system is directly related to the biorthogonal scaling and wavelet functions corresponding to the filters. For the free massless boson our construction reproduces various scaling dimensions exactly. If the system is not scale-invariant, we explain how one can still define versions of the wavelet and scaling functions which are not scale-invariant.
A natural application of a quantum computer based on bosonic variables [20] is to simulate bosonic quantum field theories [22], and wavelets are a very efficient choice of basis to discretize a quantum field theory for this purpose [23]. We explain that for any free 1+1-dimensional bosonic field theory, one can use suitably chosen biorthogonal wavelets to discretize the theory and use the corresponding wavelet decomposition to prepare its (approximate) ground state using the bosonic Gaussian entanglement renormalization circuit. The idea to use wavelets to discretize a field theory is quite natural, see for instance [23][24][25] for some recent discussions of discretizing bosonic field theories using wavelets. Our approach however fundamentally differs from these works in that we use biorthogonal wavelets (as is natural in the bosonic setting), which moreover are specifically designed to target the Hamiltonian of the field theory (rather than using off-the-shelf wavelets such as the Daubechies wavelets). We hope that our investigations can provide a potential starting point for the efficient simulation of interacting quantum field theories on quantum computers.

Organization of the paper
In Section 2 we give a brief review of biorthogonal filters and wavelet theory. In Section 3 we briefly review the formalism of quadratic bosonic Hamiltonians and Gaussian unitaries. We then explain the relation between entanglement renormalization and biorthogonal wavelet filters. In particular, in Section 3.1 we derive a relation the filters have to satisfy to disentangle the ground state of a given Hamiltonian. In Section 4 we explain how this gives rise to a circuit, and we state Section 4 which proves bounds on the accuracy of the approximation. Finally, in Section 5 we introduce continuous wavelet functions, and show that this gives a natural interpretation of entanglement renormalization in the corresponding quantum field theory. In the appendices we provide a review of the fermionic MERA/wavelet correspondence in Appendix A, an algorithm for constructing appropriate biorthogonal filters in Appendix B, an explanation of how to construct Gaussian circuits from filters in Appendix C and a precise statement and proof of Section 4 in Appendix D. Supporting code used to generate the numerical results in this work can be found at [26].

Perfect reconstruction and biorthogonal filters
Perfect reconstruction filters, or biorthogonal wavelet filters, are filters that decompose a signal into a high-frequency part and a low frequency part. This is reminiscent of the disentangling procedure of entanglement renormalization, and in this work we explain the precise connection. We first give a brief account of the theory of biorthogonal wavelet filters, see [14] for an W g W g low high low high Figure 2: Iterating the wavelet decomposition W g resolves a signal into scales. Illustration for the Haar wavelet filters [19].
introduction. We consider a pair of real-valued sequences g s , h s ∈ 2 ( ), called scaling or low-pass filters. Often they will be finite impulse response (FIR) filters, meaning that they have finite support. We demand that these filters satisfy the perfect reconstruction condition on their Fourier transforms and define corresponding wavelet or high-pass filters by We similarly define W h using the filters h s and h w in place of g s and g w , respectively. By applying W g again to f low , the original signal is recursively resolved into scales, see Fig. 2. It follows from Eq.
(1) that f can be reconstructed from its decomposition W g f by applying the transposed The roles of g and h can be exchanged in this procedure.

Entanglement renormalization and filter design
We will consider translation invariant chains of harmonic oscillators (q n , p n ), with a Hamiltonian of the form where V nm = V n−m defines a positive definite symmetric matrix. The ground state of such a quadratic Hamiltonian is a Gaussian state, determined by the dispersion relation ω(k) of the Hamiltonian. We study Gaussian circuits that map an unentangled state to the entangled ground state of a translation invariant Hamiltonian (or conversely, disentangle the ground state to an unentangled state). We consider Gaussian maps defined byq n = m A nm q m ,p n = m B nm p m . This preserves the canonical commutation relations if and only if the matrices A and B are such that B = (A T ) −1 . By a Gaussian circuit we will hence understand a sequence of Gaussian maps, each of which maps modes (q n , p n ) to a linear combination of itself and its direct neighbours. For details about quadratic bosonic Hamiltonians and Gaussian states see, for instance, [20,[27][28][29].
In Fourier space one simply maps a product state to the state with dispersion relation ω(k) by appropriately 'squeezing' each Fourier mode. However, this is a very non-local operation, whereas we are interested in a procedure that is local in real space. For more discussion of this point and variational algorithms to find Gaussian entanglement renormalization maps, see [21].
Since W −1 g = W T h the map W = W g ⊕W h defines a Gaussian map for any pair of biorthogonal wavelet filters (g, h). This has the structure of a layer of entanglement renormalization, filtering out the high frequency modes. However, we need to choose the filters g and h such that W actually disentangles the state, and the wavelet output is unentangled. If we normalize the dispersion relation such that ω(π) = 1, then the condition for the wavelet output to be disentangled is that the Fourier transforms of the filters satisfy Intuitively, what happens is that W separates the bosonic modes in high frequency and low frequency modes, and Eq. (4) makes sure that the high frequency modes are not entangled to the low frequency modes in the ground state. We derive this condition below in Section 3.1. The scaling (low frequency) modes are again mapped to a Gaussian state, possibly with a different dispersion relation. We can now recursively apply the same construction to the scaling output, as in Fig. 1, now with the renormalized dispersion relation, given by precomposed with a 'squeezing' normalization layer to ensure the normalization ω(π) = 1 before we apply the wavelet decomposition. We will later see this procedure can be decomposed as a circuit, see Fig. 3. While we motivated the procedure from the perspective of disentangling a given entangled state, the resulting circuit can also be used in the opposite direction, to prepare the ground state by applying the circuit to a product state, thus realizing the state as a bosonic MERA state. A paradigmatic example is the harmonic chain, which has dispersion relation ω(k) = m 2 + sin 2 k 2 . In particular, the massless harmonic chain is gapless and has dispersion relation ω(k) = |sin( k 2 )|. For the latter, Eq. (5) amounts to ω(k) → |sin( k 4 ) cos( k 4 )| = 1 2 sin( k 2 ), so the dispersion relation is invariant under the renormalization step if we include the subsequent normalization. Hence the state on the scaling output of the entanglement renormalization will be the same after any number of layers. This implies that we can keep iterating the same entanglement renormalization layer with identical filters at each layer, giving a scale-invariant bosonic entanglement renormalization procedure for the massless harmonic chain.
In the massive case, the mass renormalizes as This is a relevant perturbation to the massless chain [21], and with increasing number of layers the dispersion relation becomes flat; correspondingly we can let the filters at the deeper layers approach orthogonal wavelet filters.

Derivation of filter condition
We will now derive Eq. (4). The ground state of the Hamiltonian in Eq. (3) is completely determined by its covariance matrix γ = γ q ⊕ γ p , whose Fourier transform is given by , The covariance matrix of an unentangled (uncorrelated) product state is 1 2 1. Recall that any symplectic linear map S on the set of modes (q n , p n ) defines a unitary map which maps Gaussian states to Gaussian states. Under a map of the form A⊕ A T −1 the covariance matrix transforms as We first normalize such that ω(π) = 1, which can be implemented by the symplectic (squeezing) map ( ω(π)1) ⊕ (1/ ω(π)1). Suppose we have filters (g, h) satisfying Eq. (4), then W = W g ⊕ W h disentangles the ground state. To see that this is indeed true, we compute the result of applying the wavelet decomposition map to the ground state covariance matrix γ = γ q ⊕ γ p given in terms of the dispersion relation by Eq. (8). For this, we remark that from where ω (1) is the renormalized dispersion relation on the scaling output defined in Eq. (5) in the main text. This shows that Similarly, it holds that We thus see that W has unentangled the high-frequency modes to a product state, and the low frequency modes are renormalized to have a new dispersion relation ω (1) given by Eq. (5). The full entanglement renormalization circuit consists of repeated applications of such layers. To introduce some notation, we let ω (l) be the dispersion relation after l layers of to normalize the dispersion relation. This figure can be interpreted both as a linear circuit implementing a symplectic transformation, and as its second quantization, which is a bosonic Gaussian circuit. renormalization, recursively defined by (cf. Eq. (5), note that we first normalize the dispersion relation by a factor ω (l) (π)) The normalization by ω (l) (π) could also be absorbed in the filters, but we would like the filters to be such that g(0) = h(0) = 2, as is standard in the signal processing literature and convenient for the analysis. Then at the l-th layer we need filters g (l) , h (l) satisfying Eq. (4)), and we let Finally, we define the L-layer renormalization map as Then, R (L) maps the state with dispersion relation ω to a product state with covariance matrix 1 2 1 on the L high frequency levels, and a state with dispersion relation ω (L) on the remaining low frequency level.

Entanglement renormalization circuits
If g and h are FIR filters of size 2M , we show in Appendix C that W gives rise to a Gaussian circuit of depth M that maps the low-frequency modes to the odd sublattice and the highfrequency modes to the even sublattice as shown in Fig. 3. This is exactly the structure of an entanglement renormalization circuit. The converse to this construction is also true: any Gaussian entanglement renormalization circuit as described above arises in this way. When using a finite depth circuit, we may not be able to satisfy the relation in Eq. (4) exactly if ω(k) is not a ratio of trigonometric polynomials. In particular, this is the case for the harmonic chain. In this case we can still hope to approximate the dispersion relation, and correspondingly prepare a state that is close the true ground state. This raises two interesting questions. Firstly, the existence of filters that approximately satisfy Eq. (4) is not clear. In Appendix B we describe an explicit procedure for constructing such filters. Secondly, one can wonder whether a good approximation of the dispersion relation at the level of a single layer will indeed give rise to a good approximation of the ground state. Fortunately, the structure of entanglement renormalization is remarkably robust to small errors [10,30], and in [17] a robustness result for wavelet based fermionic entanglement renormalization was proven. The bosonic setting is somewhat different, as the Hilbert spaces are infinite dimensional. We will now discuss that when the family of filters has a well-defined 'continuum limit', we can nevertheless prove a rigorous approximation theorem.
We would like to bound the approximation error when using L layers of entanglement renormalization. Suppose we are given a family of filter pairs (g (l) , h (l) ) for l = 1, . . . , L, where the l-th pair represents the l-th layer such that so they approximately reproduce the dispersion relation at each layer (up to normalization). Moreover, we need these families of filters to give rise to a 'stable' wavelet decomposition, in the sense that many iterations of the decomposition maps yield a uniformly bounded map. This is a standard assumption in wavelet theory. If we are only interested in an approximation, and the theory flows to either a critical theory or a trivial theory we only need a small number of 'transition layers' and can pick fixed filters (g (l) , h (l) ) = (g, h) for large l. In Appendix D, we prove a general approximation theorem in this setting, which applies to an arbitrary quadratic Hamiltonian. We measure the error in the two-point functions 〈p i p j 〉 and 〈q i q j 〉 (or covariance matrix). In the particular case of the harmonic chain in Eq. (6), our result specializes as follows.
Theorem (Informal). For the harmonic chain with mass m, the approximation error using the MERA state resulting from L layers of entanglement renormalization is bounded by the latter assuming m > 0. In the massless case, the latter bound is replaced by In the massless case, there is an IR divergence and 〈q i q j 〉 is only defined up to a constant, so we define 〈q i q j 〉 by subtracting the divergence; see Eq. (47) in Appendix D for details.
The intuition behind the proof is that the contribution of the L-th layer to the correlation function is bounded by O(2 − L 2 ), so we need O(log 1 δ ) layers to get within error δ (even with perfect filters), while each layer contributes a factor of to the error in the filter relation. Balancing these two contributions yields the desired bound. In Appendix B we provide a construction of filters g, h satisfying Eq. (4) for the massless harmonic chain. This construction depends on two parameters K and L, where K controls the number of vanishing moments of the filters and L controls the accuracy of the approximation of the dispersion relation. This corresponds to a circuit depth of M = K = 2L for a single layer. In Fig. 4 we illustrate the approximation result by numerically computing correlation functions of the massless harmonic chain using these filters [26].
If we denote by M the circuit depth of a single layer, then we find numerically that is exponentially small as a function of M , whereas the other wavelet-dependent parameters we have suppressed above only grow polynomially. Hence, the total required depth of a single layer of entanglement renormalization for a desired error is polylogarithmic in 1 . This shows that our entanglement renormalization circuits prepare the ground state very efficiently: a circuit of depth O(polylog( 1 δ )) achieves an accuracy δ on the correlation functions.

The continuum limit
The discrete wavelet transform has a natural continuum limit, in terms of continuous scaling and wavelet functions. This gives a way to interpret the continuous limit of the entanglement renormalization maps. In the following, we demonstrate that the scaling functions are a natural UV cut-off that is compatible with the entanglement renormalization circuits, and in the critical case we find that we can reproduce certain conformal data exactly from a single layer of renormalization. We consider the free boson, described by bosonic fields φ(x), π(x) and Hamiltonian We are particularly interested in the massless case, which gives rise to a conformal field theory.

Scaling and wavelet functions
The continuum limit of the discrete biorthogonal wavelet transform is determined by the scaling functions. Given biorthogonal wavelet filters g, h the associated scaling functions are defined in Fourier space for a = g, h byφ and the associated wavelet functions bŷ Both have compact support if the filters are finite, an example is shown in Fig. 5 . Moreover, we can define rescaled and shifted versions It then follows that the sets {ψ where we find the coefficients w[l, n] ands[n] precisely by applying the discrete wavelet transformation W h to the signal s. Moreover, if we let then f h l and f g l converge in norm to f as l goes to minus infinity. See [14] for an introduction to scaling and wavelet functions and their properties.

Entanglement renormalization for the massless boson
We now suppose that the biorthogonal wavelet filters g, h are related by the dispersion relation of the massless harmonic chain that was discussed in Section 3, that is, We claim that, in this case, the scaling functions defined in Eq. (12) are related aŝ To verify this claim, we note that as a consequence of Eq. (17) and the relation in Eq. (2) we have h s (k) = |cos( k 2 )|g s (k). Next, from Eq.
This expression implies that γ(k) has to satisfy γ(k) = |cos( k 4 )|γ( k 2 ), and we can easily verify that γ(k) = 2|sin( k 2 )| |k| , which has the right normalization γ(0) = 1. This proves Eq. (18), which in turn, using Eqs. (13) and (17), also implies that Equation (19) shows that wavelet functions are related precisely by the linear dispersion relation of the massless free boson. We consider correlation functions of smeared fields precisely by applying the discrete wavelet transform, and the factor of 2 −l derives from our normalization of the dispersion relation (the 'squeezing layer' in Fig. 3). In other words, the correlation functions will be given precisely by applying the entanglement renormalization circuit to the operators n s[n]q n and ns [n]p n . For general functions f , we may approximate them with scaling functions as in Eq. (15) and thus map (the inner product here is again the L 2 ( ) inner product). The scale l corresponds to a choice of UV cut-off. If l is sufficiently small, then by Eq. (15) on can then compute correlation functions using the (discrete) entanglement renormalization circuit to good approximation. This approach is completely analogous to the fermionic construction described in detail in [19].
It yields a natural way to interpret the continuous limit of bosonic entanglement renormalization as quantum field theory.
As an application, we can consider the entanglement renormalization superoperator Φ, which coarse-grains operators by conjugating with a single layer of the renormalization circuit. For critical lattice models, Φ has been proposed to approximately encode the conformal data of the continuum limit of the theory [9]. For instance, for a primary field in the conformal field theory with scaling dimension λ, there should be a local operator O such that Φ(O) ≈ 2 −λ O. We will now verify that the entanglement renormalization superoperator reproduces exactly the scaling dimensions of the φ and π fields in the massless case, as well as the scaling dimension of a number of descendants (equal to the number of vanishing moments of the wavelet filters), similar as for the fermionic wavelet MERA [16,19]. This is seen by considering the operators O φ (x) = n φ h (x − n)q n and O π (x) = n φ g (x − n)p n for any x ∈ , which are the discretizations of the operators φ(x) and π(x). It can be easily seen that the entanglement renormalization superoperator maps these operators where we use that for the scaling function Eq. (12) it holds that [14] Similarly one finds O π (x) → 1 2 O π ( x 2 ). This corresponds, as expected, to scaling dimensions 0 and 1. If the scaling function is differentiable, we see that by differentiating Eq. (21) we get that 1 , which leads similarly to a descendent field O φ (1) = ∂ x φ h (x − n)q n with the right scaling dimension. It turns out that if φ h has K vanishing moments (or equivalently, a factor (1 + e ik ) K in the scaling filter h s [14]), then there exists a vector φ h,l [n] with l = 1, . . . , K and n taking integer values on the support of the wavelet such that even if φ h is not l time differentiable (note that φ h,l is only defined at integer values), see Theorem 7.1 in [32], and similarly for φ g,l . This shows that computing the eigenvalues of the entanglement renormalization superoperator Φ will yield the eigenvalues of K descendants of the φ and π fields. At this point we observe that a wavelet filter leading to K vanishing moments must have support at least 2K, so one needs (as expected) a larger circuit depth to capture more descendent scaling dimensions. For example, in our explicit constructions in Appendix B the filter size is 2K + 4L where L controls the accuracy of the approximation of the dispersion relation.

The massive bosonic field
The free massive boson with mass m can be approached similarly. In that case, we suppose we have two families of filters g (l) and h (l) , now with l ∈ and such that (m (l) ) 2 + 1 g (l) where m (0) = m and m (l) is the mass after l layers of renormalization, as defined by Eq. (7). If these filters are chosen in a way that they converge to a fixed orthonormal filter as l goes to infinity, and to a fixed pair of biorthogonal filters as in the massless case for l to −∞, it makes sense to define a new type of scaling and wavelet functions which are different at each level l as a generalization of of the scaling and wavelet functions: for a = g, h as a generalization of of the scaling and wavelet functions defined in Eq. (12) and Eq. (13). Again, the wavelet functions ψ a l,n (x) = 2 − l 2 ψ l (2 −l x − n) for a = g, h form a dual basis (provided they exist). The behaviour for l → ±∞ is consistent with the fact that the mass term is a relevant perturbation of the conformal field theory and the theory flows from a critical massless boson to a trivial theory. As before, we can now discretize the theory using the scaling functions at some given scale and use the discrete circuit to compute correlation functions.

Other perspectives
The idea that wavelet theory should be a natural tool to discretize a field theory in order to perform renormalization has a long history [33]. As mentioned in the introduction, our approach differs from other works such as [23][24][25] which investigate the use of wavelets to discretize quantum field theories, in that we use biorthogonal wavelets, which moreover are specifically designed to target the Hamiltonian of the field theory. There is also a different approach to entanglement renormalization for quantum field theories, known as cMERA [34,35]. This takes a different perspective by formulating a variational class of states directly in the continuum, rather than considering a discretization. In both cases, the correlation functions of the theory are accurately reproduced up to some cut-off. The precise relation between MERA and cMERA is not very well understood, for instance it is not clear that discretizing a cMERA state could yield a MERA. Intriguingly, cMERA is formally strongly reminiscent of the continuous wavelet transform (CWT). The continuous wavelet transform [14] can be defined for a much broader class of wavelet functions ψ, and if ψ is a biorthogonal wavelet the CWT can be discretized to a discrete wavelet transform. Reformulating cMERA as the second quantization of a CWT would therefore give a clear relationship between MERA and cMERA for free bosonic systems. A starting point could be the cMERA in [36], which reproduces some scaling dimensions exactly. However, the CWT appears to break some of the symplectic properties of the discrete biorthogonal wavelet transform and it remains an open problem to make this connection more explicit. Finally, another reason why the field theory limit of entanglement renormalization is of interest is its tentative relation to holography in theories of quantum gravity, as conjectured in [37]. The entanglement renormalization circuit can be thought of as mapping a system into one higher dimension by adding an additional 'scale' direction. An interpretation in terms of wavelets was proposed in [38] for fermions and extended to bosonic systems in [39].

Conclusion
In this work we have explained how Gaussian entanglement renormalization circuits can be naturally contructed from (and are in fact equivalent to) the second quantizations of biorthogonal wavelet transforms. There are a few technical aspects that would be interesting to study in more detail. First, one could carry out a fully rigorous analysis of the continuum limit discussed in Section 5, as in [19]. This poses some mathematical challenges. For example, if the system is not scale invariant then our notion of wavelet functions goes beyond the standard framework of wavelet theory, and one would have to identify suitable conditions on the filters that ensure that the wavelet and scaling functions as defined in Eq. (22) are well-behaved functions and that standard wavelet theory generalizes. Second, it would be desirable to identify conditions under which the procedure outlined in Appendix B is rigorously guaranteed to find good approximate solutions of Eq. (11). We note that even for Hilbert pair wavelets, which are relevant in the fermionic setting and which inspired our construction, this is not known and a subject of recent research in the signal processing community [40,41].
Overall, we believe that this work, together with [16,17] for the fermionic case, completes our conceptual understanding of Gaussian entanglement renormalization for free theories as the second quantization of wavelet decompositions. We hope that this offers a path towards constructing and analyzing entanglement renormalization circuits for interacting models. One clear direction is to apply perturbation theory in the wavelet basis. A similar approach has already been taken for cMERA in [42], where one can also do perturbation theory around a Gaussian cMERA. Another interesting direction is to investigate integrable models, where we know explicit solutions for the ground state, and try to formulate these in terms of wavelet modes. Finally, continuous wavelet transforms might also help illuminate the relation between MERA and cMERA as we discussed in Section 5.4.

Acknowledgements
We thank Adrián Franco Rubio for inspiring discussions during the inception of this project.

A Review of the fermionic wavelet-MERA correspondence
In this appendix we briefly review the construction of entanglement renormalization circuits for massless free fermions, as worked out in [16,17,19]. While not strictly needed to understand the results of the current work, which deals with free bosons, it is instructive to contrast the construction and state of the art with the fermionic setting. We closely follow the exposition in [17]. Let a n for n ∈ be fermionic operators, satisfying the anticommutation relations {a † n , a m } = δ n,m , {a n , a m } = {a † n , a † m } = 0. We work in the framework of Gaussian or free fermions. We consider Hamiltonians of the form where h is a Hermitian matrix. The ground state of such a Hamiltonian is given by where |Ω〉 is the Fock vacuum and a † ( f ) := n f [n]a † n and where the f i are all eigenvectors of h with negative eigenvalue (so m h n,m f i [m] = λ i f i [n] with λ i < 0). In other words, precisely the negative energy modes are occupied in the ground state. The set of number-preserving Gaussian operations is given by all evolutions along Hamiltonians of the form in Eq. (23). Such transformations are the fermionic second quantization of unitaries acting on the single-particle space. That is, to every unitary operator U acting on 2 ( ) we associate the unitary which maps a n toã n = m u n,m a m .
We restrict to the nearest neighbor hopping Hamiltonian As opposed to the current work on bosonic models, so far no general wavelet construction for arbitrary free fermion models is known, but only for the Hamiltonian in Eq. (24) and its higher dimensional generalizations. To prepare the ground state of this Hamiltonian we separately consider the even and the odd sublattice and let a 1,n = a 2n and a 2,n = a 2n+1 , and we apply a phase gate, writing b i,n = (−1) n a i,n . The Hamiltonian in Eq. (24) is then transformed to This Hamiltonian can be written in Fourier space as The ground state |ψ〉 is now given by filling the negative energy modes (the Fermi sea), which can be found by observing that for each k the matrix 0 e −ik − 1 e ik − 1 0 has a negative eigenvalue with eigenvector which means that the space of negative energy modes consists of all functions f = ( f 1 , f 2 ) in 2 ( ) ⊗ 2 which are such that their Fourier transforms satisfy f 2 (k) = − sgn(k)ie i k 2 f 1 (k). That is, if f is such a function, then (b 1 ( f 1 ) † + b 2 ( f 2 ) † )|ψ〉 = 0. We now consider a pair of orthogonal wavelet filters (note that in contrast to the bosonic case, this is not a pair of biorthogonal filters, but two filters which are each orthogonal) g w and h w , and we let W g and W h denote the corresponding wavelet decomposition maps, which are now orthogonal, i.e. W −1 a = W T a for a = g, h. In particular, this means that we can apply the fermionic second quantization of W h to the b 1 fermions and W g to the b 2 fermions. If the filters are such that for −π < k < π Gaussian unitaries a n → m U n,m a m for U unitary q n → A n,m q m and p n → A n,m p m with A −1 = B T Hamiltonian H = − n a † n a n+1 + a † n+1 a n H = 1 2 n∈ p 2 n + n,m∈ q n V n−m q m with dispersion relation ω(k) Wavelet filters g, h orthogonal wavelet filters (g, h) pair of biorthogonal wavelet filters Application of wavelet transform b 1,n = (−1) n a 2n , b 2,n = (−1) n a 2n+1 , and apply W h to the b 1 fermions and W g to the b 2 fermions

Disentangling circuit
Apply wavelet decomposition, then H on wavelet modes.
Apply squeezing to normalize dispersion relation, then apply the wavelet decomposition.

Continuum theory
Free Dirac fermion Free bosonic scalar field (for this will allow us to renormalize the ground state. To see this, consider any mode f = ( f 1 , f 2 ) in the Fermi sea, so f 2 (k) = − sgn(k)ie i k 2 f 1 (k). When we apply the wavelet transforms f i is mapped to ( f i,w , f i,s ), a wavelet and a scaling component. Then one can show that from Eq. (26) it follows that f 1,w = f 2,w and f 2,s (k) = − sgn(k)ie i k 2 f 1,s (k). We may now apply (the fermionic second quantization of) a Hadamard gate to the wavelet component, which then disentangles the wavelet component of the negative energy mode. This shows that if we take the the ground state |ψ〉 and we first apply the wavelet transforms W g and W h followed by H on the wavelet output, we map to |ψ〉 to itself on the scaling output and the product state n b † 1,n |Ω〉 (that is, the state in which all even sublattice modes b 1,n are filled and all odd sublattice modes b 2n are not filled) on the wavelet output. Thus we have implemented a layer of entanglement renormalization. One can write the fermionic second quantization of an orthogonal wavelet transform with a finite filter as a finite depth fermionic Gaussian circuit [18]. We may iteratively apply the renormalization to the scaling component to completely completely disentangle the ground state, and if we consider the circuit in the opposite direction then it maps layers of product states to the ground state of Eq. (25). One may think of this way of preparing the ground state as filling the Fermi sea layer by layer, now choosing a wavelet basis for the Fermi sea instead of the usual Fourier basis.
The relation Eq. (26) can not be satisfied exactly by a pair of finite filters, but it can be approximated. In [16] it was shown that a construction using Daubechies D4 filters already gives a good result, and in [17] a general method using known constructions from signal processing applications [40] was suggested, and it was also shown that a good approximation of Eq. (26) leads to a good approximation of the ground state (which inspired our Theorem 1). The Hamiltonian in Eq. (25) is in fact the Kogut-Susskind discretization of the Dirac fermion [43].
In [19] it has been worked out how the wavelet and scaling functions corresponding to the filters g and h have a natural interpretation in the quantum field theory, analogous to the discussion in Section 5. In Table 1 we give an overview of the analogies between the fermionic and bosonic case.

B Construction of filters
Next we will explain how to construct filter pairs that yield a good approximation of a given dispersion relation. Suppose we are given a dispersion relation ω(k). Let us assume that ω(k) = ω(−k). We would like to construct a biorthogonal filter pair (g, h) such that or equivalently h s (k) ≈ ω(k + π)g s (k) .
We will describe a general approach to this problem inspired by the Daubechies wavelet construction, similar to the construction of Hilbert pair wavelets due to Selesnick [40] which were previously used in the construction of fermionic MERAs [17]. For this, we start with a rational approximation where a and b are real finite symmetric sequences on [−L, L]. The approximation only has to be accurate around k = 0. We will make the following ansatz for the Fourier transform of the scaling filters where f (k) is the Fourier transform of a real finite sequence f [n] that still needs to be determined. The parameter K determines the number of vanishing moments of the biorthogonal wavelets, just as in the Daubechies wavelet construction. By construction, g s (k) and h s (k) are small near k = π, and Eq. (28) is satisfied. In order for Eq. (29) to generate biorthogonal wavelet filters, they need to satisfy the condition in Eq. (1) which translates to where s(k) = a(k)b(k)(2 cos( k 2 )) 2K . One may try to solve this by letting r(k) = f (k) f (−k). Then r should be taken as a solution to the linear system Now, if possible, we perform a spectral factorization r(k) = f (k) f (−k). A necessary and sufficient condition for this is that r(k) ≥ 0 for all k. Unfortunately, we do not know of a condition on a and b that guarantees this. The resulting filters (g, h) will have support of size 2M where M = K + 2L. Finally, in the scale-invariant case, the stability condition that will be required in Theorem 1 can be checked explicitly for compactly supported filters by looking at the operators P g and P h defined by on the space of polynomials of degree at most 2M with zero mean. The filters yield square integrable scaling functions and uniformly bounded wavelet decomposition maps if and only if the eigenvalues of P g and P h are smaller in absolute value than 2 (see [44] or Theorem 4.2 in [45]).
For the massless harmonic chain, one particular choice for a and b is given by where d[n] is a maximally flat all-pass filter with delay 1 4 of degree L [40], so it has the property that e −i Lk d(−k)/d(k) ≈ e −i k 2 on k ∈ (−π, π). In Fig. 6 we show the goodness of the approximation in Eqs. (27) and (28) as a function of K and L. The resulting filters and wavelets for K = 2, L = 4 are shown in Fig. 5. We remark that the construction in Eq. (30) is not necessarily optimal. From numerical evidence in Fig. 6 it appears that the accuracy of the approximation improves exponentially with increasing support [26]. An interesting open problem is to rigorously prove the existence of approximate solutions to Eq. (27) with (exponentially) improving approximation accuracy as the filter size increases.

C Construction of circuits from filters
We now discuss how to explicitly construct a Gaussian circuit from a given pair of biorthogonal wavelet filters, and show that any translation-invariant Gaussian circuit of the form of Fig. 3 always arises from such a filter pair. Motivated by the fermionic setting it has been extensively discussed in [18] how to construct unitary local circuits from orthogonal wavelet filters. The construction for biorthogonal wavelet filters is very similar and the symmetric case has already been discussed in [18], but for completeness we provide it here. Given a pair of biorthogonal filters (g, h) of support 2M we will construct a binary circuit of depth M that implements the wavelet decomposition map. We will assume that g and h are supported on [−M + 1, M ], which we can always achieve by a as illustrated for g in Fig. 7. Now A ⊕ (A T ) −1 is a binary circuit, and its second quantization gives a Gaussian bosonic quantum circuit. It remains to construct the matrices a i given the filters g and h. We will need the perfect reconstruction condition Eq. (1) which becomes which follows directly from Eq. (2). First suppose that M = 1. In that case we let this case we may choose the sequence of filters to eventually become constant. For a = g, h and L ∈ , we define the L-layer decomposition maps We assume that the family is stable in the sense that the corresponding (generalized) scaling functions φ a l defined in Eq. (22) exist, are square integrable, and bounded in L ∞ -norm. We can also define the wavelet decomposition maps starting at layer L ≥ 0, that is, For L = 0 we recover W (L) a as defined earlier. We assume that the wavelet decomposition maps are bounded. Finally, we shall assume that the filters have finite support. Then the same is true for the scaling functions. In the case that the filters are independent of l, the above notion of stability is equivalent to the familiar notion from wavelet theory. For finitely supported filters there exists an easy criterion to determine this, see [31].
For the entanglement renormalization circuit, we also insert a squeezing operation between each wavelet decomposition layer, defining R a (l) , R (L) a and R L for a = g, h as in Eq. (10). Our approximation to the covariance matrix is then given by , with dispersion relation ω(k) such that ω (l) (π) ≤ 1 and ω (l) (k) ≤ Ω for l = 1, . . . , L for some Ω ≥ 1. Suppose we have a sequence of filters such that Eq. (34) holds for ≤ 1, with finite support of size at most M and scaling functions that are uniformly bounded by φ a l ∞ ≤ B for a = g, h and l = 1, . . . , L. Assume moreover that the wavelet decomposition maps are uniformly bounded by W (l ,l) a ≤ D for all a = g, h,g and 1 ≤ l ≤ l ≤ L, where D ≥ 1. Then the approximation error of the covariance matrices can be bounded as follows: where C := 4B 2 M 3 2 Ω.
To interpret the error bounds, we note that Now we can proceed as in the proof of Lemma 1 in [17] and estimate where in the second line we use that φ b is compactly supported on [x 0 , x 0 + M − 1] for some x 0 , in the third inequality we use Cauchy-Schwarz, and for the final inquality we use that at most 2M terms in the sum have nonzero overlap. Finally we may estimate φ b 0 2 ≤ M B 2 and φ a L 2 ∞ ≤ B 2 , which yields Eq. (41).
The following lemma bounds the approximation error for L layers as a function of an intermediate layer L that will later be chosen appropriately.
To interpret these bounds, we note that γ p,(L) = max k ω (L) (k) 2 , which is typically O(1). As a remark, for the critical harmonic chain we that ω (l) (π) = 1 2 , in which case it is not hard to see that the scaling of Eq.