aainstitutetext: Theoretical Division, Los Alamos National Laboratory, Los
Alamos, NM 87544, USA

# Entanglement Renormalization for Quantum Field Theories with Discrete
Wavelet Transforms

Daniele S. M. Alves<EMAIL_ADDRESS>
###### Abstract

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

††preprint: LA-UR-24-23540

##  1 Background: Entanglement Renormalization

Ongoing advances in quantum technologies have opened up the possibility that
quantum field theories (QFTs) might become simulatable in quantum computers in
the medium-to-long term future, offering new methods to address long-standing
and intractable problems that we have no hope to solve with classical
computers. With such prospects, much thought and research efforts have been
recently directed towards developing tools that are deliberately bypassed in
the traditional, Euclidean space treatments of QFTs. Such tools include the
construction of explicit Hilbert space representations of ground and excited
states of QFTs—including their entanglement patterns—and real-time Hamiltonian
evolution.

In building and understanding these tools, it pays off to borrow some of the
wisdom developed in condensed matter physics, which has traditionally used the
framework of wavefunctions, entanglement, and coarse graining to understand
the properties and dynamics of condensed matter systems. In particular, a
significant insight that emerged in the early 90s was that entanglement plays
a significantly role in how fundamental degrees of freedom reorganize
themselves as effective degrees of freedom under Renormalization Group (RG)
evolution.

Until the early 90s, there were notorious difficulties in numerically
implementing real-space block-spin renormalization of quantum spin chains. In
block-spin renormalization, a coarse-graining step consists of mapping a block
of sites of the fine-grained lattice into a single site of the coarse-grained
lattice via a decimation procedure that reduces the dimensionality of the
Hilbert space. In Wilson’s numerical RG Wilson:1974mb , this is implemented
via diagonalization of the Hamiltonian restricted to the block’s Hilbert space
and its subsequent truncation, keeping only the lowest energy eigenstates of
the block. This procedure, however, generically fails to retain the whole
support of the ground state, and provides rather unreliable results for a
variety of lattice systems PhysRevB.19.4876 ; XIANG1992861 ;
PhysRevLett.42.1492 , with the exception of Hamiltonians which can be put into
a special form. In 1992, White and Noack White:1992zza  understood the
fundamental difficulty with this approach: by neglecting interactions of the
block with the rest of the lattice, Wilson’s numerical RG fails to capture
entanglement between neighboring blocks, and therefore cannot preserve ground
state properties, regardless of any reasonable increase in the number of
states kept.

The shortcomings of the standard numerical RG approach were overcome by
White’s Density Matrix Renormalization Group (DMRG) White:1992zz . Instead of
the block’s Hamiltonian spectrum, this approach considered the spectrum of the
block’s reduced density matrix, and truncated its Hilbert space by discarding
eigenstates of the reduced density matrix with eigenvalues smaller than a
given threshold, determined by the desired level of accuracy. In this way, the
entanglement properties of the ground state could be retained to any desired
precision and a proper RG flow could be generated.

Eventually, the DMRG was understood in terms of the quantum information
theoretic concept of tensor networks (TNs), and recast as a variational
optimization method over the class of TNs known as Matrix Product States (MPS)
Schollwock11 ; Verstraete:2008 . Essentially, TNs represent quantum many-body
states in terms of networks of interconnected tensors, which in turn capture
the relevant entanglement properties of the wave function. TNs provide an
unbiased approach to numerical quantum many-body physics; they have been shown
to be highly efficient and accurate tools to compute physical observables,
while also being free from the sign problem that plagues Monte Carlo sampling
methods in Euclidean space (such as in, e.g., Lattice QFT), and also suitable
to simulate real-time dynamics.

The success of TNs in describing properties of quantum many-body systems can
be attributed to the efficiency of optimized TNs in identifying and keeping
track of the relevant degrees of freedom. This, in turn, is directly related
to the fact that the low energy states of realistic, gapped local Hamiltonians
are not generic states in the Hilbert space: they are heavily constrained by
locality, and their entanglement is limited by area laws. In particular, a key
insight of TNs is that, for systems described by local Hamiltonians, coarse
graining must be preceded by the removal of short-range entanglement by local
unitary transformations, so that short-distance degrees of freedom can be
safely traced out while preserving ground state entanglement properties
White:1992zza ; Vidal:2007hda . A particular class of TNs in which this is
self evident and intuitive is called MERA, an acronym for Multiscale
Entanglement Renormalization Ansatz Evenbly:2007hxg .

In MERA, short range entanglement across the boundaries of adjacent blocks is
first removed by local unitary transformations, or disentanglers, and then a
coarse-graining transformation can be applied to trace out short distance
degrees of freedom while safely preserving ground state entanglement
properties. An important consequence of the use of disentanglers is that the
dimension of the effective sites no longer needs to grow with each coarse
graining step, and therefore an efficient description of critical fixed points
becomes feasible 2011arXiv1109.5334E ; Pfeifer09 . MERA has been successfully
applied to study strongly interacting systems in low dimensions, including the
emergence of symmetry breaking order and topological order Aguado08 ;
PhysRevB.79.195123 , quantum critical points 2008PhRvL.101r0503G ; Pfeifer09 ,
real time evolution, fermionic and anyonic systems (which are intractable by
quantum Monte Carlo) 2010PhRvA..81a0303C ; 2010PhRvA..81e0303P , as well as a
lattice realization of the holographic principle Swingle:2009bg ;
Swingle:2012wq .

Figure 1: A (non-technical) depiction of the MERA architecture. From bottom to
top, it can be viewed as an RG flow, whereas from top to bottom, it can be
viewed as the construction of the ground state ansatz described in the text.

The construction of a MERA to represent the ground state wavefunction of a
system essentially amounts to a reverse RG flow. Consider, as a matter of
example, a lattice system governed by a local Hamiltonian. The starting point
of the MERA is an unentangled state of $n$ lattice sites, e.g., the product
state $|P\,\rangle=|0\rangle^{\otimes n}$. The reverse RG flow starts by
acting on $|P\,\rangle$ with quasi-local unitary transformations111Meaning
transformations that act on neighboring sites separated by a maximum finite
distance. to create entanglement between neighboring sites. These quasi-local
unitaries are called entanglers, and their collective action over the entire
system is represented by an operator $U_{0}$ such that

$|\Psi_{0}\rangle~{}\equiv~{}U_{0\,}|P\,\rangle~{}=~{}U_{0}\Big{[}|0\rangle^{\otimes
n}\Big{]}.$ (1)

A _“fine-graining”_ mapping then follows by the introduction of $n$ new
lattice sites intercalated with the original sites, effectively reducing the
lattice spacing by half. These new sites are also initialized as a product
state, so that the augmented, finer-lattice state is given by

$|\Psi_{0}^{\prime}\rangle~{}\equiv~{}W_{1}|\Psi_{0}\rangle~{}\equiv~{}\Big{[}|0\rangle^{\otimes
n}\Big{]}\otimes|\Psi_{0}\rangle.$ (2)

The next step is to entangle the newly introduced sites with the rest of the
lattice via new quasi-local entanglers represented by $U_{1}$, yielding

$|\Psi_{1}\rangle~{}=~{}U_{1}|\Psi_{0}^{\prime}\rangle.$ (3)

New layers of fine-graining and entangling operations can then be successively
applied to flow deeper into the UV, until a ground state ansatz at the desired
UV resolution (or until a fixed point) is reached, namely,

$|\Psi_{\text{\tiny UV}}\rangle~{}=\,\left[\prod_{j=0}^{j_{\text{\tiny
UV}}}U_{\\!j}W_{\\!j}\,\right]\\!\\!|P\,\rangle.$ (4)

The parameters of the entanglers $U_{\\!j}$ are variationally optimized (by
minimizing the expectation value of the Hamiltonian) to recover the ground
state of the lattice system. A (non-technical) depiction of the MERA described
above is shown in figure 1.

In the example described above, the entanglement introduced between adjacent
sites at scale $j$ gives rise to entanglement between $2^{(j_{\text{\tiny
UV}}-j)}$ sites in the final ansatz $|\Psi_{\text{\tiny UV}}\rangle$. Hence,
this iterative approach of introducing quasi-local entanglement at each scale
allows one to capture long distance correlations with a series of manageable,
quasi-local unitary operations. In essence, MERA works because (i)
entanglement is quasi-local in position space (this stems from interactions
being local) and (ii) the entanglement between effective degrees of freedom at
a given resolution is insensitive to shorter-range entanglement at much higher
resolutions (this stems from the decoupling principle Appelquist:1974tg ,
i.e., the insensitivity of a low energy effective theory to its UV
completion).

The general properties aforementioned are not exclusive to condensed matter
systems described by lattices; indeed, locality of interactions and decoupling
are also properties of continuum relativistic systems such as local QFTs.
Hence, one would expect that the MERA principles could be generalized to build
variational classes for QFTs. And indeed, in 2011, Haegeman et al.
Haegeman:2011uy  proposed such a generalization, which they coined
_“continuum_ MERA,” or cMERA for short. The elements of the original, discrete
MERA were conceptually modified in their adaptation to the continuum. In
particular, Haegeman et al. did not propose a discretization of the field
theory to more closely parallel the lattice structure and countable layers of
the original MERA; instead, cMERA was formulated directly in the continuum,
with the system’s resolution determined by a Wilsonian UV cut-off
$\Lambda_{\text{\tiny UV}}$. In lieu of a discrete layered architecture,
continuous scaling transformations were used to flow the system from one scale
to another. In particular, the role of the fine-graining mapping was formally
implemented by infinitesimal dilation operations to increase the UV cut-off
$\Lambda_{\text{\tiny UV}}$ under RG flow. Finally, the action of the unitary
entangler operator was implemented in a basis of momentum-space degrees of
freedom, effectively generating correlations at all scales up to the shortest
distance cut-off, $1/\Lambda_{\text{\tiny UV}}$. While formally attractive and
initially promising, the cMERA formulation turned out to be unwieldy and did
not lead to any notable results, except in examples of free field theories
Franco-Rubio:2017tkt ; zou2019magic . Specifically, the path-ordered dilation
transformations required to implement changes in scale proved to be
impractical for any non-trivial (i.e., interacting) field theory. More
problematic, however, was the fact that the implementation of entangling
operations in momentum space obscured the locality of entanglement and failed
to explicitly capture a key feature of MERA, namely, the hierarchical pattern
of entanglement between localized, neighboring degrees of freedom at each
scale.

The shortcomings of the cMERA formulation mentioned above were never addressed
in the many attempts to implement and/or generalize it to interacting field
theories Hu:2017rsp ; Cotler:2018ehb ; Franco-Rubio:2019nne ; Fernandez-
Melgarejo:2019sjo ; Fernandez-Melgarejo:2020fzw . This hindered any meaningful
progress in the study of QFTs with cMERA since it was first proposed over a
decade ago. In particular, no useful/practical entangler ansatzes for
interacting QFTs have been fully demonstrated to date.

In this article, we re-examine the generalization of the MERA principles for
field theories, and propose an alternative to cMERA which preserves, in
explicit form, key features of Entanglement Renormalization, namely, (i)
locality of entanglement and (ii) the layered architecture of the original
MERA. We coin this alternative MERA adaptation to QFTs “wMERA.” A key
underlying principle of wMERA is the use of Discrete Wavelet Transforms (DWT)
to discretize continuum quantum fields into countable degrees of freedom
localized in position space222We are not the first to use DWTs in the context
of MERA (see also Brennen:2014iqu ; Evenbly:2016cly ; Haegeman:2017vrx ;
Witteveen:2019lsk ; Witteveen:2020hdp ; Stottmeister:2020ezd ; George:2022ysn
; Sewell:2022nwi ). While there are some overlapping elements of our proposed
wMERA framework with previous works in the literature, our framework is built
with the explicit goal of making the generalization of MERA to interacting
QFTs more straightforward. In addition, to the best of our knowledge, we are
the first to use this framework to study _dynamics_ , i.e., to perform real-
time evolution and compute correlation functions at finite timelike
separations.. Intuitively, DWTs encode low- and high-pass filters that, upon a
choice of scale, separate the degrees of freedom of a quantum field into
short-range fluctuations and long-range fluctuations. This simple basis choice
provides a natural framework for RG evolution, while keeping locality and
separation of scales in explicit form, which makes the MERA generalization
quite straightforward. In particular, the fine-graining transformations are
implemented via Inverse Wavelet Transforms, which add shorter-range (wavelet)
degrees of freedom at each inverse RG step to flow the system to the UV.

The outline of this article is as follows. In section 2, we provide a brief
overview of cMERA and highlight its main shortcomings. In section 3, we
introduce the concept and basic machinery of discrete wavelet transforms, and
the specific class of Daubechies wavelets. In section 4, we provide the
generic formulation of wMERA using DWTs, the concrete implementation of which
is then illustrated in section 5 for two concrete examples, namely, a free
scalar field and a free Dirac fermion in (1+1)d. Our concluding remarks and
outlook are presented in section 6.

##  2 Brief overview of cMERA

In 2011, Haegeman et al. Haegeman:2011uy  introduced cMERA as a formal
adaptation of MERA to the continuum. Briefly, in cMERA the ground state
wavefunction of a system at its infrared (IR) fixed point,
$|\Omega_{\text{IR}}\rangle$, is related to its microscopic counterpart at an
ultraviolet (UV) scale $s_{\text{UV}}$, $|\Psi(s_{\text{UV}})\rangle$, via the
following unitary transformation:

$|\Psi(s_{\text{UV}})\rangle~{}=~{}\mathcal{P}\left[e^{-i\int_{{}_{s_{\text{\tiny
IR}}}}^{{}^{s_{\text{\tiny
UV}}}}~{}ds~{}\left(\hat{K}(s)\,+\,\hat{L}\;\right)}\right]\;|\Omega_{\text{IR}}\rangle\,.$
(5)

Above, $\mathcal{P}$ denotes path-ordering in scale $s$, $\hat{L}$ is the
generator of scaling transformations, and $\hat{K}(s)$ is the generator of
entanglement at scale $s$. Imposing renormalization group (RG) invariance of
physical observables,

$\langle\Psi(s_{\text{UV}})|\,\hat{O}(s_{\text{UV}})\,|\Psi(s_{\text{UV}})\rangle~{}=~{}\langle\Omega_{\text{IR}}|\,\hat{O}(s_{\text{IR}})\,|\Omega_{\text{IR}}\rangle\,,$
(6)

they obtained the RG equation for a generic hermitian operator $\hat{O}(s)$:

$\frac{\;d\hat{O}(s)}{ds}~{}=~{}-i\,\big{[}\hat{K}(s)+\hat{L}\,,\;\hat{O}(s)\big{]}.$
(7)

Haegeman et al. Haegeman:2011uy  illustrated the cMERA construction for a few
examples of free field theories333These examples are shown in detail in the
appendix of _version 1_ of their preprint posted on the arXiv online
repository.. For instance, for a free scalar field theory in $(d+1)$
dimensions, the entangling exponent in (5), $\hat{K}(s)+\hat{L}$, was given by
the following Gaussian generators:

$\displaystyle\hat{L}$ $\displaystyle=$
$\displaystyle\int_{0}^{\Lambda}d^{d}k~{}\frac{1}{2}\Bigg{[}\,\hat{\pi}(-\vec{k})\,\big{(}\vec{k}\\!\cdot\\!\vec{\nabla}_{k\,}\hat{\phi}(\vec{k})\big{)}\,+\,\big{(}\vec{k}\\!\cdot\\!\vec{\nabla}_{k\,}\hat{\phi}(\vec{k})\big{)}\,\hat{\pi}(-\vec{k})$
(8)
$\displaystyle~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}+\frac{d+1}{2}\,\left(\,\hat{\phi}(\vec{k})\,\hat{\pi}(-\vec{k})\,+\,\hat{\pi}(-\vec{k})\,\hat{\phi}(\vec{k})\,\right)\,\Bigg{]}\,,$

and

$\hat{K}(s)~{}=~{}-\int_{0}^{\Lambda}d^{d}k~{}\frac{1}{4}\,\left[1+\left(\frac{\Lambda}{m_{\phi}}e^{s}\right)^{2}\right]^{-1}\\!\left(\,\hat{\phi}(\vec{k})\,\hat{\pi}(-\vec{k})\,+\,\hat{\pi}(-\vec{k})\,\hat{\phi}(\vec{k})\,\right),$
(9)

where $\hat{\phi}(\vec{k})$ and $\hat{\pi}(\vec{k})$ are the momentum-space
modes of the scalar field and its canonical conjugate momentum, respectively;
$m_{\phi}$ is the mass of the scalar field; and $\Lambda$ is an _ad hoc_ UV
cut-off.

While formally elegant and concise, the cMERA formulation turned out to be
cumbersome even for the simplest case of free field theories. As we shall see
in section 5, the _ad hoc_ cut-off $\Lambda$ and the path-ordered architecture
of scaling transformations in (5) is completely unnecessary for free field
theories. More importantly, however, they are impractical for any nontrivial
(i.e., interacting) field theory, for which the entangler generator
$\hat{K}(s)$ must be generalized beyond quadratic order in the canonical
fields to be able to capture nontrivial $n$-point correlation functions beyond
the two-point function.

Compounded by these obstacles is the fact that all incarnations of cMERA
attempted so far in the literature were, to the best of this author’s
knowledge, implemented in momentum space. Since interactions are not local in
momentum space, neither is entanglement. More problematic is the fact that
interactions entangle momentum-space degrees of freedom across vastly
different scales. Hence, at the conceptual level, cMERA constructions in
momentum space obscure the hierarchical pattern of entanglement between
localized degrees of freedom in position space at each scale. At the practical
level, momentum-space cMERA is difficult to implement numerically, since it
does not share the original MERA’s layered architecture, nor the quasi-
locality of its entanglers.

To address these obstacles, we propose an alternative adaptation of MERA to
the continuum that is directly formulated in position space, and that does
away with the use of continuous path-ordered scaling transformations. In our
proposed formulation, which we coin “wMERA,” the layers of the MERA are
obtained via a Multiresolution Analysis, whereby the quantum fields are
expanded in a wavelet basis. There is several practical advantages of wMERA
over cMERA: (i) the wavelet discretization of quantum fields automatically
regulates the theory in the UV; (ii) the locality of degrees of freedom
preserves the quasi-locality of entangler ansatzes; (iii) separation of scales
is explicitly built-in, enabling the layered architecture of MERA to capture
the hierarchical pattern of entanglement across different scales; and (iv)
with the degrees of freedom and scales automatically discretized, numerical
methods developed for the original MERA can be adapted to wMERA.

Before describing wMERA in detail, we first we provide a brief introduction to
Discrete Wavelet Transforms and, in particular, to Multiresolution Analysis
with Daubechies wavelets, in section 3.

##  3 Brief overview of DWTs with Daubechies wavelets

Wavelet transforms offer a flexible framework to deconstruct the features of a
function, image, or sequence in a hierarchy of resolutions. They are
extensively used in signal processing and data compression.

In quantum field theories, on the other hand, Fourier Transforms (FTs) have
been traditionally used to decompose the degrees of freedom of a quantum field
into momentum modes. While such decompositions are optimal for problems
dealing with scattering of plane waves, the delocalized nature of Fourier
modes makes FTs highly inefficient in representing spatially localized
features, and, by extension, the pattern of ground state entanglement of
local, interacting field theories. A suitable alternative to FTs, which we
will use in our wMERA formulation, are discrete wavelet transforms (DWTs),
which will allow us to decompose the degrees of freedom of a quantum field
into localized modes in a hierarchy of resolutions444The discussion in this
section is heavily based on Bulut & Polyzou Bulut:2013bg ..

In order to intuitively understand some of the key features of DWTs, we begin
by describing Haar wavelets, which are arguably the simplest wavelet bases
there are.

###  3.1 Haar wavelets

The Haar wavelet bases are generated from discrete translations and rescalings
of the Haar scaling function, which is essentially the unit box function
centered around $x=1/2$:

$s(x)=\begin{cases}1~{}~{}\text{if}~{}&0\leq x<1\\\
0~{}&\text{otherwise.}\end{cases}$ (10)

Suppose one wishes to represent square integrable functions on the real line,
$L^{2}(\mathbb{R})$, with a resolution of $\Delta x=1$. An orthonormal basis
for such decomposition can be constructed from discrete translations of the
Haar scaling function $s(x)$:

$\\{s_{n}(x)\\}\big{|}_{n\in\mathbb{Z}}~{}~{}~{}~{}\text{with}~{}~{}~{}~{}~{}s_{n}(x)~{}\equiv~{}{\mathcal{T}}^{\,n}[s(x)]~{}\equiv~{}s(x-n).$
(11)

Above, ${\mathcal{T}}$ is the unit translation operator.

Figure 2: Illustration of the decomposition of a generic function $f(x)$ in
three Haar scaling function bases with different resolutions: $\Delta x=a_{0}$
(left), $\Delta x=a_{0}/2$ (middle), and $\Delta x=a_{0}/4$ (right).

Now, suppose one wishes to improve the resolution of the decomposition basis
by a factor of 2, i.e., to $\Delta x=1/2$. In this case, the appropriate Haar
basis would consist of rescaled scaling functions which are narrower by a
factor of 2 (to improve resolution) and taller by a factor of $\sqrt{2}$ (to
preserve unit normalization):

$\displaystyle\\{s_{n}^{1}(x)\\}\big{|}_{n\in\mathbb{Z}}~{}~{}~{}~{}\text{with}~{}~{}~{}~{}~{}s_{n}^{1}(x)~{}\equiv~{}\mathcal{D}[s_{n}(x)]$
$\displaystyle\equiv$ $\displaystyle\sqrt{2}\;s_{n}(2x)$ (12) $\displaystyle=$
$\displaystyle\sqrt{2}\;s(2x-n).$

Above, ${\mathcal{D}}$ is the dyadic rescaling operator. An example in which a
generic function is decomposed in this basis in a sequence of improved
resolutions is shown in figure 2.

It is easy to see that the space $\mathcal{S}_{0}$ spanned by the original
basis $\\{s_{n}(x)\\}$ is a subset of the space $\mathcal{S}_{1}$ spanned by
the improved-resolution basis $\\{s^{1}_{n}(x)\\}$:

$\mathcal{S}_{0}~{}\subset~{}\mathcal{S}_{1}.$ (13)

The orthogonal complement of $\mathcal{S}_{0}$ in $\mathcal{S}_{1}$ is a
nonempty set called the wavelet subspace $\mathcal{W}_{0}$:

$\mathcal{S}_{1}~{}=~{}\mathcal{S}_{0}\oplus\mathcal{W}_{0}~{},~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\mathcal{W}_{0}\neq\emptyset\,.$
(14)

$\mathcal{W}_{0}$ is spanned by the Haar wavelet functions $w_{n}(x)$, which
are generated from discrete translations of the _mother wavelet_ $w(x)$:

$w_{n}(x)\equiv{\mathcal{T}}^{\,n}[w(x)]=w(x-n)~{}~{}~{}~{}~{}~{}\text{with}~{}~{}~{}~{}~{}~{}w(x)=\begin{cases}~{}~{}\,1~{}~{}\text{if}~{}&0\leq
x<\frac{1}{2}\\\ -1~{}~{}\text{if}~{}&\frac{1}{2}\leq x<1\\\
~{}~{}\,0~{}&\text{otherwise.}\end{cases}$ (15)

In other words, the bases

$\\{s_{n}^{1}(x)\\}\big{|}_{n\in\mathbb{Z}}~{}~{}~{}~{}~{}~{}\text{and}~{}~{}~{}~{}~{}~{}\\{s_{n}(x)\\}\big{|}_{n\in\mathbb{Z}}\,\cup\;\\{w_{n}(x)\\}\big{|}_{n\in\mathbb{Z}}$
(16)

span the same subspace $\mathcal{S}_{1}$ in $L^{2}(\mathbb{R})$. Their
elements are related by the following orthogonal transformation:

$\displaystyle s_{n}(x)$ $\displaystyle=$
$\displaystyle~{}~{}\frac{1}{\sqrt{2}}\Big{(}s_{2n}^{1}(x)+s_{2n+1}^{1}(x)\Big{)}$
(17a) $\displaystyle w_{n}(x)~{}$ $\displaystyle=$
$\displaystyle~{}~{}\frac{1}{\sqrt{2}}\Big{(}s_{2n}^{1}(x)-s_{2n+1}^{1}(x)\Big{)}$
(17b)

whose inverse is:

$\displaystyle s^{1}_{2n}(x)$ $\displaystyle=$
$\displaystyle~{}~{}\frac{1}{\sqrt{2}}\Big{(}s_{n}(x)+w_{n}(x)\Big{)}$ (18a)
$\displaystyle s^{1}_{2n+1}(x)~{}$ $\displaystyle=$
$\displaystyle~{}~{}\frac{1}{\sqrt{2}}\Big{(}s_{n}(x)-w_{n}(x)\Big{)}.$ (18b)

The transformation in (17) defines the Haar _wavelet transform_ (WT), which
can be recast in the following form:

$\displaystyle s_{n}(x)$ $\displaystyle=$
$\displaystyle~{}~{}\mathbf{L}_{nm}\,s_{m}^{1}(x),$ (19a) $\displaystyle
w_{n}(x)~{}$ $\displaystyle=$
$\displaystyle~{}~{}\mathbf{H}_{nm}\,s_{m}^{1}(x)$ (19b)

with

$\displaystyle\mathbf{L}_{nm}~{}\equiv~{}h_{m-2n}~{}~{}~{}~{}~{}\text{and}~{}~{}~{}~{}~{}h_{n}=\begin{cases}~{}\frac{1}{\sqrt{2}}~{}~{}\text{if}~{}&n=0\\\
~{}\frac{1}{\sqrt{2}}~{}~{}\text{if}~{}&n=1\\\
~{}0~{}&\text{otherwise.}\end{cases}$ (20)
$\displaystyle\mathbf{H}_{nm}~{}\equiv~{}g_{m-2n}~{}~{}~{}~{}~{}\text{and}~{}~{}~{}~{}~{}g_{n}=\begin{cases}~{}~{}\frac{1}{\sqrt{2}}~{}~{}\text{if}~{}&n=0\\\
-\frac{1}{\sqrt{2}}~{}~{}\text{if}~{}&n=1\\\
~{}~{}\,0~{}&\text{otherwise.}\end{cases}$ (21)

The projection operators $\mathbf{L}$ and $\mathbf{H}$ in (20) and (21) are,
respectively, the Haar low-pass filter and the Haar high-pass filter. They are
defined in terms of the finite array of weight coefficients $\\{h_{n}\\}$ and
$\\{g_{n}\\}$. With this notation, the _inverse wavelet transform_ (IWT) in
(18) can be rewritten as:

$s^{1}_{n}(x)~{}=~{}\mathbf{L}^{T}_{nm}\,s_{m}(x)+\mathbf{H}^{T}_{nm}\,w_{m}(x).$
(22)

The Haar wavelet transform and its inverse are depicted in figure 3.

Figure 3: _Left panel:_ a pictorial illustration of the Haar wavelet
transform, which decomposes a signal into low-resolution “collective” features
(green scaling functions) and high-resolution features (orange wavelet
functions). _Right panel:_ the inverse Haar wavelet transform, which
recombines this decomposition into a basis of finer-resolution collective
features.

If one wishes to represent functions with resolution $\Delta x=1/2^{r}$ (where
the resolution index $r$ is always an integer), a Haar basis can be
constructed from $r$ dyadic rescalings of $s_{n}(x)$, yielding basis function
elements with width equal to $1/2^{r}$:

$\displaystyle s_{n}^{r}(x)~{}\equiv~{}\mathcal{D}^{r}[s_{n}(x)]$
$\displaystyle=$ $\displaystyle 2^{r/2}\;s_{n}(2^{r}x)$ (23) $\displaystyle=$
$\displaystyle 2^{r/2}\;s(2^{r}x-n).$

In particular, the subspace $\mathcal{S}_{r+1}$ spanned by the basis
$\\{s_{n}^{r+1}(x)\\}\big{|}_{n\in\mathbb{Z}}$ contains all lower resolution
subspaces $\mathcal{S}_{r},\,\mathcal{S}_{r-1},\,...,\,\mathcal{S}_{0}$ in a
nested form:

$\mathcal{S}_{0}~{}\subset~{}\mathcal{S}_{1}~{}\subset~{}\mathcal{S}_{2}~{}\subset~{}...\subset~{}\mathcal{S}_{r}~{}\subset~{}\mathcal{S}_{r+1}~{}\subset~{}...\,.$
(24)

In full analogy with (14), the orthogonal complement of $\mathcal{S}_{r}$ in
$\mathcal{S}_{r+1}$ is the wavelet subspace $\mathcal{W}_{r}$:

$\mathcal{S}_{r+1}~{}=~{}\mathcal{S}_{r}\oplus\mathcal{W}_{r},$ (25)

which is spanned by a basis of rescaled wavelets:

$\displaystyle w_{n}^{r}(x)~{}\equiv~{}\mathcal{D}^{r}[w_{n}(x)]$
$\displaystyle=$ $\displaystyle 2^{r/2}\;w_{n}(2^{r}x)$ (26) $\displaystyle=$
$\displaystyle 2^{r/2}\;w(2^{r}x-n).$

The generalization of the (inverse) wavelet transform relations between the
two equivalent bases

$\\{s_{n}^{r+1}(x)\\}\big{|}_{n\in\mathbb{Z}}~{}~{}~{}~{}~{}~{}\text{and}~{}~{}~{}~{}~{}~{}\\{s^{r}_{n}(x)\\}\big{|}_{n\in\mathbb{Z}}\,\cup\;\\{w^{r}_{n}(x)\\}\big{|}_{n\in\mathbb{Z}}$
(27)

spanning $\mathcal{S}_{r+1}$ is, expectedly, similar to the $r=0$ expressions
in (19) and (22):

$\displaystyle s^{r}_{n}(x)$ $\displaystyle=$
$\displaystyle~{}~{}\mathbf{L}_{nm}\,s_{m}^{r+1}(x),$ (28a) $\displaystyle
w^{r}_{n}(x)~{}$ $\displaystyle=$
$\displaystyle~{}~{}\mathbf{H}_{nm}\,s_{m}^{r+1}(x),$ (28b)

and

$s^{r+1}_{n}(x)~{}=~{}\mathbf{L}^{T}_{nm}\,s^{r}_{m}(x)+\mathbf{H}^{T}_{nm}\,w^{r}_{m}(x).$
(29)

Combining (24) and (25), we have a multiresolution decomposition of
$L^{2}(\mathbb{R})$:

$L^{2}(\mathbb{R})~{}=~{}\mathcal{S}_{r}\oplus\mathcal{W}_{r}\oplus\mathcal{W}_{r+1}\oplus\mathcal{W}_{r+2}\oplus...\oplus\mathcal{W}_{r+j}\oplus...,$
(30)

i.e., for a base resolution choice $r$, the set of functions

$\\{s_{n}^{r}(x)\\}\big{|}_{n\in\mathbb{Z}}\,\cup~{}\\{w_{n}^{r}(x)\\}\big{|}_{n\in\mathbb{Z}}\,\cup~{}\\{w_{n}^{r+1}(x)\\}\big{|}_{n\in\mathbb{Z}}\,\cup~{}\dots\cup\\{w_{n}^{r+j}(x)\\}\big{|}_{n\in\mathbb{Z}}\,\cup~{}\dots$
(31)

form an orthonormal basis for $L^{2}(\mathbb{R})$ with arbitrarily high
resolution555Note that scaling functions $s_{n}^{r}(x)$ and
$s_{n}^{r^{\prime}}(x)$ with different resolutions (i.e., $r\neq r^{\prime}$)
are not orthogonal to each other and are not elements of the same basis.
Likewise, wavelets $w_{n}^{\,r^{\prime}}(x)$ with $r^{\prime}<r$ are not
orthogonal to $s_{n}^{r}(x)$ and are not elements of the basis in (31).. The
scaling functions $s_{n}^{r}(x)$ encode structure up to scale $1/2^{r}$, and
the wavelets $w_{n}^{r+j}(x)$ capture structure at finer resolutions
$1/2^{r+j+1}$.

It is easy to see that the multiresolution decomposition in (30) and its
associated basis in (31) offer the desired architecture upon which to build
the layers of a MERA for a QFT ground state. The discretized degrees of
freedom of a quantum field $\phi^{r}_{i}$, obtained from smearing $\phi(x)$
over the scaling functions $s_{i}^{r}(x/a_{0})$, heuristically describe
collective field fluctuations over distances of order
$\mathcal{O}(a_{0}/2^{r})$. Shorter-distance fluctuations of the field at
higher resolutions $r^{\prime}>r$ are captured by $\phi(x)$ projections over
the wavelet functions $w_{n}^{r^{\prime}}(x/a_{0})$. Inverse wavelet
transforms are a means to recombine scaling- and wavelet-degrees of freedom
into shorter-distance collective modes, whereas wavelet transforms are a means
to tease out UV fluctuations that need to be disentangled to enable coarse-
graining.

Unfortunately, the Haar multiresolution decomposition has shortcomings that
make it not completely suitable for this intended MERA application. In
particular, the Haar scaling and wavelet functions are not differentiable,
rendering the kinetic term in the Hamiltonian for discretized field degrees of
freedom ill-defined. The naive lattice QFT discretization, which strongly
parallels the Haar discretization, goes around this problem by adopting the ad
hoc dictionary

$\frac{1}{2}\int
dx\;(\partial_{x}\phi)^{2}~{}~{}\to~{}~{}~{}\frac{1}{2}\,\sum_{ij}\,\phi^{r}_{i}\,\mathds{K}_{ij}^{r}\,\phi^{r}_{j}$
(32)

at resolution $\Delta x=a_{0}/2^{r}$, where $a_{0}$ is the reference lattice
spacing for $r=0$, and the discretized kinetic term coefficients are given by

$\mathds{K}_{ij}^{r}~{}\equiv~{}\frac{1}{(a_{0}/2^{r})^{2}}~{}\big{(}2\delta_{ij}-\delta_{i+1,j}-\delta_{i-1,j}\big{)}.$
(33)

The problem is that the “Haar smearing” interpretation of $\phi^{r}_{i}$ as

$\phi^{r}_{i}\equiv\int dx\;\phi(x)\,s_{i}^{r}(x/a_{0})$ (34)

is inconsistent with the scaling dependence of the discretized kinetic term in
(33). Specifically, consistency requires that coarse graining $\mathds{K}^{r}$
(through low-pass filter projections) should recover the kinetic term at lower
resolutions. However, that is not the case. One can verify by inspection that

$\big{(}\mathbf{L}\,\mathds{K}^{r}\,\mathbf{L}^{\mathsf{T}}\big{)}_{ij}~{}=~{}2\times\,\mathds{K}_{ij}^{(r-1)}$
(35)

which is the incorrect scaling.

This problem is specific to wavelet classes whose basis elements are not
differentiable—not a fatal flaw of Multiresolution Analysis in general. It can
be bypassed by adopting a class of wavelets with continuous first derivatives.
Indeed, Haar wavelets belong to a broader family of wavelets called Daubechies
wavelets, named after their inventor Ingrid Daubechies 10.1002/cpa.3160410705
. Other types of wavelets within the Daubechies family retain many of the
simple and attractive properties of Haar wavelets while satisfying our
differentiability requirement, as we shall discuss shortly. In what follows,
briefly review of the definition and properties of Daubechies wavelets.

###  3.2 Daubechies D$N$ wavelets

Daubechies wavelets provide a family of orthonormal bases for all square
integrable functions on the real line, $L^{2}(\mathbb{R})$. Each type of
Daubechies wavelets is characterized their number $N$ of non-vanishing weights
$\\{h_{j}\\}$—hence the nomenclature Daubechies “D$N$” wavelets. The Haar
wavelets, which have $N=2$ non-vanishing weights $h_{0}=h_{1}=1/\sqrt{2}$ (see
(20)), are also called Daubechies D2 wavelets.

The Daubechies D$N$ basis elements are functions with compact support at
different spatial resolutions, including basis elements which vanish outside
an arbitrarily small open set. They are fractal functions, since any element
of the basis can be generated from a single scaling function $s(x)$ by
applying a finite number of discrete translations and rescalings to $s(x)$,
and forming linear superpositions of these rescaled and translated functions.
In particular, the D$N$ scaling function satisfies a linear renormalization
group equation:

$s(x)~{}=~{}\mathcal{D}\\!\left[\sum_{j=0}^{N-1}~{}h_{j}\,\mathcal{T}^{j}[s(x)]\right]\\!,$
(36)

where the weights $\\{h_{j}\\}$ are a sequence of $N$ non-vanishing numbers
(to be discussed shortly), and the unit translation ($\mathcal{T}$) and dyadic
rescaling ($\mathcal{D}$) operators have been previously introduced in (11)
and (12). Essentially, $s(x)$ is the fixed point of an operation which takes
the weighted average of a finite number of translated copies of itself and
rescales it to half of its original support.

The scaling function $s(x)$ has compact support on the finite interval
$[0,N-1]$. Further scaling functions with finer resolution are obtained by
rescalings and translations of $s(x)$:

$s_{n}^{r}(x)~{}\equiv~{}\mathcal{D}^{r}\big{[}\mathcal{T}^{n}[s(x)]\big{]}\,.$
(37)

Indeed, the subspace spanned by $\\{s^{r}_{n}(x)\\}$, $\mathcal{S}_{r}$, has a
finest resolution scale of $1/2^{r}$ and contains all lower resolution
subspaces $\mathcal{S}_{r-1},\,\mathcal{S}_{r-2},\,...,\,\mathcal{S}_{0}$ in
the same nested form as in (24).

The wavelet subspaces are generated from the mother wavelet $w(x)$, which also
descends from the scaling function $s(x)$ via:

$w(x)~{}=~{}\sum_{j=0}^{N-1}~{}g_{j}\,\mathcal{D}\left[\mathcal{T}^{j}[s(x)]\right].$
(38)

Above, the weights $\\{g_{j}\\}$ are related to the weights $\\{h_{j}\\}$ in
(36) by:

$g_{j}~{}=~{}(-1)^{j}h_{N-1-j}.$ (39)

In full analogy with the case of Haar wavelets, a basis of translated and
rescaled mother wavelets $\\{w_{n}^{r}(x)\\}$ will span the D$N$ wavelet
subspace $\mathcal{W}_{r}$. A multiresolution decomposition of
$L^{2}(\mathbb{R})$ with Daubechies D$N$ wavelets then follows with the same
structure of (30) and (31).

Additionally, the action of the low- and high-pass filters in (inverse)
wavelet transforms for Daubechies D$N$ wavelets is formally the same as for
Haar wavelets:

WT
$\displaystyle\begin{cases}\,s^{r}_{n}(x)=~{}~{}\mathbf{L}_{nm}\,s_{m}^{r+1}(x),\\\
\,w^{r}_{n}(x)~{}=~{}\mathbf{H}_{nm}\,s_{m}^{r+1}(x),\end{cases}$ (40) IWT
$\displaystyle\begin{cases}\,s^{r+1}_{n}(x)~{}=~{}\mathbf{L}^{T}_{nm}\,s^{r}_{m}(x)+\mathbf{H}^{T}_{nm}\,w^{r}_{m}(x),\end{cases}$
(41)

with the difference that, above, the low-pass and high-pass filters
$\mathbf{L}$ and $\mathbf{H}$ are determined by the D$N$ weights that enter in
(36) and (39):

$\mathbf{L}_{nm}~{}\equiv~{}h_{m-2n}~{},~{}~{}~{}~{}~{}\mathbf{H}_{nm}~{}\equiv~{}g_{m-2n}.$
(42)

The weights $\\{h_{j}\\}$ for Daubechies D$N$ wavelets are determined by three
conditions:

(1)

Unit area for $s(x)$:

$\int dx~{}s(x)=1~{}~{}~{}\iff~{}~{}\sum_{j=0}^{N-1}h_{j}=\sqrt{2}\,,$ (43)

(2)

Orthonormality of the basis $\\{s_{n}(x)\\}$:

$\int
dx\,s_{n}(x)s_{m}(x)=\delta_{nm}~{}~{}~{}\iff~{}~{}\sum_{j=0}^{N-1}h_{j}\,h_{j-2n}=\delta_{n,0}\,,$
(44)

(3)

Ability of the basis $\\{s_{n}(x)\\}$ to locally represent $k$-degree
polynomials with $0\leq k\leq(N/2-1)$. In other words, there must exist a
sequence $\\{c_{n}\\}$ such that:

$x^{k}=\sum_{n=-\infty}^{\infty}\,c_{n}\,s_{n}(x).$ (45)

Given that wavelets and scaling functions are orthogonal to each other, the
property above is equivalent to the requirement that the wavelets possess
$N/2$ $-$$1$ vanishing moments:

$\int dx\;x^{k}\,w(x)=0~{}~{},~{}~{}~{}~{}0\leq k\leq(N/2-1).$ (46)

It then follows from (38) and (39) that this last condition on the weights
$\\{h_{j}\\}$ can be written as:

$\sum_{j=0}^{N-1}j^{n}\,g_{j}~{}=~{}\sum_{j=0}^{N-1}j^{n}(-1)^{j}\,h_{N-1-n}~{}=~{}0\,.$
(47)

There are two solutions to conditions (43), (44), and (47) above, namely,
$\\{h_{j}\\}$ and $\\{h_{j}^{\prime}\\}$. They are related by reverse
ordering, i.e., $h_{j}^{\prime}=h_{N-1-j}$, and their corresponding $s(x)$ and
$s^{\prime}(x)$ are mirror images of each other.

For our wMERA formulation, we would like to choose a Daubechies D$N$ basis
with a minimum requirement of a continuous first derivative. The smallest $N$
value for which this condition is satisfied is $N=6$ doi:10.1137/0522089 ;
doi:10.1137/0523059 . Bases of smoother functions with higher
differentiability could be chosen at the expense of additional computational
cost, since such bases would have a larger number of weights $\\{h_{j}\\}$.
For our present purposes, we will consider Hamiltonians exhibiting at most
first derivatives, and therefore, for the remainder of this paper, we will
adopt Daubechies D6 wavelets for our wMERA formulation.

Figure 4: The Daubechies D6 scaling function $s(x)$ and mother wavelet $w(x)$.

The weights for Daubechies D6 wavelets are given by:

$\displaystyle h_{0}$ $\displaystyle~{}=~{}$
$\displaystyle\frac{1}{16\sqrt{2}}\;\bigg{(}1+\sqrt{10}+\sqrt{5+2\sqrt{10}}\;\bigg{)}$
(48a) $\displaystyle h_{1}$ $\displaystyle~{}=~{}$
$\displaystyle\frac{1}{16\sqrt{2}}\;\bigg{(}5+\sqrt{10}+3\sqrt{5+2\sqrt{10}}\;\bigg{)}$
(48b) $\displaystyle h_{2}$ $\displaystyle~{}=~{}$
$\displaystyle\frac{1}{16\sqrt{2}}\;\bigg{(}10-2\sqrt{10}+2\sqrt{5+2\sqrt{10}}\;\bigg{)}$
(48c) $\displaystyle h_{3}$ $\displaystyle~{}=~{}$
$\displaystyle\frac{1}{16\sqrt{2}}\;\bigg{(}10-2\sqrt{10}-2\sqrt{5+2\sqrt{10}}\;\bigg{)}$
(48d) $\displaystyle h_{4}$ $\displaystyle~{}=~{}$
$\displaystyle\frac{1}{16\sqrt{2}}\;\bigg{(}5+\sqrt{10}-3\sqrt{5+2\sqrt{10}}\;\bigg{)}$
(48e) $\displaystyle h_{5}$ $\displaystyle~{}=~{}$
$\displaystyle\frac{1}{16\sqrt{2}}\;\bigg{(}1+\sqrt{10}-\sqrt{5+2\sqrt{10}}\;\bigg{)}$
(48f) $\displaystyle h_{j}$ $\displaystyle~{}=~{}$ $\displaystyle
0~{}~{}\text{if}~{}~{}j<0~{}\text{or}~{}j>5$ (48g)

Figure 4 shows the D6 scaling function $s(x)$ and mother wavelet $w(x)$.

In anticipation of our wMERA study of free scalar and fermion theories in
section 5, we end this section with a compilation of overlap integrals of D6
scaling and wavelet functions and their first derivatives, which will
determine the kinetic-term coefficients of the discretized (1+1)-d
Hamiltonians we will consider. For a derivation of the numerical values of the
integrals listed below, see, e.g., Polyzou:2020ifj .

For scalar fields, there will be three types of kinetic-term coefficients,
namely:

$\displaystyle\mathds{K}_{ij}^{{ss}\,(r)}~{}~{}~{}$ $\displaystyle\equiv$
$\displaystyle~{}\int\,dx~{}\nabla s_{i}^{r}(x)\,\nabla s^{r}_{j}(x)$ (49)
$\displaystyle\mathds{K}_{ij}^{{sw}\,(r)}~{}~{}$ $\displaystyle\equiv$
$\displaystyle~{}\int\,dx~{}\nabla s^{r}_{i}(x)\,\nabla w^{r}_{j}(x)$ (50)
$\displaystyle\mathds{K}_{ij}^{{ww}\,(r)}~{}$ $\displaystyle\equiv$
$\displaystyle~{}\int\,dx~{}\nabla w^{r}_{i}(x)\,\nabla w^{r}_{j}(x).$ (51)

Discrete translational invariance implies that these kinetic coefficients
depend only on the separation between lattice sites, $i-j$. In addition, the
compact support of $s_{n}^{r}(x)$ and $w_{n}^{r}(x)$ on the finite interval
$n/2^{r}\leq x\leq(n+5)/2^{r}$ implies that the kinetic coefficients vanish
for separations $|i-j|\geq 5$. Defining

$\mathds{K}_{ij}^{{ss}\,(r=0)}~{}\equiv~{}\mathds{K}_{ij}^{{ss}}~{}\equiv~{}K_{|i-j|},$
(52)

the nonzero ${K_{j}}$ coefficients are given by:

$K_{0}=\frac{295}{56}~{},~{}~{}~{}~{}K_{1}=-\frac{356}{105}~{},~{}~{}~{}~{}K_{2}=\frac{92}{105}~{},~{}~{}~{}~{}K_{3}=-\frac{4}{35}~{},~{}~{}~{}~{}K_{4}=-\frac{3}{560}\,,$
(53)

and illustrated in figure 5 (left).

Figure 5: The kinetic term coefficients for scalars (left) and fermions
(right).

The other kinetic coefficients can be obtained by applications of low- and
high-pass filters (defined through (42), (48), and (39)):

$\displaystyle\mathds{K}^{{ss}\,(r)}~{}~{}$ $\displaystyle\equiv~{}$
$\displaystyle 2^{2r}\,\mathds{K}^{{ss}}$ (54)
$\displaystyle\mathds{K}^{{sw}\,(r)}~{}$ $\displaystyle\equiv~{}$
$\displaystyle
2^{2(r+1)}\,\mathds{K}^{{sw}}~{},~{}~{}~{}~{}~{}~{}~{}\,\mathds{K}^{{sw}}_{ij}\\!~{}~{}\equiv~{}\big{(}\mathbf{L}\,\mathds{K}^{ss}\,\mathbf{H}^{\mathsf{T}}\big{)}_{ij}$
(55) $\displaystyle\mathds{K}^{{ww}\,(r)}$ $\displaystyle\equiv~{}$
$\displaystyle
2^{2(r+1)}\,\mathds{K}^{{ww}}~{},~{}~{}~{}~{}~{}~{}\,\mathds{K}^{{ww}}_{ij}~{}\equiv~{}\big{(}\mathbf{H}\,\mathds{K}^{ss}\,\mathbf{H}^{\mathsf{T}}\big{)}_{ij}.$
(56)

For a fermion field, the four types of kinetic coefficients will be:

$\displaystyle\mathds{D}_{ij}^{{ss}\,(r)}~{}~{}~{}$ $\displaystyle\equiv$
$\displaystyle~{}\int\,dx~{}s_{i}^{r}(x)\,\nabla s^{r}_{j}(x)$ (57)
$\displaystyle\mathds{D}_{ij}^{{sw}\,(r)}~{}~{}$ $\displaystyle\equiv$
$\displaystyle~{}\int\,dx~{}s^{r}_{i}(x)\,\nabla w^{r}_{j}(x)$ (58)
$\displaystyle\mathds{D}_{ij}^{{ws}\,(r)}~{}~{}$ $\displaystyle\equiv$
$\displaystyle~{}\int\,dx~{}w^{r}_{i}(x)\,\nabla
s^{r}_{j}(x)~{}=-\mathds{D}_{ji}^{{sw}\,(r)}$ (59)
$\displaystyle\mathds{D}_{ij}^{{ww}\,(r)}~{}$ $\displaystyle\equiv$
$\displaystyle~{}\int\,dx~{}w^{r}_{i}(x)\,\nabla w^{r}_{j}(x).$ (60)

Unlike the scalar field case with symmetric $\mathds{K}^{{ss}}$ kinetic
coefficients, in the fermionic case the $\mathds{D}^{{ss}}$ kinetic
coefficients are antisymmetric. Defining

$\mathds{D}_{ij}^{{ss}\,(r=0)}~{}\equiv~{}\mathds{D}^{{ss}}_{ij}~{}\equiv~{}\text{sign}(i-j)\,D_{|i-j|},$
(61)

the non-zero ${D_{j}}$ coefficients are given by

$D_{1}=-\frac{272}{365}~{},~{}~{}~{}~{}D_{2}=\frac{53}{365}~{},~{}~{}~{}~{}D_{3}=-\frac{16}{1095}~{},~{}~{}~{}~{}D_{4}=-\frac{1}{2920}\,,$
(62)

and illustrated in figure 5 (right). The other fermionic kinetic coefficients
can be obtained by applications of low- and high-pass filters:

$\displaystyle\mathds{D}^{{ss}\,(r)}~{}~{}$ $\displaystyle\equiv~{}$
$\displaystyle 2^{r}\,\mathds{D}^{{ss}}$ (63)
$\displaystyle\mathds{D}^{{sw}\,(r)}~{}$ $\displaystyle\equiv~{}$
$\displaystyle
2^{r+1}\,\mathds{D}^{{sw}}~{},~{}~{}~{}~{}~{}~{}~{}\,\mathds{D}^{{sw}}_{ij}\\!~{}~{}\equiv~{}\big{(}\mathbf{L}\,\mathds{D}^{ss}\,\mathbf{H}^{\mathsf{T}}\big{)}_{ij}$
(64) $\displaystyle\mathds{D}^{{ww}\,(r)}$ $\displaystyle\equiv~{}$
$\displaystyle
2^{r+1}\,\mathds{D}^{{ww}}~{},~{}~{}~{}~{}~{}~{}\,\mathds{D}^{{ww}}_{ij}~{}\equiv~{}\big{(}\mathbf{H}\,\mathds{D}^{ss}\,\mathbf{H}^{\mathsf{T}}\big{)}_{ij}.$
(65)

##  4 Generic wMERA formulation

The decomposition of the Hilbert space of a QFT into localized degrees of
freedom organized by a hierarchy of scales makes the connection between wMERA
and its discrete counterpart, MERA, more straightforward than the cMERA
approach. In this section, we outline the generic architecture of wMERA using
DWTs. To make the notation concrete, we assume a scalar field theory in (1+1)
spacetime dimensions, but this assumption can be easily generalized to
theories with more fields, either bosonic or fermionic. The generalization to
higher ($d$+1) spacetime dimensions can also be made by considering wavelet
bases in $L^{2}(\mathbb{R}^{d})$; this, however, is beyond the scope of this
paper.

In the previous section, we introduced the resolution index $r$ associated
with a lattice spacing $a_{0}/2^{r}$ for the scaling and wavelet bases. In
particular, under our unspecified choice for unit convention, the reference
lattice spacing for $r=0$ was set to $a_{0}=1$. Although the choice of $a_{0}$
is completely arbitrary, it is useful to identify it with some natural scale
of the theory. For the sake of specificity, we will adopt $a_{0}=1/m_{\phi}$,
where $m_{\phi}$ is the mass of the scalar field in our toy example. Given
this choice, we have that for an arbitrary resolution $r$, the “lattice site
index” $n$ of the scaling and wavelet functions represents a displacement from
the origin by $n$ integer steps of $1/(2^{r}\,m_{\phi})$. In addition, in
order to preserve the dimensionless inner product of the scaling functions, we
will also explicitly pull out a dimensionful normalization factor in their
definition:

$s_{n}^{r}\big{(}x\big{)}~{}\equiv~{}\sqrt{2^{r}\,m_{\phi}}\;s\big{(}2^{r}m_{\phi}\,x\,-\,n\big{)}\,.$
(66)

The discretized degrees of freedom of a scalar field $\hat{\Phi}(x)$ are then
obtained from smearing the field in our chosen wavelet basis:

$\displaystyle\hat{\phi}_{n}^{\,r}$ $\displaystyle~{}=~{}$ $\displaystyle\int
dx~{}\hat{\Phi}(x)\,s_{n}^{r}(x)$ (67a)
$\displaystyle\hat{\varphi}_{n}^{r^{\prime}}$ $\displaystyle~{}=~{}$
$\displaystyle\int
dx~{}\hat{\Phi}(x)\,w_{n}^{r^{\prime}}(x)~{}~{}~{}~{}~{}~{}(r^{\prime}\geq
r).$ (67b)

The corresponding discretized degrees of freedom of the canonical conjugate
momentum $\hat{\Pi}(x)$ are:

$\displaystyle\hat{\pi}_{n}^{r}$ $\displaystyle~{}=~{}$ $\displaystyle\int
dx~{}\hat{\Pi}(x)\,s_{n}^{r}(x)$ (68a)
$\displaystyle\hat{\varpi}_{n}^{r^{\prime}}$ $\displaystyle~{}=~{}$
$\displaystyle\int
dx~{}\hat{\Pi}(x)\,w_{n}^{r^{\prime}}(x)~{}~{}~{}~{}~{}~{}(r^{\prime}\geq r).$
(68b)

They obey the usual canonical commutation relations,

$\displaystyle\Big{[}\hat{\phi}_{n}^{\,r},\,\hat{\pi}_{m}^{r}\Big{]}~{}=~{}i\,\delta_{nm}~{}~{},~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\Big{[}\hat{\varphi}_{n}^{r^{\prime}},\,\hat{\varpi}_{m}^{r^{\prime\prime}}\Big{]}~{}=~{}i\,\delta_{r^{\prime}r^{\prime\prime}}\,\delta_{nm},$
(69)

all other commutators vanishing. With the definitions above, the DWT
representation of the scalar field $\hat{\Phi}(x)$ at a fixed resolution
($r$+1) can be written as:

$\displaystyle\hat{\Phi}^{r+1}(x)$ $\displaystyle~{}=$
$\displaystyle~{}\sum_{n}\;\hat{\phi}_{n}^{\,r+1}\,s_{n}^{r+1}(x)$ (70a)
$\displaystyle~{}=$
$\displaystyle~{}\sum_{n}\;\hat{\phi}_{n}^{\,r}\,s_{n}^{r}(x)~{}+~{}\hat{\varphi}_{n}^{r}\,w_{n}^{\,r}(x),$
(70b)

and similarly for the canonical conjugate momentum $\hat{\Pi}(x)$:

$\displaystyle\hat{\Pi}^{r+1}(x)$ $\displaystyle~{}=$
$\displaystyle~{}\sum_{n}\;\hat{\pi}_{n}^{r+1}\,s_{n}^{r+1}(x)$ (71a)
$\displaystyle~{}=$
$\displaystyle~{}\sum_{n}\;\hat{\pi}_{n}^{r}\,s_{n}^{r}(x)~{}+~{}\hat{\varpi}_{n}^{r}\,w_{n}^{\,r}(x).$
(71b)

Alternatively, the discretized degrees of freedom can be recast in terms of
creation and annihilation operators by defining:

$\displaystyle\hat{\phi}_{n}^{\,r}$ $\displaystyle~{}\equiv$
$\displaystyle~{}\frac{1}{\sqrt{2\Delta}}\Big{(}\hat{\sigma}_{n}^{r\,\dagger}+\hat{\sigma}_{n}^{r}\Big{)}~{},~{}~{}~{}~{}~{}~{}\hat{\pi}_{n}^{r}~{}\equiv~{}i\,\sqrt{\frac{\Delta}{2}}\;\Big{(}\hat{\sigma}_{n}^{r\,\dagger}-\hat{\sigma}_{n}^{r}\Big{)},$
(72a) $\displaystyle\hat{\varphi}_{n}^{r}$ $\displaystyle~{}\equiv$
$\displaystyle~{}\frac{1}{\sqrt{2\Delta}}\Big{(}\hat{\omega}_{n}^{r\,\dagger}+\hat{\omega}_{n}^{r}\Big{)}~{},~{}~{}~{}~{}~{}~{}\hat{\varpi}_{n}^{r}~{}\equiv~{}i\,\sqrt{\frac{\Delta}{2}}\;\Big{(}\hat{\omega}_{n}^{r\,\dagger}-\hat{\omega}_{n}^{r}\Big{)}.$
(72b)

Above, $\hat{\sigma}^{r}$ and $\hat{\sigma}^{r\,\dagger}$ are the creation and
annihilation operators for “collective excitations” at scales
$1/(2^{r}\,m_{\phi})$, whereas $\hat{\omega}^{r}$ and
$\hat{\omega}^{r\,\dagger}$ are the creation and annihilation operators for
the corresponding “fluctuating” wavelet modes. They obey the commutation
relations:

$\displaystyle\big{[}\hat{\sigma}_{n}^{r},\;\hat{\sigma}_{m}^{r\,\dagger}\big{]}~{}=~{}\delta_{nm}~{},~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\big{[}\hat{\omega}_{n}^{r},\;\hat{\omega}_{m}^{r^{\prime}\dagger}\big{]}~{}=~{}\delta_{rr^{\prime}}\,\delta_{nm}\,,$
(73)

with all other commutators vanishing. The width $\Delta$ appearing in (72a),
(72b) is arbitrary; it will be adjusted by the entangler when the wMERA
parameters are variationally optimized. An obvious initial choice is
$\Delta=m_{\phi}$.

Note that the discretized degrees of freedom also obey (inverse) wavelet
transform relations. Specifically, let $\\{\hat{\text{s}},\,\hat{\text{w}}\\}$
denote any one of the following pairs: $\\{\hat{\phi},\,\hat{\varphi}\\}$,
$\\{\hat{\pi},\,\hat{\varpi}\\}$, $\\{\hat{\sigma},\,\hat{\omega}\\}$, or
$\\{\hat{\sigma}^{\dagger},\,\hat{\omega}^{\dagger}\\}$. The (inverse) wavelet
transform relations for each of these pairs are given by:

WT
$\displaystyle\begin{cases}\,\hat{\text{s}}^{r}_{n}=~{}~{}\mathbf{L}_{nm}\,\hat{\text{s}}_{m}^{r+1},\\\
\,\hat{\text{w}}^{r}_{n}~{}=~{}\mathbf{H}_{nm}\,\hat{\text{s}}_{m}^{r+1},\end{cases}$
(74) IWT
$\displaystyle\begin{cases}\,\hat{\text{s}}^{r+1}_{n}~{}=~{}\mathbf{L}^{\mathsf{T}}_{nm}\,\hat{\text{s}}^{r}_{m}+\mathbf{H}^{\mathsf{T}}_{nm}\,\hat{\text{w}}^{r}_{m}.\end{cases}$
(75)

The D6-discretization of quantum fields also leads to discretized Hamiltonians
at different resolutions. Consider, for the sake of illustration, a generic
Hamiltonian for our toy example of a scalar field in (1+1)d:

$\hat{H}=\int
dx\,\Bigg{[}\,\frac{1}{2}\,\hat{\Pi}(x)^{2}~{}+~{}\frac{m_{\phi}^{2}}{2}\;\hat{\Phi}(x)^{2}~{}+~{}\frac{1}{2}\,\nabla\hat{\Phi}(x)\,\nabla\hat{\Phi}(x)~{}+~{}\mathcal{V}\big{[}\hat{\Phi}(x)\big{]}\Bigg{]}.$
(76)

Plugging expressions (70a) and (71a) into (76) and integrating over $x$, we
obtain the discretized system’s Hamiltonian at resolution $r$ in terms of
scaling modes:

$\hat{H}^{r}_{\phi}~{}=~{}\frac{1}{2}\,\sum_{i}\bigg{[}\,\big{(}\hat{\pi}_{i}^{r}\big{)}^{2}~{}+~{}m_{\phi}^{2}\big{(}\hat{\phi}_{i}^{\,r}\big{)}^{2}~{}+~{}(2^{r}m_{\phi})^{2}\,\sum_{j}\hat{\phi}_{i}^{\,r}\;\mathds{K}_{ij}^{{ss}}\;\hat{\phi}_{j}^{\,r}\bigg{]}~{}+~{}\mathcal{V}_{\phi}^{r}\big{[}\hat{\phi}^{\,r}\big{]},$
(77)

where the kinetic term coefficients $\mathds{K}^{{ss}}$ have been defined in
(52) and (53).

In (76), we left the interaction potential
$\mathcal{V}\big{[}\hat{\Phi}\big{]}$ unspecified. Assuming
$\mathcal{V}\big{[}\hat{\Phi}\big{]}=\lambda\,\hat{\Phi}^{n}$, the
corresponding discretized interaction potential in (77) will have the generic
form:

$\mathcal{V}_{\phi}^{r}\big{[}\hat{\phi}^{\,r}\big{]}~{}=~{}\lambda\\!\\!\\!\\!\sum_{~{}~{}i_{1},i_{2},\dots,i_{n}}\\!\\!\\!\\!I_{i_{1}i_{2}\dots
i_{n}}\;\hat{\phi}^{\,r}_{i_{1}}\;\hat{\phi}^{\,r}_{i_{2}}\,\dots\,\hat{\phi}^{\,r}_{i_{n}}\,,$
(78)

where

$I_{i_{1}i_{2}\dots i_{n}}~{}\equiv~{}\int
dx~{}s^{r}_{i_{1}}(x)\,s^{r}_{i_{2}}(x)\,\dots\,s^{r}_{i_{n}}(x).$ (79)

The discretized Hamiltonian can also be written in terms of creation and
annihilation operators of scaling modes by using (72a). Assuming
$\Delta=m_{\phi}$ for convenience, we have:

$\displaystyle\hat{H}^{r}_{\sigma}$ $\displaystyle~{}=~{}$
$\displaystyle\frac{m_{\phi}}{2}\,\sum_{i}\left(\hat{\sigma}^{r\,\dagger}_{i}\hat{\sigma}^{r}_{i}+\hat{\sigma}^{r}_{i}\,\hat{\sigma}^{r\,\dagger}_{i}\right)\,\;+\;~{}m_{\phi}\,2^{2r-2}\,\sum_{ij}\big{(}\hat{\sigma}^{r\,\dagger}_{i}+\hat{\sigma}^{r}_{i}\big{)}\,\mathds{K}_{ij}^{{ss}}\,\big{(}\hat{\sigma}^{r\,\dagger}_{j}+\hat{\sigma}^{r}_{j}\big{)}$
(80)
$\displaystyle+~{}\mathcal{V}_{\sigma}^{r}\big{[}\hat{\sigma}^{r(\dagger)}\big{]}.$

Finally, through the use of IWTs defined in (75)—or, alternatively, by using
the field decompositions in (70b) and (71b) combined with (72)—the discretized
Hamiltonian can also be cast in terms of scaling and wavelet modes.
Concretely,

$\displaystyle\hat{H}_{\sigma\omega}^{(r+1)}$ $\displaystyle~{}=~{}~{}$
$\displaystyle\hat{H}_{\sigma}^{r}~{}+~{}\frac{m_{\phi}}{2}\,\sum_{i}\left(\hat{\omega}^{r\,\dagger}_{i}\hat{\omega}^{r}_{i}+\hat{\omega}^{r}_{i}\,\hat{\omega}^{r\,\dagger}_{i}\right)$
(81)
$\displaystyle+~{}m_{\phi}\,2^{2r-2}\,\sum_{ij}\big{(}\hat{\omega}^{r\,\dagger}_{i}+\hat{\omega}^{r}_{i}\big{)}\,\mathds{K}_{ij}^{{ww}}\,\big{(}\hat{\omega}^{r\,\dagger}_{j}+\hat{\omega}^{r}_{j}\big{)}$
$\displaystyle+~{}m_{\phi}\,2^{2r-1}\,\sum_{ij}\big{(}\hat{\sigma}^{r\,\dagger}_{i}+\hat{\sigma}^{r}_{i}\big{)}\,\mathds{K}_{ij}^{{sw}}\,\big{(}\hat{\omega}^{r\,\dagger}_{j}+\hat{\omega}^{r}_{j}\big{)}$
$\displaystyle+~{}\mathcal{V}_{\sigma\omega}^{(r+1)}\big{[}\hat{\sigma}^{r(\dagger)}\\!,\,\hat{\omega}^{r(\dagger)}\big{]},$

where the kinetic term coefficients $\mathds{K}^{{sw}}$ and
$\mathds{K}^{{ww}}$ have been defined in (55) and (56), respectively.

The wavelet decomposition of the field degrees of freedom also implies a
decomposition of the Hilbert space of this theory into orthogonal subspaces at
different scales:

$\mathcal{H}~{}=~{}\mathcal{H}_{r}~{}\oplus~{}\mathfrak{h}_{r}~{}\oplus~{}\mathfrak{h}_{r+1}~{}\oplus...\oplus~{}\mathfrak{h}_{r+j}~{}\oplus...\,.$
(82)

Above, $\mathcal{H}_{r}$ is the Hilbert subspace containing the “long-range”
scaling modes describing collective excitations of the scalar field over
scales $1/(2^{r}m_{\phi})$. It is identified with the Fock space generated by
creation operators $\hat{\sigma}^{r\,\dagger}$ acting on the product state
defined by:

$\big{|}{\mathbf{P}_{\\!\sigma}\big{\rangle}}_{\\!r}~{}\equiv~{}\prod_{n}\,\big{|}0_{\sigma_{n}}\big{\rangle}_{\\!r}~{},~{}~{}~{}\text{with}~{}~{}~{}~{}~{}~{}\hat{\sigma}_{n}^{r}\,\big{|}0_{\sigma_{n}}\big{\rangle}_{\\!r}~{}=~{}0\,.$
(83)

Likewise, the finer resolution Hilbert subspaces $\mathfrak{h}_{r^{\prime}}$
($r^{\prime}\geq r$) contain modes describing short-range quantum fluctuations
of the scalar field at scales $1/(2^{r^{\prime}+1}m_{\phi})$. Each of these
subspaces is identified with the Fock space generated by the creation
operators $\hat{\omega}^{r^{\prime}\dagger}$ acting on the product state
defined by:

$\big{|}{\mathbf{P}_{\\!\omega}\big{\rangle}}_{\\!r^{\prime}}~{}=~{}\prod_{n}\,\big{|}0_{\omega_{n}}\big{\rangle}_{\\!r^{\prime}}~{},~{}~{}~{}\text{with}~{}~{}~{}~{}~{}~{}\hat{\omega}_{n}^{r^{\prime}}\,\big{|}0_{\omega_{n}}\big{\rangle}_{\\!r^{\prime}}~{}=~{}0\,.$
(84)

With the Hilbert space decomposition described in (82)-(84), we can sketch the
main ingredients of wMERA:

(1)

At very low resolutions $r_{\text{\tiny IR}}\ll 0$, the discretized degrees of
freedom describe fluctuations at lengths much longer than the correlation
length of the system (that is, $L\gg 1/m_{\phi}$ in our toy example). At such
scales, the nonlocal terms in the Hamiltonian—including the kinetic term—are
suppressed, and the ground state exhibits exponentially negligible spatial
entanglement. Hence, the far infrared ground state is well-approximated by a
product state ansatz666If the (non-quadratic) Hamiltonian interactions are
non-negligible in the far IR, the ground state might not coincide with the
state in (83). Nevertheless, it should be well-approximated by a product state
which can be related to the state in (83) via an ultra-local entangler.,

$\big{|}{\mathbf{0}\big{\rangle}}_{\\!r_{{}_{\text{IR}}}}~{}\simeq~{}\big{|}{\mathbf{P}_{\\!\sigma}\big{\rangle}}_{\\!r_{{}_{\text{IR}}}}.$
(85)

At resolutions close to the correlation length of the system, (that is, $r\sim
0$ and $L\sim 1/m_{\phi}$ in our toy example), entanglement is still
relatively localized (i.e., it is only non-negligible between neighboring
sites, and decays exponentially for site separations $|n-m|\gg 1$). Hence, any
resolution close to $r=0$ is a good starting point for building the wMERA.
Choosing, for concreteness, $r=0$ as the initial resolution, the ground state
can be described by the following ansatz state in the Hilbert space
$\mathcal{H}_{0}$:

$\big{|}{\mathbf{0_{\sigma}}\big{\rangle}}_{\\!0}~{}\simeq~{}\hat{U}_{0}\big{|}{\mathbf{P}_{\\!\sigma}\big{\rangle}}_{\\!0}\,,$
(86)

where $\hat{U}_{0}$ is a quasi-local unitary entangler ansatz. The (_a priori_
unkown) parameters of $\hat{U}_{0}$ are obtained variationally by minimizing
the expectation value of the discretized Hamiltonian at $r=0$:

${\big{\langle}}\mathbf{0_{\sigma}}\big{|}\,\hat{H}_{\sigma}^{0}\,\big{|}{\mathbf{0_{\sigma}}}\big{\rangle}_{\\!0}~{}=~{}{\big{\langle}}{\mathbf{P}_{\\!\sigma}}\big{|}\,\hat{U}_{0}^{\dagger}\,\hat{H}_{\sigma}^{0}\,\hat{U}_{0}\,\big{|}{{\mathbf{P}_{\\!\sigma}}}\big{\rangle}_{\\!0}.$
(87)

This constitutes the first wMERA layer.

(2)

To build the second layer, the first step is to enlarge the Hilbert space to
include the finer resolution wavelet subspace $\mathfrak{h}_{0}$ via the fine-
graining mapping ${W}_{1}$:

${W}_{1}:~{}~{}\mathcal{H}_{0}~{}~{}\mapsto~{}~{}\mathcal{H}_{0}\,\oplus\,\mathfrak{h}_{0}\;\cong\;\mathcal{H}_{1}\,.$
(88)

Accordingly, the ansatz is augmented by wavelet states initialized in a
product state with the rest of the system. By abuse of notation, we write:

$\hat{W}_{1}\,\big{|}{\mathbf{0_{\sigma}}\big{\rangle}}_{\\!0}~{}\equiv~{}\big{|}{\mathbf{0_{\sigma}}\big{\rangle}}_{\\!0}\otimes\big{|}{\mathbf{P}_{\\!\omega}\big{\rangle}}_{\\!0}\,.$
(89)

(3)

The second step is to entangle the newly introduced wavelet states with the
original, lower-resolution scaling states. This is parameterized by a quasi-
local entangler $\hat{V}_{1}$:

$\big{|}{\mathbf{0_{\sigma\omega}}\big{\rangle}}_{\\!1}~{}=~{}\hat{V}_{1}\,\Big{[}\big{|}{\mathbf{0_{\sigma}}\big{\rangle}}_{\\!0}\otimes\big{|}{\mathbf{P}_{\\!\omega}\big{\rangle}}_{\\!0}\Big{]}.$
(90)

Similarly to the first layer, the parameters of the entangler ansatz
$\hat{V}_{1}$ are variationally optimized by minimizing the expectation value
of the discretized Hamiltonian at $r$ $=$ $1$:

${\big{\langle}}\mathbf{0_{\sigma\omega}}\big{|}\hat{H}^{1}_{\sigma\omega}\big{|}{\mathbf{0_{\sigma\omega}}}\big{\rangle}_{\\!1}~{}=~{}{\big{\langle}}{\mathbf{P}_{\\!\sigma\omega}}\big{|}\,\big{(}\hat{V}_{1}\,\hat{U}_{0}\big{)}^{\dagger}\hat{H}^{1}_{\sigma\omega}\big{(}\hat{V}_{1}\,\hat{U}_{0}\big{)}\,\big{|}{{\mathbf{P}_{\\!\sigma\omega}}}\big{\rangle}_{\\!1},$
(91)

where, above, the $r=1$ Hamiltonian $\hat{H}^{1}_{\sigma\omega}$ is expressed
in terms of scaling and wavelet modes as in (81), and

$\big{|}{\mathbf{P}_{\\!\sigma\omega}\big{\rangle}}_{\\!1}~{}\equiv~{}\big{|}{\mathbf{P}_{\\!\sigma}}\big{\rangle}_{\\!0}\otimes\big{|}{\mathbf{P}_{\\!\omega}}\big{\rangle}_{\\!0}.$
(92)

(4)

The state in (90) is the ansatz for the ground state at resolution $r=1$
defined in the Fock space of $r=0$ scaling and wavelet states. As we will
justify shortly, at this point it is convenient to perform an inverse wavelet
transform $\big{(}\hat{T}_{1}\big{)}^{-1}$ to re-express this ansatz as a
state in the Fock space of $r=1$ scaling states:

$\big{|}{\mathbf{0_{\sigma}}\big{\rangle}}_{\\!1}~{}=~{}\big{(}\hat{T}_{1}\big{)}^{\\!-1}\,\big{|}{\mathbf{0_{\sigma\omega}}\big{\rangle}}_{\\!1}\,.$
(93)

Note that the variational optimization of the entangler $\hat{V}_{1}$ can be
alternatively performed _after_ the IWT step is applied. Noting that the $r$ =
1 Hamiltonians $\hat{H}^{1}_{\sigma}$ and $\hat{H}^{1}_{\sigma\omega}$
(generically defined in (80) and (81)) are related by

$\hat{H}^{1}_{\sigma}=\big{(}\hat{T}_{1}\big{)}^{\\!-1}\hat{H}^{1}_{\sigma\omega}\;\hat{T}_{1},$
(94)

and defining

$\hat{U}_{1}=\big{(}\hat{T}_{1}\big{)}^{\\!-1}\big{(}\hat{V}_{1}\,\hat{U}_{0}\big{)}\,\hat{T}_{1},$
(95)

we can alternatively perform the variational optimization of the $r$ = 1
entangler parameters by minimizing the following expectation value:

${\big{\langle}}\mathbf{0_{\sigma}}\big{|}\,\hat{H}_{\sigma}^{1}\,\big{|}{\mathbf{0_{\sigma}}}\big{\rangle}_{\\!1}~{}=~{}{\big{\langle}}{\mathbf{P}_{\\!\sigma}}\big{|}\,\hat{U}_{1}^{\dagger}\,\hat{H}_{\sigma}^{1}\,\hat{U}_{1}\,\big{|}{{\mathbf{P}_{\\!\sigma}}}\big{\rangle}_{\\!1}.$
(96)

(5)

Finally, these operations can be performed iteratively to construct the layers
of the wMERA until one reaches either a fixed point, or the desired UV
resolution $r_{\text{\tiny UV}}$:

$\boxed{~{}\big{|}{\mathbf{0_{\sigma}}\big{\rangle}}_{\\!r_{{}_{\text{\tiny
UV}}}}~{}=~{}\left[\\!\\!\prod_{~{}r~{}=~{}r_{{}_{\text{\tiny
IR}}}}^{r_{{}_{\text{\tiny
UV}}}}\big{(}\hat{T}_{r}\big{)}^{\\!-1}\hat{V}_{r}\,\hat{W}_{r}\right]\\!\\!\big{|}{\mathbf{P}_{\\!\sigma}\big{\rangle}}_{\\!r_{{}_{\text{\tiny
IR}}}}\,}\,.$ (97)

This state describes the vacuum of the theory at scale $r_{\text{\tiny UV}}$
provided that the parameters of the transformations $\hat{V}_{r}$ are, at each
layer, chosen so as to minimize the expectation value of the system’s
Hamiltonian $\hat{H}^{r}$. If these conditions are not sufficient to uniquely
determine the ansatz parameters, one can also demand that expectation values
of other physical observables are preserved under RG flow, such as $n$-point
correlation functions with $n>2$.

We expect that the variational optimization of the wMERA parameters in (97)
can be efficiently implemented numerically with unitary ansatzes $\hat{V}_{r}$
that are _quasi-local_ and _sparse_. If this is indeed the case, then
computational techniques for the variational optimization of MERA tensor
networks developed for discrete systems should be adaptable to wMERA for
studying QFTs.

In the next section, we will show that quasi-local and sparse wMERA unitaries
are indeed sufficient to describe ground state entanglement in two simple
examples: a free scalar field and a free Dirac fermion in (1+1) spacetime
dimensions.

We end this section by justifying the need for inverse wavelet transforms at
every resolution layer in the wMERA formulation. Note that if these
transformations were not performed, the description of the system would incur
two problematic features: (1) a non-uniform lattice structure, and (2)
discretized Hamiltonians $\hat{H}^{r}$ becoming increasingly nonlocal as the
number of wMERA layers is increased. To illustrate this situation, consider a
resolution layer $r=r_{0}+n$, where $n\gg 1$ and $r_{0}$ is the initial wMERA
layer. In the absence of IWTs, the degrees of freedom at this layer would
consist of the scaling modes, $\phi^{(r_{0})}$, and $n+1$ types of wavelet
modes:
$\varphi^{(r_{0})},\,\varphi^{(r_{0}+1)},\,\dots,\,\varphi^{(r_{0}+n)}$. These
$n+2$ different types of modes would have different lattice densities, with
$\phi^{(r_{0})}$ and $\varphi^{(r_{0})}$ having the sparsest lattices (with
lattice spacing $a_{0}$), and $\varphi^{(r_{0}+n)}$ having the densest lattice
(with lattice spacing $a_{0}/2^{n}$). In addition, any given scaling mode
$\phi_{j}^{(r_{0})}$ would couple to a total of
$\big{(}5n^{2}+23n+36\big{)}/2$ other modes, namely,

$\begin{cases}&\phi_{j-4}^{(r_{0})}\;\dots\;\phi_{j+4}^{(r_{0})}\\\
&\varphi_{j-4}^{(r_{0})}\;\dots\;\varphi_{j+4}^{(r_{0})}\\\
&\varphi_{2j-4}^{(r_{0}+1)}\;\dots\;\varphi_{2j+9}^{(r_{0}+1)}\\\
&~{}~{}~{}\vdots\\\
&\varphi_{2^{n}j-4}^{(r_{0}+n)}\;\dots\;\varphi_{2^{n}j+5n+4}^{(r_{0}+n)}.\end{cases}$

In other words, the nonlocality of the Hamiltonians $\hat{H}^{r_{0}+n}$ would
grow quadratically with the number of wMERA layers. This situation is depicted
in figure 6. It is easy to see that not only this lattice structure becomes
quickly impractical, but the ansatz would require unitary entanglers that more
and more nonlocal at each layer. The IWTs address these difficulties by
redefining the degrees of freedom when moving from one layer to the next so
that the uniformity of the lattice structure is preserved, and the quasi-
locality of the discretized Hamiltonians $\hat{H}^{r}$ is restored at each
layer.

Figure 6: Depiction of the lattice describing the system at a deep UV layer if
inverse wavelet transforms were not incorporated in the wMERA architecture.
The degrees of freedom of the first (bottom layer) consist of scaling modes
(depicted in yellow) and wavelet modes (depicted in dark blue). The next
layers would iteratively introduce shorter- and shorter-distance wavelet modes
(depicted in dark blue) with successively denser lattices. In this lattice
architecture, a scaling mode in the first layer would couple to
$\mathcal{O}(n^{2})$ wavelet modes, where $n$ is the number of wMERA layers
that support this particular resolution.

##  5 Concrete wMERA Examples

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

###  5.1 Free theory of a scalar field in (1+1)d

Consider a free scalar field theory in (1+1)d. Its Hamiltonian is given by

$\hat{H}=\frac{1}{2}\int
dx\,\Bigg{[}\,\hat{\Pi}(x)^{2}~{}+~{}m_{\phi}^{2}\,\hat{\Phi}(x)^{2}~{}+~{}\nabla\hat{\Phi}(x)\,\nabla\hat{\Phi}(x)\Bigg{]},$
(98)

where the scalar field $\hat{\Phi}(x)$ its canonical conjugate momentum
$\hat{\Pi}(x)$ obey the equal-time commutation relations

$\Big{[}\hat{\Phi}(x),\,\hat{\Pi}(y)\Big{]}~{}=~{}i\,\delta(x-y);$ (99)

all other equal-time commutators vanishing.

Since free theories only have Gaussian correlations, it should come as no
surprise that their entanglers are Gaussian777For an excellent reference on
multimode Gaussian entanglers, see PhysRevA.41.4625 . (i.e., the entanglers’
exponents are quadratic functions of $\hat{\Phi}$ and $\hat{\Pi}$).
Specifically, writing the field and its canonical conjugate momentum in terms
of _position-space_ creation and annihilation operators:

$\hat{\Phi}(x)~{}\equiv~{}\frac{1}{\sqrt{2m_{\phi}}}\Big{(}\hat{a}(x)^{\dagger}+\hat{a}(x)\Big{)}~{},~{}~{}~{}~{}~{}~{}\hat{\Pi}(x)~{}\equiv~{}i\,\sqrt{\frac{m_{\phi}}{2}}\;\Big{(}\hat{a}(x)^{\dagger}-\hat{a}(x)\Big{)},$
(100)

we would like to find the entangler $\hat{U}$ that relates the product state
$\big{|}{\mathbf{P}\big{\rangle}}$ to the ground state
$\big{|}{\mathbf{0}\big{\rangle}}$:

$\big{|}{\mathbf{0}\big{\rangle}}~{}=~{}\hat{U}\big{|}{\mathbf{P}\big{\rangle}},$
(101)

where the product state $\big{|}{\mathbf{P}\big{\rangle}}$ has, by definition,
no entanglement in position space and can be thought of as the generalization
of the state in (83) to the continuum:

$\hat{a}(x)\,\big{|}\mathbf{P}\big{\rangle}~{}=~{}0~{}~{}~{}\forall~{}~{}~{}x.$
(102)

In this case, the entangler ansatz is a generator of Bogoliubov
transformations (a.k.a. multimode squeezing), and has the generic form:

$\hat{U}~{}=~{}\exp\left[-\frac{i}{2}\int
dx\,dy~{}\zeta(x-y)~{}\Big{(}\hat{\Phi}(x)\hat{\Pi}(y)+\hat{\Pi}(x)\hat{\Phi}(y)\Big{)}\,\right],$
(103)

or, alternatively in terms of position-space creation and annihilation
operators,

$\boxed{~{}\hat{U}~{}=~{}\exp\left[~{}\frac{1}{2}\int
dx\,dy~{}\zeta(x-y)~{}\Big{(}\hat{a}^{\dagger}(x)\hat{a}^{\dagger}(y)-\hat{a}(x)\hat{a}(y)\Big{)}\,\right]\,}.$
(104)

The entangler exponent $\zeta(x-y)$ defines the entangling parameter that
needs to be variationally optimized. The easiest way to proceed is to go to
Fourier space by defining:

$\displaystyle\hat{\Phi}_{k}$ $\displaystyle\;\equiv\,\int
dx~{}e^{-ikx}\,\hat{\Phi}(x)~{},~{}~{}~{}~{}~{}~{}~{}$
$\displaystyle\hat{\Pi}_{k}\;\equiv\,\int dx~{}e^{-ikx}\,\hat{\Pi}(x)$ (105a)
$\displaystyle\,\hat{a}_{k}$ $\displaystyle\;\equiv\,\int
dx~{}e^{-ikx}\,\hat{a}(x)~{},~{}~{}~{}~{}~{}~{}~{}$
$\displaystyle\;\hat{a}^{\dagger}_{k}\;\equiv\,\int
dx~{}e^{ikx}\,\hat{a}^{\dagger}(x)$ (105b)

and

$\zeta_{k}~{}\equiv~{}\int dx~{}e^{-ikx}\,\zeta(x).$ (106)

Note that $\hat{a}^{\dagger}_{k}$ and $\hat{a}_{k}$ in (105b) are _not_ the
creation and annihilation operators of momentum-mode Hamiltonian eigenstates.
In particular, they do _not_ annihilate the ground state, but, because of
(102), they do annihilate the initial product state:

$\hat{a}_{k}\,\big{|}\mathbf{P}\big{\rangle}~{}=~{}0~{}~{}~{}\forall~{}~{}~{}k.$
(107)

In Fourier space, the entangler and the Hamiltonian are given by

$\hat{U}~{}=~{}\exp\left[\,\frac{1}{2}\,\int\frac{dk}{2\pi}~{}\zeta_{k}~{}\Big{(}\hat{a}^{\dagger}_{k}\,\hat{a}^{\dagger}_{-k}-\hat{a}_{k}\,\hat{a}_{-k}\Big{)}\,\right]$
(108)

and

$\hat{H}~{}=~{}\frac{1}{4m_{\phi}}\int\frac{dk}{2\pi}~{}~{}\Bigg{[}\big{(}\omega_{k}^{2}+m_{\phi}^{2}\big{)}\Big{(}{\hat{a}_{k}}^{\dagger}\,\hat{a}_{k}+\hat{a}_{k}\,\hat{a}^{\dagger}_{k}\Big{)}+k^{2}\,\Big{(}\hat{a}^{\dagger}_{k}\,\hat{a}^{\dagger}_{-k}+\hat{a}_{k}\,\hat{a}_{-k}\Big{)}\Bigg{]},$
(109)

where $\omega_{k}^{2}$ has the usual definition $\omega_{k}^{2}\equiv
k^{2}+m_{\phi}^{2}$. Recalling the commutation relations
$\big{[}\hat{a}_{k},\,\hat{a}_{q}^{\dagger}\big{]}=2\pi\,\delta(k-q)$, one can
easily derive the action of the entangler on the creation and annihilation
operators,

$\displaystyle\hat{U}^{\dagger\,}\hat{a}_{k}\,\hat{U}$ $\displaystyle~{}=~{}$
$\displaystyle\text{cosh}(\zeta_{k})\,\hat{a}_{k}~{}+~{}\text{sinh}(\zeta_{k})\,\hat{a}_{-k}^{\dagger}$
(110a) $\displaystyle\hat{U}^{\dagger\,}\hat{a}_{k}^{\dagger}\,\hat{U}$
$\displaystyle~{}=~{}$
$\displaystyle\text{cosh}(\zeta_{k})\,\hat{a}_{k}^{\dagger}~{}+~{}\text{sinh}(\zeta_{k})\,\hat{a}_{-k},$
(110b)

and the entangler-transformed Hamiltonian,

$\displaystyle\hat{U}^{\dagger}\hat{H}\,\hat{U}~{}=~{}\frac{1}{4m_{\phi}}\int\frac{dk}{2\pi}$
$\displaystyle\Bigg{[}\Big{(}\big{(}\omega_{k}^{2}+m_{\phi}^{2}\big{)}\,\text{cosh}(2\zeta_{k})+k^{2}\,\text{sinh}(2\zeta_{k})\Big{)}\,\Big{(}{\hat{a}_{k}}^{\dagger}\,\hat{a}_{k}+\hat{a}_{k}\,\hat{a}^{\dagger}_{k}\Big{)}$
$\displaystyle~{}~{}+~{}\Big{(}\big{(}\omega_{k}^{2}+m_{\phi}^{2}\big{)}\,\text{sinh}(2\zeta_{k})+k^{2}\,\text{cosh}(2\zeta_{k})\Big{)}\,\Big{(}\hat{a}^{\dagger}_{k}\,\hat{a}^{\dagger}_{-k}+\hat{a}_{k}\,\hat{a}_{-k}\Big{)}\Bigg{]}.$

Obtaining the vacuum expectation value of the Hamiltonian is now
straightforward. Using (5.1) and (107), we have:

$\displaystyle\big{\langle}\mathbf{0}\big{|}\hat{H}\big{|}\mathbf{0}\big{\rangle}~{}$
$\displaystyle=$
$\displaystyle~{}{\big{\langle}}\mathbf{P}\big{|}\hat{U}^{\dagger}\hat{H}\,\hat{U}\big{|}\mathbf{P}\big{\rangle}$
(112) $\displaystyle=$ $\displaystyle~{}\frac{1}{4m_{\phi}}\int
dk\,\,\delta(k-k)\,\Big{(}\big{(}\omega_{k}^{2}+m_{\phi}^{2}\big{)}\,\text{cosh}(2\zeta_{k})+k^{2}\,\text{sinh}(2\zeta_{k})\Big{)}.$

The final step is to minimize the vacuum expectation value of the Hamiltonian
with respect to entangler exponent $\zeta_{k}$:

$\displaystyle\frac{\partial}{\partial\zeta_{k}}\big{\langle}\mathbf{0}\big{|}\hat{H}\big{|}\mathbf{0}\big{\rangle}=0~{}~{}$
$\displaystyle\Rightarrow$
$\displaystyle~{}~{}\big{(}\omega_{k}^{2}+m_{\phi}^{2}\big{)}\,\text{sinh}(2\zeta_{k})+k^{2}\,\text{cosh}(2\zeta_{k})~{}=~{}0$
(113) $\displaystyle\Rightarrow$
$\displaystyle~{}~{}~{}\boxed{~{}\zeta_{k}\,=\,\frac{1}{2}\,\text{log}\Bigg{(}\frac{m_{\phi}}{\omega_{k}}\Bigg{)}\,}\,.$

We can now use (106) and (113) to obtain the _exact_ entangler exponent in
position space,

$\boxed{~{}\zeta(x-y)~{}=~{}\frac{\,e^{-m_{\phi}|x-y|}\,}{4|x-y|}\;}\,.$ (114)

It is also informative to look at the entangler-transformed Hamiltonian using
the exact entangler in (113) and (114):

$\displaystyle\hat{U}^{\dagger}\hat{H}\,\hat{U}~{}~{}$ $\displaystyle=$
$\displaystyle~{}\int\frac{dk}{2\pi}\,\omega_{k}\,{\hat{a}_{k}}^{\dagger}\,\hat{a}_{k}~{}~{}+~{}~{}\text{(zero-
point energy)}$ (115a) $\displaystyle=$
$\displaystyle~{}\int\\!dx\,dy~{}\,\Omega(x-y)\,\hat{a}^{\dagger}(x)\hat{a}(y)~{}~{}+~{}~{}\text{(zero-
point energy)},$ (115b)

where, above,

$\Omega(x-y)~{}\equiv~{}-\frac{\,m_{\phi}^{2}\,}{\pi}\,\frac{K_{1}\big{(}m_{\phi}|x-y|\big{)}}{m_{\phi}|x-y|}~{}~{}~{}~{}\xrightarrow[~{}m_{\phi}\to\,0~{}]{}~{}~{}~{}~{}-\frac{1}{\pi|x-y|^{2}}\,.$
(116)

While the use of FTs in this exercise was convenient to obtain the exact
entangler, this trick unfortunately only works well in the case of free
theories for which the entangler is Gaussian (as we argued extensively in
sections 1 and 2). For interacting theories, we expect that working directly
in position space through the use of DWTs will make calculations much more
tractable. For the remainder of this subsection, we show how this actually
works by solving this theory using wMERA.

We work with the D6-discretized degrees of freedom defined in (67), (68), and
(72), whose dynamics are governed by the discretized Hamiltonians (80) and
(81) with vanishing interaction potential (i.e., $\mathcal{V}=0$). To recall
our notation, at any given resolution $r$ the product state
$\big{|}\mathbf{P}_{\\!\sigma}\big{\rangle}_{\\!r}$ is related to the ground
state by a unitary entangler $\hat{U}_{r}$:

$\big{|}\mathbf{0}_{\sigma}\big{\rangle}_{\\!r}~{}=~{}\hat{U}_{r}\,\big{|}\mathbf{P}_{\\!\sigma}\big{\rangle}_{\\!r}.$
(117)

In our free theory example, the ansatz for $\hat{U}_{r}$ is a generator of
Bogoliubov transformations between scaling modes:

$\hat{U}_{r}~{}=~{}\exp\left[\,\frac{1}{2}\,\sum_{i,j}\,\zeta_{ij}^{(r)}\;\Big{(}\hat{\sigma}_{i}^{\,r\,\dagger}\hat{\sigma}_{j}^{\,r\,\dagger}\,-\,\hat{\sigma}_{i}^{\,r}\hat{\sigma}_{j}^{\,r}\Big{)}\,\right]\\!.$
(118)

The entangling exponent ${\zeta}^{(r)}$ is a Hermitian matrix whose entries
satisfy the following relation

$\zeta_{ij}^{(r)}~{}=~{}\zeta_{|i-j|}^{(r)}$ (119)

due to discrete translational invariance. The entangler $\hat{U}_{r}$
transforms scaling modes in the following way:

$\displaystyle\hat{U}_{r}^{\dagger}\,\hat{\sigma}_{i}^{\,r}\,\hat{U}_{r}$
$\displaystyle~{}=~{}$
$\displaystyle\big{(}\text{Cosh}\zeta^{(r)}\big{)}_{ij}\,\hat{\sigma}_{j}^{\,r}~{}~{}+~{}\big{(}\text{Sinh}\zeta^{(r)}\big{)}_{ij}\,\hat{\sigma}_{j}^{\,r\,\dagger}$
(120a)
$\displaystyle\hat{U}_{r}^{\dagger}\,\hat{\sigma}_{i}^{\,r\,\dagger}\,\hat{U}_{r}$
$\displaystyle~{}=~{}$
$\displaystyle\big{(}\text{Cosh}\zeta^{(r)}\big{)}_{ij}\,\hat{\sigma}_{j}^{\,r\,\dagger}~{}\\!+~{}\big{(}\text{Sinh}\zeta^{(r)}\big{)}_{ij}\,\hat{\sigma}_{j}^{\,r}\,.$
(120b)

The other type of entangler discussed in section 4 was the unitary operator
$\hat{V}_{r}$ that entangled scaling and wavelet modes at the wMERA layer with
resolution $r$. While $\hat{V}_{r}$ is still a Gaussian entangler, it has a
more general parametrization than $\hat{U}_{r}$. Concretely,

$\displaystyle\hat{V}_{r}$ $\displaystyle=$
$\displaystyle\exp\Bigg{\\{}\sum_{i,j}\Bigg{[}\,\nu^{ww(r-1)}_{ij}\Big{(}\hat{\omega}_{i}^{\,r-1\,\dagger}\hat{\omega}_{j}^{\,r-1\,\dagger}\\!-\hat{\omega}_{i}^{\,r-1}\hat{\omega}_{j}^{\,r-1}\Big{)}\\!+\nu^{sw(r-1)}_{ij}\Big{(}\hat{\sigma}_{i}^{\,r-1\,\dagger}\hat{\omega}_{j}^{\,r-1\,\dagger}\\!-\hat{\sigma}_{i}^{\,r-1}\hat{\omega}_{j}^{\,r-1}\Big{)}$
(121)
$\displaystyle~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}+~{}\rho^{sw(r-1)}_{ij}\Big{(}\hat{\sigma}_{i}^{\,r-1\,\dagger}\hat{\omega}_{j}^{\,r-1}-\hat{\omega}_{j}^{\,r-1\,\dagger}\hat{\sigma}_{i}^{\,r-1}\Big{)}\,+\,\dots~{}\Bigg{]}\,\Bigg{\\}}.$

Above, the ellipsis denotes higher order terms in a BCH expansion which are
proportional to $\hat{\sigma}_{i}^{\dagger}\hat{\sigma}_{j}$,
$\hat{\omega}_{i}^{\dagger}\hat{\omega}_{j}$, and
$\hat{\sigma}_{i}^{\dagger}\hat{\sigma}_{j}^{\dagger}-\hat{\sigma}_{i}\hat{\sigma}_{j}$.

As mentioned in section 4, at each wMERA layer one can either variationally
optimize the $\hat{V}_{r}$ ansatz parameters in (121) by minimizing

${\big{\langle}}\mathbf{0_{\sigma\omega}}\big{|}\hat{H}^{r}_{\sigma\omega}\big{|}{\mathbf{0_{\sigma\omega}}}\big{\rangle}_{\\!r}~{}=~{}\prescript{}{(r-1)\\!}{\Big{[}}\big{\langle}{\mathbf{P}_{\\!\omega}\big{|}}\otimes{\big{\langle}}\mathbf{0_{\sigma}}\big{|}\Big{]}\,\hat{V}_{r}^{\dagger}\,\hat{H}^{r}_{\sigma\omega}\,\hat{V}_{r}\,\Big{[}\big{|}{\mathbf{0_{\sigma}}\big{\rangle}}\otimes\big{|}{\mathbf{P}_{\\!\omega}\big{\rangle}}\Big{]}_{\\!(r-1)},$
(122)

or one can perform the optimization _after_ the IWT reorganizes the scaling
and wavelet modes $\hat{\sigma}^{\,(r-1)(\dagger)}$ and
$\hat{\omega}^{\,(r-1)(\dagger)}$ into $\hat{\sigma}^{\,r(\dagger)}$ scaling
modes. While this choice is completely up to the developer of the optimization
algorithm, this author found the second option easier to implement
numerically, for two main reasons. First, there is a significant amount of
redundancy in the $\hat{V}_{r}$ ansatz parameters $\nu_{ij}^{ww(r)}$,
$\nu_{ij}^{sw(r)}$, and $\rho_{ij}^{sw(r)}$ (they secretly depend on a smaller
number of parameters, as it should became clear shortly). Second, while we
have not shown explicitly how the $\hat{V}_{r}$ ansatz in (121) transforms the
scaling and wavelet modes $\hat{\sigma}_{i}^{\,r(\dagger)}$ and
$\hat{\omega}_{i}^{\,r(\dagger)}$, one can easily anticipate that these
transformations will have a more complicated form compared to (120).

Hence, we choose to proceed with the wMERA optimization by instead minimizing
the parameters of each layer _after_ the IWT is performed. That is, we will
optimize

${\big{\langle}}\mathbf{0_{\sigma}}\big{|}\,\hat{H}_{\sigma}^{r}\,\big{|}{\mathbf{0_{\sigma}}}\big{\rangle}_{\\!r}~{}=~{}{\big{\langle}}{\mathbf{P}_{\\!\sigma}}\big{|}\,\hat{U}_{r}^{\dagger}\,\hat{H}_{\sigma}^{r}\,\hat{U}_{r}\,\big{|}{{\mathbf{P}_{\\!\sigma}}}\big{\rangle}_{\\!r}.$
(123)

Note, however, that this does not mean we intend minimize all
$\zeta_{ij}^{(r)}$ parameters of the ansatz in (118) at once, since this would
not take any input from previous wMERA layers, and would lead to challenges
with numerical efficiency and accuracy. Instead, the crucial point to note is
that the “long range” entanglement encoded in $\zeta_{ij}^{(r)}$ has already
been optimized in the previous wMERA layers. The way to take advantage of that
is through the following: using WTs, we can recast $\hat{U}_{r}$ in terms of
scaling and wavelet degrees of freedom888The first equality in (5.1) is
essentially a generalization of the $r=1$ relation in (95).,

$\displaystyle\hat{T}_{r}\,\hat{U}_{r\,}\big{(}\hat{T}_{r}\big{)}^{\\!-1}\,$
$\displaystyle=$ $\displaystyle~{}\hat{V}_{r}\,\hat{U}_{r-1}$ $\displaystyle=$
$\displaystyle~{}\exp\Bigg{\\{}\,\frac{1}{2}\,\sum_{i,j}\Bigg{[}~{}\xi^{ss{(r-1)}}_{ij}\,\Big{(}\hat{\sigma}_{i}^{{}_{(r-1)}\dagger}\hat{\sigma}_{j}^{{}_{(r-1)}\dagger}-\hat{\sigma}_{i}^{{}_{(r-1)}}\hat{\sigma}_{j}^{{}_{(r-1)}}\Big{)}\,+$
$\displaystyle~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\xi^{ww(r-1)}_{ij}\Big{(}\hat{\omega}_{i}^{{}_{(r-1)}\dagger}\hat{\omega}_{j}^{{}_{(r-1)}\dagger}-\hat{\omega}_{i}^{{}_{(r-1)}}\hat{\omega}_{j}^{{}_{(r-1)}}\Big{)}\,+$
$\displaystyle~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\xi^{sw(r-1)}_{ij}\Big{(}\hat{\sigma}_{i}^{{}_{(r-1)}\dagger}\hat{\omega}_{j}^{{}_{(r-1)}\dagger}-\hat{\sigma}_{i}^{{}_{(r-1)}}\hat{\omega}_{j}^{{}_{(r-1)}}\Big{)}~{}\Bigg{]}\,\Bigg{\\}}.\qquad\qquad$

Above, the $\xi^{(r-1)}$ exponents are the WTs of the entangling exponent
$\zeta^{(r)}$,

$\displaystyle\xi^{ss(r-1)}_{ij}~{}$ $\displaystyle\equiv$
$\displaystyle~{}~{}\big{(}\mathbf{L}\,\zeta^{(r)}\mathbf{L}^{\mathsf{T}}\big{)}_{ij}$
(125a) $\displaystyle\xi^{ww(r-1)}_{ij}~{}$ $\displaystyle\equiv$
$\displaystyle~{}~{}\big{(}\mathbf{H}\,\zeta^{(r)}\mathbf{H}^{\mathsf{T}}\big{)}_{ij}$
(125b) $\displaystyle\xi^{sw(r-1)}_{ij}~{}$ $\displaystyle\equiv$
$\displaystyle~{}~{}2\,\big{(}\mathbf{L}\,\zeta^{(r)}\mathbf{H}^{\mathsf{T}}\big{)}_{ij}\,.$
(125c)

Note that they too have a restricted form due to discrete translational
invariance:

$\xi^{ss(r)}_{ij}\,=~{}\xi_{|i-j|}^{ss(r)},~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\xi^{ww(r)}_{ij}\,=~{}\xi_{|i-j|}^{ww(r)},~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\xi^{sw(r)}_{ij}\,=~{}\xi_{i-j}^{sw(r)}.$
(126)

The IWT relating $\zeta^{(r)}$ to $\xi^{(r-1)}$ is given by

$\zeta^{(r)}_{ij}~{}=~{}\big{(}\mathbf{L}^{\mathsf{T}}\xi^{ss(r-1)}\mathbf{L}\big{)}_{ij}\,+\,\big{(}\mathbf{H}^{\mathsf{T}}\xi^{ww(r-1)}\mathbf{H}\big{)}_{ij}\,+\,\frac{1}{2}\big{(}\mathbf{L}^{\mathsf{T}}\xi^{sw(r-1)}\mathbf{H}\,+\,\mathbf{H}^{\mathsf{T}}\xi^{sw(r-1)\mathsf{T}}\mathbf{L}\big{)}_{ij}\,.$
(127)

This expression is key for our algorithm to minimize
$\big{\langle}\hat{H}_{\sigma}^{r}\big{\rangle}$. Specifically, at a fixed
wMERA layer with resolution $r$, we treat $\xi_{i-j}^{ww(r)}$ and
$\xi_{i-j}^{sw(r)}$ as free parameters, and fix the condition

$\xi^{ss(r-1)}_{ij}~{}=~{}\zeta^{(r-1)}_{ij}\\!,$ (128)

that is, we use as input the numerical values for $\zeta^{(r-1)}$ obtained
from the optimization of the previous layer. While in the exact solution the
relation (128) is only approximately true, it is quite accurate for lattice
site separations $|i-j|\gtrsim 5$. For shorter lattice site separations
$|i-j|\leq 4$, the optimization of the free parameters $\xi^{ww(r)}$ and
$\xi^{sw(r)}$ can compensate for any inaccuracy introduced by (128).

With the parametrization of the entanglers now fully specified, we can proceed
to write the vacuum expectation value to the Hamiltonian in (123) in terms of
$\zeta^{(r)}_{ij}$. First, we obtain the entangler-transformed Hamiltonian in
(80) (with $\mathcal{V}=0$) using (120):

$\displaystyle\hat{U}_{r}^{\dagger}\hat{H}_{\sigma}^{r}\,\hat{U}_{r}~{}=~{}\,\frac{m_{\phi}}{2}\,\sum_{i,j}$
$\displaystyle\\!\\!\\!\Bigg{[}\Big{(}\text{C}_{2\zeta}^{(r)}+2^{(2r-1)}\mathds{K}^{ss}\,\text{E}_{2\zeta}^{(r)}\Big{)}_{\\!ij}\,\Big{(}{\hat{\sigma}_{i}}^{\,r\,\dagger}\hat{\sigma}_{j}^{\,r}+\hat{\sigma}_{j}^{\,r}\,\hat{\sigma}^{\,r\,\dagger}_{i}\Big{)}$
(129)
$\displaystyle~{}~{}+~{}\Big{(}\text{S}_{2\zeta}^{(r)}+2^{(2r-1)}\mathds{K}^{ss}\,\text{E}_{2\zeta}^{(r)}\Big{)}_{\\!ij}\,\Big{(}\hat{\sigma}^{\,r\,\dagger}_{i}\hat{\sigma}^{\,r\,\dagger}_{j}+\hat{\sigma}_{i}^{\,r}\,\hat{\sigma}_{j}^{\,r}\Big{)}\Bigg{]},~{}~{}~{}~{}~{}~{}~{}~{}~{}$

where, above, we have used the short-hand notation

$\displaystyle\big{(}\text{C}_{2\zeta}^{(r)}\big{)}_{ij}$
$\displaystyle\equiv$
$\displaystyle\big{(}\text{Cosh}\,2\zeta^{(r)}\big{)}_{ij}~{}~{}~{}~{}~{}~{}~{}~{}\big{(}\text{S}_{2\zeta}^{(r)}\big{)}_{ij}\equiv\big{(}\text{Sinh}\,2\zeta^{(r)}\big{)}_{ij}~{}~{}~{}~{}~{}~{}~{}~{}\big{(}\text{E}_{2\zeta}^{(r)}\big{)}_{ij}\equiv\big{(}\text{Exp}\,2\zeta^{(r)}\big{)}_{ij}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}$
$\displaystyle=$
$\displaystyle\big{(}\text{Cosh}\,2\zeta^{(r)}\big{)}_{i-j}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\;=\big{(}\text{Sinh}\,2\zeta^{(r)}\big{)}_{i-j}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}=\big{(}\text{Exp}\,2\zeta^{(r)}\big{)}_{i-j\,.}$

Using (129) and (83), we then have

$\displaystyle{\big{\langle}}\mathbf{0_{\sigma}}\big{|}\,\hat{H}_{\sigma}^{r}\,\big{|}{\mathbf{0_{\sigma}}}\big{\rangle}_{\\!r}~{}$
$\displaystyle=$
$\displaystyle~{}{\big{\langle}}\mathbf{P}_{\\!\sigma}\big{|}\,\hat{U}_{r}^{\dagger}\,\hat{H}_{\sigma}^{r}\,\hat{U}_{r}\,\big{|}{\mathbf{P}_{\\!\sigma}}\big{\rangle}_{\\!r}$
(131) $\displaystyle=$
$\displaystyle~{}\frac{m_{\phi}}{2}\;\text{Tr}\left[\text{C}_{2\zeta}^{(r)}+2^{(2r-1)}\mathds{K}^{ss}\,\text{E}_{2\zeta}^{(r)}\right].$

As it turns out, the minimization of this vacuum expectation value can be
performed analytically. Differentiating (131) with respect $\zeta^{(r)}$ and
setting it to zero, we obtain

$\text{S}_{2\zeta}^{(r)}+2^{(2r-1)}\mathds{K}^{ss}\,\text{E}_{2\zeta}^{(r)}~{}=~{}0,$
(132)

whose solution is given by

$\zeta^{(r)}_{ij}\,=\,-\frac{1}{4}\,\Big{(}\text{Log}\big{[}2^{2r}\,\mathds{K}^{ss}+\mathds{1}\big{]}\Big{)}_{\\!ij}\,.$
(133)

However, as previously stated, our actual goal in this demonstration of the
wMERA algorithm is to minimize (131) _numerically_ using the parametrization
in (127) and the input in (128); we will only use the exact solution for
$\zeta^{(r)}_{ij}$ in (133) to validate our numerical results. We can further
manipulate (131) using (5.1), (52), and (53) to obtain:

${\big{\langle}}\mathbf{0_{\sigma}}\big{|}\,\hat{H}_{\sigma}^{r}\,\big{|}{\mathbf{0_{\sigma}}}\big{\rangle}_{\\!r}~{}=~{}\frac{m_{\phi}}{2}\,\sum_{i=\infty}^{\infty}\Bigg{[}\big{(}\text{Cosh}\,2\zeta^{(r)}\big{)}_{\\!0}+\,2^{(2r-1)}\\!\sum_{j=-4}^{4}K_{|j|}\,\big{(}\text{Exp}\,2\zeta^{(r)}\big{)}_{\\!j}\Bigg{]}.$
(134)

The expression above has two remarkable aspects. First, even though it is a
trace over a square matrix of infinite dimension, the infinite summation over
one of the indices factors out (this is due to discrete translational
invariance). Second, the r.h.s. of (134) depends only on a small, finite
number of elements of the matrix exponential of $\zeta^{(r)}$ (this is due to
the Daubechies discretization, whose scaling and wavelet functions have
compact support in position space). In particular, because of the form of
(134), there is no need to perform a finite volume truncation of the system in
order to minimize $\big{\langle}\hat{H}_{\sigma}^{r}\big{\rangle}$ with
respect to the variational parameters. However, a truncation is still needed
in the number of variational parameters, of which there is in principle an
infinite number. One may question whether any such truncation would preserve
the properties of the ground state. To argue that this is the case, note that
the matrix $\zeta^{(r)}_{i-j}$ is the discrete counterpart of the entangler in
the continuum, $\zeta(x-y)$, whose exact solution is given in (114). In
particular, $\zeta(x-y)$ becomes exponentially suppressed for separations
larger than the correlation length of the system, i.e., for $|x-y|\gtrsim
m_{\phi}^{-1}$. Therefore, we expect that the entries $\zeta^{(r)}_{i-j}$ will
likewise be exponentially suppressed for site separations $|i-j|\gtrsim
2^{r}$.

Recalling that the free parameters of $\zeta^{(r+1)}$ are given by
$\xi^{ww(r)}$ and $\xi^{sw(r)}$ as defined in (127), in practice our
truncation of free parameters will amount to setting

$\displaystyle\xi^{ww(r)}_{i}$ $\displaystyle=\;$ $\displaystyle
0~{}~{}~{}\text{for}~{}~{}|i|>i^{(r)}_{\text{max,}}$ (135a)
$\displaystyle\xi^{sw(r)}_{j}$ $\displaystyle=\;$ $\displaystyle
0~{}~{}~{}\text{for}~{}~{}|j|>j^{(r)}_{\text{max,}}$ (135b)

where $i^{(r)}_{\text{max}}$ and $j^{(r)}_{\text{max}}$ are finite integers
which can be increased to improve numerical accuracy. Remarkably, we will see
shortly that $i^{(r)}_{\text{max}}$ and $j^{(r)}_{\text{max}}$ do not need to
grow with resolution $r$, and can be kept small (i.e., in the range of $\sim$
3–5) for all wMERA layers while still retaining excellent accuracy999The
parallel of this feature in the original MERA is the fact that the bond
dimension of the tensor ansatzes do not need to grow with the resolution of
the MERA layer Vidal:2007hda .. This is in contrast with the truncation of
entangling parameters $\zeta^{(r)}_{i}$, which must grow with $r$ to allow the
same desired degree of accuracy to be retained, as we just argued in the
previous paragraph.

Figure 7: The optimized values of the entangling parameters
$\zeta^{(r)}_{ij}$, $\xi^{ww(r)}_{ij}$, and $\xi^{sw(r)}_{ij}$, defined in
(118) and (125), for a few illustrative choices of resolution $r$. We recall
that the lattice spacing is related to the resolution index $r$ via
$a_{r}=1/(2^{r}m_{\phi})$.

In appendix A, we discuss a method to obtain the entries of
$\big{(}\text{Exp}\,2\zeta\big{)}_{\\!j}$ as a function of a truncated array
$\zeta_{i}$. This, combined with (127), can plugged into (134) for the
numerical variational optimization.

The exact solutions for $\zeta^{(r)}_{i-j}$, $\xi^{ww(r)}_{i-j}$, and
$\xi^{sw(r)}_{i-j}$ are shown in figure 7 for a few illustrative choices of
resolution $r$. The plots clearly show, as previously stated, that the range
of the entangling parameters $\xi^{ww(r)}$ and $\xi^{sw(r)}$ does not grow
exponentially with increasing resolution. Indeed, their optimal values not
only quickly saturate at resolutions close to $r=0$, but their range decays
exponentially with increasing lattice site separations, which justifies the
truncations in (135). This is in stark contrast to the scaling entangling
parameters $\zeta^{(r)}$, whose range grows exponentially with resolution.

The results of our numerical wMERA optimization show excellent agreement with
the exact solutions within the truncated range of parameters in our algorithm.
Specifically, in obtaining our results, we truncated the $\xi^{ww}_{i}$,
$\xi^{sw}_{j}$, and $\zeta_{k}$ arrays at $i_{\text{max}}=3$,
$j_{\text{max}}=4$, and $k_{\text{max}}=21$, respectively, in all the layers
of our calculation, which ranged from $r=-$3 to 10.

For a fixed resolution layer $r$, our numerical results for $\zeta^{(r)}_{ij}$
allow us to recover the continuum result for the entangling exponent
$\zeta(x-y)$ restricted to the range
$|x-y|\sim\mathcal{O}(1-10)\times(2^{r}m_{\phi})^{-1}$:

$\displaystyle\zeta^{(r)}(x-y)~{}~{}$ $\displaystyle=$
$\displaystyle~{}~{}\sum_{i,j}\,s^{r}_{i}(x)s^{r}_{j}(y)\;\zeta^{(r)}_{ij}$
$\displaystyle=$
$\displaystyle~{}~{}\sum_{i,j}\,m_{\phi}\;2^{r}\,s\big{(}2^{r}m_{\phi}\,x-i\big{)}\,s\big{(}2^{r}m_{\phi}\,y-j\big{)}\;\zeta^{(r)}_{i-j}\,.$

Piecing together the functions $\zeta^{(r)}(x-y)$ at different resolutions, we
can reconstruct the continuum entangling exponent $\zeta(x-y)$ over several
orders of magnitude in $|x-y|$, as shown in figure 8.

Our results for $\zeta^{(r)}_{ij}$ also allow us to reconstruct the scalar
field two-point function in the continuum. For spacelike separations, it
suffices to compute the equal-time two-point function; the answer for any
other spacelike interval can be obtained by a Lorentz transformation.
Recalling the DWT representation of the scalar field $\hat{\Phi}^{r}(x)$ at a
fixed resolution $r$ in (70a), we have:

$\displaystyle\big{\langle}\hat{\Phi}^{r}(x)\hat{\Phi}^{r}(y)\big{\rangle}~{}~{}$
$\displaystyle=$
$\displaystyle~{}~{}\sum_{i,j}\,s^{r}_{i}(x)s^{r}_{j}(y)~{}{\big{\langle}}\mathbf{0_{\sigma}}\big{|}\,\hat{\phi}_{i}^{r}\,\hat{\phi}_{j}^{r}\,\big{|}{\mathbf{0_{\sigma}}}\big{\rangle}_{\\!r}$
$\displaystyle=$
$\displaystyle~{}~{}\sum_{i,j}\,s^{r}_{i}(x)s^{r}_{j}(y)~{}{\big{\langle}}\mathbf{P_{\sigma}}\big{|}\big{(}\hat{U}_{r\,}^{\dagger}\hat{\phi}_{i}^{r}\,\hat{U}_{r}\big{)}\,\big{(}\hat{U}_{r\,}^{\dagger}\hat{\phi}_{j}^{r}\,\hat{U}_{r}\big{)}\big{|}{\mathbf{P_{\sigma}}}\big{\rangle}_{\\!r}$
$\displaystyle=$
$\displaystyle~{}~{}\sum_{i,j}\,2^{(r-1)}\,s\big{(}2^{r}m_{\phi}\,x-i\big{)}\,s\big{(}2^{r}m_{\phi}\,y-j\big{)}\;\big{(}\text{Exp}\,2\zeta^{(r)}\big{)}_{i-j\,,}$

where, in the last equality, we used (66), (72a), (83), and (120).

Similarly to the case of $\zeta^{(r)}(x-y)$, the two-point function in (5.1)
is restricted to the range
$|x-y|\sim\mathcal{O}(1-10)\times(2^{r}m_{\phi})^{-1}$. Piecing together
$\big{\langle}\hat{\Phi}^{r}(x)\hat{\Phi}^{r}(y)\big{\rangle}$ at different
resolutions, we recover the continuum result for
$\big{\langle}\hat{\Phi}(x)\hat{\Phi}(y)\big{\rangle}$ over several orders of
magnitude in $|x-y|$, as shown in figure 8.

Figure 8: _Top:_ the reconstructed continuum entangling exponent $\zeta(x-y)$
(see (104)), calculated using the numerical output of our wMERA optimization.
The solid black curve shows the exact continuum solution for $\zeta(x-y)$ (see
(114)). _Bottom:_ the reconstructed continuum scalar two-point function at
equal times, $\big{\langle}\hat{\Phi}(x)\hat{\Phi}(y)\big{\rangle}$,
calculated using the numerical output of our wMERA optimization. The exact
continuum solution for $\big{\langle}\hat{\Phi}(x)\hat{\Phi}(y)\big{\rangle}$
is given by the solid black curve (see (144)). In both plots, the functions
were calculated for different ranges in $m_{\phi}|x-y|$ using different wMERA
layers, indicated by the resolution index $r$. We recall that the lattice
spacing is related to the resolution index $r$ via $a_{r}=1/(2^{r}m_{\phi})$.

Our final task in solving this free scalar theory with wMERA is to obtain the
scalar two-point function at finite-time separations. This requires
implementing real-time evolution, which is also straightforward in the wMERA
framework. The first step is to obtain the entangler-transformed Hamiltonians
at different resolutions using the variational solution for the entangler
parameters. Plugging (133) into (129), we obtain:

$\hat{U}_{r}^{\dagger}\hat{H}_{\sigma}^{r}\,\hat{U}_{r}~{}=~{}m_{\phi}\,\sum_{i,j}\;{\hat{\sigma}_{i}}^{\,r\,\dagger\,}\Omega^{(r)}_{ij}\,\hat{\sigma}_{j}^{\,r}~{}+~{}\text{(zero-
point energy)}$ (138)

where

$\Omega^{(r)}_{ij}~{}\equiv~{}\Big{(}e^{-2\zeta^{(r)}}\Big{)}_{ij\,.}$ (139)

We can also define the entangler-transformed time-evolution operator101010On
the r.h.s. of (140), we omitted a constant phase stemming from the zero-point
energy.:

$\hat{\mathds{U}}_{r}(t)~{}\equiv~{}\hat{U}_{r}^{\dagger}\Big{(}e^{-i\,t\,\hat{H}_{\sigma}^{r}}\Big{)}\hat{U}_{r}~{}=~{}\text{Exp}\bigg{[}\\!\\!-i\,(m_{\phi}\,t)\sum_{i,j}{\hat{\sigma}_{i}}^{\,r\,\dagger\,}\Omega^{(r)}_{ij}\,\hat{\sigma}_{j}^{\,r}\bigg{]}.$
(140)

Its action on the scaling modes at resolution $r$ is easily derived:

$\displaystyle\,\hat{\mathds{U}}_{r}^{\dagger}(t)\,\hat{\sigma}_{i}^{\,r}\,\hat{\mathds{U}}_{r}(t)$
$\displaystyle~{}=~{}$
$\displaystyle\big{(}e^{-i\,(m_{\phi}t)\,\Omega^{(r)}}\big{)}_{ij}\,\hat{\sigma}_{j}^{\,r}$
(141a)
$\displaystyle\hat{\mathds{U}}_{r}^{\dagger}(t)\,\hat{\sigma}_{i}^{\,r\dagger}\;\hat{\mathds{U}}_{r}(t)$
$\displaystyle~{}=~{}$
$\displaystyle\big{(}e^{+i\,(m_{\phi}t)\,\Omega^{(r)}}\big{)}_{ij}\,\hat{\sigma}_{j}^{\,r\,\dagger}.$
(141b)

We now have all the ingredients to obtain the scalar field two-point function
at resolution $r$ for finite-time separations :

$\displaystyle\big{\langle}\hat{\Phi}^{r}(x,t)\hat{\Phi}^{r}(y,0)\big{\rangle}~{}$
$\displaystyle=$
$\displaystyle~{}{\big{\langle}}\mathbf{0_{\sigma}}\big{|}\big{(}e^{i\,t\,\hat{H}_{\sigma}^{r}}\big{)}\hat{\Phi}^{r}(x)\big{(}e^{-i\,t\,\hat{H}_{\sigma}^{r}}\big{)}\hat{\Phi}^{r}(y)\big{|}{\mathbf{0_{\sigma}}}\big{\rangle}_{\\!r}$
$\displaystyle=$
$\displaystyle~{}\sum_{i,j}\,s^{r}_{i}(x)s^{r}_{j}(y)~{}{\big{\langle}}\mathbf{P_{\sigma}}\big{|}\hat{\mathds{U}}_{r}^{\dagger}(t)\,\big{(}\hat{U}_{r\,}^{\dagger}\hat{\phi}_{i}^{r}\,\hat{U}_{r}\big{)}\,\hat{\mathds{U}}_{r}(t)\,\big{(}\hat{U}_{r\,}^{\dagger}\hat{\phi}_{j}^{r}\,\hat{U}_{r}\big{)}\big{|}{\mathbf{P_{\sigma}}}\big{\rangle}_{\\!r}$
$\displaystyle=$
$\displaystyle~{}\sum_{i,j}\,2^{(r-1)}\,s\big{(}2^{r}m_{\phi}\,x-i\big{)}\,s\big{(}2^{r}m_{\phi}\,y-j\big{)}\Big{(}e^{2\zeta^{(r)}-i\,(m_{\phi}t)\,\Omega^{(r)}}\Big{)}_{i-j\,.}$

The expression above, which generalizes (5.1), is restricted to the range
$\sqrt{|t^{2}-(x-y)^{2}|}\sim\mathcal{O}(1-10)\times(2^{r}m_{\phi})^{-1}$.
Piecing together $\big{\langle}\hat{\Phi}^{r}(x,\Delta
t)\hat{\Phi}^{r}(x,0)\big{\rangle}$ at different resolutions, we recover the
continuum result for $\big{\langle}\hat{\Phi}(x,\Delta
t)\hat{\Phi}(x,0)\big{\rangle}$ at timelike separations over several orders of
magnitude in $\Delta t$, as shown in figure 9.

Figure 9: The reconstructed continuum scalar two-point function for timelike
separations, $\big{\langle}\hat{\Phi}(x,t)\hat{\Phi}(x,0)\big{\rangle}$,
calculated using the numerical output of our wMERA optimization. The real and
imaginary parts are shown in the top and bottom plots, respectively. The
apparent derivative discontinuity at $m_{\phi\,}t=1$ in the bottom plot is an
artifact of the hybrid scale of the horizontal axis: note that the scale is
logarithmic for $m_{\phi\,}t<1$, and linear for $m_{\phi\,}t>1$. This choice
was made so that the comparison between the numerical results and the exact
continuum solution (solid black curves) is displayed with better clarity. In
both plots, the functions were calculated for different ranges in
$m_{\phi\,}t$ using different wMERA layers, indicated by the resolution index
$r$. We recall that the lattice spacing is related to the resolution index $r$
via $a_{r}=1/(2^{r}m_{\phi})$.

It is important to note that our wMERA results for
$\big{\langle}\hat{\Phi}^{r}(x,t)\hat{\Phi}^{r}(y,0)\big{\rangle}$ should
track the continuum solution over scales longer than the resolution $r$, but
otherwise show some degree of unphysical structure at scales of order
$|x-y|\sim(2^{r}m_{\phi})^{-1}$. Indeed, in figures 8 and 9, we chose to plot
our results for a discrete array of points $\\{x_{i},y_{i}\\}$ satisfying

$\begin{cases}\;2^{r}m_{\phi}\,x_{i}\,\equiv\,\text{floor}(2^{r}m_{\phi}\,x_{i})+\epsilon\\\
\;2^{r}m_{\phi}\,y_{i}\;\equiv\,\text{floor}(2^{r}m_{\phi}\,y_{i})+\epsilon\end{cases}\\!\\!\\!\\!,~{}~{}~{}~{}~{}~{}~{}~{}2^{r}m_{\phi}(x_{i}-y_{i})~{}\equiv~{}\delta_{xy}\in\mathbb{Z}.$
(143)

With these definitions, we can rewrite (5.1) as:

$\displaystyle\big{\langle}\hat{\Phi}^{r}(x_{i},t)\hat{\Phi}^{r}(y_{i},0)\big{\rangle}~{}$
$\displaystyle=$
$\displaystyle~{}\sum_{j,k}\,2^{(r-1)}\,s\big{(}\delta_{xy}-k-j+\epsilon\big{)}\,s\big{(}j+\epsilon\big{)}\Big{(}e^{2\zeta^{(r)}-i\,(m_{\phi}t)\,\Omega^{(r)}}\Big{)}_{\\!k\,.}$

The unphysical structure of this reconstructed two-point function at scales
$|x_{i}-y_{i}|\sim(2^{r}m_{\phi})^{-1}$ can be seen by varying the parameter
$\epsilon$ above over its range of $0\leq\epsilon<1$ (in other words, varying
$\epsilon$ will probe the lack of smoothness of the D6 scaling function).
Hence, for every plot in this paper111111Namely, those in figures 8, 9, 10,
11, and 12. showing continuum functions reconstructed from numerical wMERA
results, we select the value of $\epsilon$ that minimizes these unphysical
fluctuations.

We can compare our numerical wMERA results for the scalar two-point function
with the exact analytical expression for the Feynman propagator in position
space, given by:

$\big{\langle}\hat{\Phi}(x,t)\hat{\Phi}(0,0)\big{\rangle}~{}=~{}\Theta\big{(}s^{2}\big{)}\,\frac{(-i)}{4}\,H_{0}^{(2)}\big{(}m_{\phi}\sqrt{s^{2}}\big{)}~{}+~{}\Theta\big{(}\\!\\!-\\!s^{2}\big{)}\,\frac{1}{2\pi}\,K_{0}\big{(}m_{\phi}\sqrt{-s^{2}}\big{)},$
(144)

where $t>0$, $\Theta$ is the Heaviside step function, $H_{0}^{(2)}$ is a
Hankel function of the second kind, $K_{0}$ is a modified Bessel function of
the second kind, and $s$ is the spacetime interval, i.e., $s^{2}\equiv
t^{2}-x^{2}$.

As shown in figures 8 and 9, the agreement between our numerical results and
the exact continuum solution is quite satisfactory considering our aggressive
truncation choices. That said, one can see that the accuracy and precision of
our results slowly degrade with increasing $r$. This is not a shortcoming of
wMERA itself, but an expected artifact of the approximations made in our
numerical code. Specifically: (i) we kept the $\zeta^{(r)}_{k}$ array
truncation fixed at $k^{(r)}_{\text{max}}=21$ for all layers; (ii) we
restricted the order of the power expansion of the matrix exponential
$\text{Exp}\,2\zeta^{(r)}$ in the $\zeta^{(r)}_{|k|}$ coefficients to 3 (see
appendix A for details); and (iii) we truncated the
$\big{(}\text{Exp}\,2\zeta^{(r)}\big{)}_{l}$ array to
$l^{(r)}_{\text{max}}=36$ for all layers. These choices, which were not
optimized, were made so we could run our code on a personal laptop with a dual
Intel i5 core processor in a couple of hours (the code was written in Wolfram
Mathematica Mathematica ). However, with more computational resources, one can
relax the truncation choices to overcome the degradation in accuracy and
precision with increasing resolution.

###  5.2 Free theory of a Dirac fermion in (1+1)d

Our second and final wMERA example is a free theory of a Dirac fermion in
(1+1) spacetime dimensions. Many steps of this exercise (and their
justifications) have a strong parallel with the scalar field theory case just
considered. Therefore, in this section, we will skip over details that have
already been dissected in the previous section.

We start by laying out our notation. We choose conventions such that the Dirac
spinor and $\gamma$-matrices are given by:

$\hat{\Psi}=\begin{pmatrix}\\!\hat{\chi}\\\
\,\hat{\psi}^{\dagger}\end{pmatrix}~{},~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\gamma^{0}=\begin{pmatrix}\,1&~{}\,0\,\\\
\,0&-1\,\end{pmatrix}~{},~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\gamma^{1}=\begin{pmatrix}\\!~{}\,0&~{}1\,\\\
\\!-1&~{}0\,\end{pmatrix}.$ (145)

The system’s dynamics is determined by the Hamiltonian

$\displaystyle\hat{H}~{}$ $\displaystyle=$ $\displaystyle~{}\int
dx\,\Big{[}\\!-i\,\hat{\overline{\Psi}}(x)\gamma^{1}\partial_{x}\hat{\Psi}(x)\;+\;m_{\psi}\,\hat{\overline{\Psi}}(x)\hat{\Psi}(x)\Big{]}$
$\displaystyle=$ $\displaystyle~{}\int
dx\,\Big{[}\\!-i\,\hat{\chi}^{\dagger}(x)\,\partial_{x}\hat{\psi}^{\dagger}(x)-i\,\hat{\psi}^{\dagger}(x)\,\partial_{x}\hat{\chi}^{\dagger}(x)+m_{\psi}\,\big{(}\hat{\chi}^{\dagger}(x)\hat{\chi}(x)-\hat{\psi}(x)\hat{\psi}^{\dagger}(x)\big{)}\Big{]}$

and equal-time anticommutation relations

$\\{\hat{\chi}(x),\,\hat{\chi}^{\dagger}(y)\\}=\delta(x-y),~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\\{\hat{\psi}(x),\,\hat{\psi}^{\dagger}(y)\\}=\delta(x-y);$
(147)

all other equal-time anticommutators vanishing.

We define the initial product state upon which the wMERA will be built as the
unentangled state in position space that is annihilated by $\hat{\chi}$ and
$\hat{\psi}$:

$\hat{\chi}(x)\,\big{|}\mathbf{P}\big{\rangle}\,=\;\hat{\psi}(x)\,\big{|}\mathbf{P}\big{\rangle}~{}=~{}0~{}~{}~{}~{}\forall~{}~{}~{}x.$
(148)

Our goal is to find the entangler $\hat{U}$ that relates the product state
$\big{|}{\mathbf{P}\big{\rangle}}$ to the Hamiltonian ground state
$\big{|}{\mathbf{0}\big{\rangle}}$:

$\big{|}{\mathbf{0}\big{\rangle}}~{}=~{}\hat{U}\big{|}{\mathbf{P}\big{\rangle}}.$
(149)

The ansatz for $\hat{U}$ is a Gaussian operator (serving the analogous role
for fermions that the generator of Bogoliubov transformations does for scalars
Khanna_2007 ). Its parametric form in generic ($d$+1)-spacetime dimensions
(with $d$ odd) is given by:

$\hat{U}~{}=~{}\exp\left[-i\int
dx^{d}\,dy^{d}~{}\,\zeta(|\vec{r}_{\,}|)~{}\,\hat{\Psi}^{\dagger\\!}(\vec{x}_{\,})\,(\hat{r}\cdot\vec{\gamma})\,\hat{\Psi}(\vec{y}_{\,})\,\right]\\!,$
(150)

where $\vec{r}\equiv\vec{x}-\vec{y}$. For $d$ = 1 and in terms of the Dirac
spinor components in (145), (150) becomes

$\boxed{~{}\hat{U}~{}=~{}\exp\left[-i\int
dx\,dy~{}\zeta(x-y)~{}\Big{(}\hat{\chi}^{\dagger}(x)\hat{\psi}^{\dagger}(y)+\hat{\psi}(y)\hat{\chi}(x)\Big{)}\,\right]\,}.$
(151)

A key difference compared to the scalar case is that, for fermions, the
entangling exponent in (151) is anti-symmetric, that is,
$\zeta(y-x)=-\zeta(x-y)$.

To proceed with the variational optimization of the entangler, we go to
Fourier space. Defining

$\displaystyle\hat{\chi}_{k}$ $\displaystyle\;\equiv\,\int
dx~{}e^{-ikx}\,\hat{\chi}(x)~{},~{}~{}~{}~{}~{}~{}~{}$
$\displaystyle\hat{\psi}_{k}\;\equiv\,\int dx~{}e^{-ikx}\,\hat{\psi}(x)$
(152a) $\displaystyle\hat{\chi}^{\dagger}_{k}$ $\displaystyle\;\equiv\,\int
dx~{}e^{ikx}\,\hat{\chi}^{\dagger}(x)~{},~{}~{}~{}~{}~{}~{}~{}$
$\displaystyle\hat{\psi}^{\dagger}_{k}\;\equiv\,\int
dx~{}e^{ikx}\,\hat{\psi}^{\dagger}(x)$ (152b)

and

$\zeta_{k}~{}\equiv~{}i\\!\int dx~{}e^{-ikx}\,\zeta(x),$ (153)

we can write the Hamiltonian as:

$\hat{H}~{}=~{}\int\frac{dk}{2\pi}~{}~{}\Big{[}k\,\big{(}\hat{\chi}^{\dagger}_{k}\,\hat{\psi}^{\dagger}_{-k}+\hat{\psi}_{-k}\,\hat{\chi}_{k}\big{)}\;+\;m_{\psi}\,\big{(}\hat{\chi}^{\dagger}_{k}\hat{\chi}_{k}\,-\,\hat{\psi}_{k}\hat{\psi}^{\dagger}_{k}\big{)}\Big{]}.$
(154)

The entangler in Fourier space is given by:

$\hat{U}~{}=~{}\exp\left[-\,\int\frac{dk}{2\pi}~{}\zeta_{k}~{}\Big{(}\hat{\chi}^{\dagger}_{k}\,\hat{\psi}^{\dagger}_{-k}-\hat{\psi}_{-k}\,\hat{\chi}_{k}\Big{)}\,\right],$
(155)

where, above, we used the fact that $\zeta_{-k}=-\,\zeta_{\,k}$. The action of
the entangler on the spinor components is easily derived121212Recalling that,
in Fourier space, the anticommutation relations are given by
$\big{\\{}\hat{\chi}_{k},\,\hat{\chi}_{q}^{\dagger}\big{\\}}=2\pi\,\delta(k-q)$
and
$\big{\\{}\hat{\psi}_{k},\,\hat{\psi}_{q}^{\dagger}\big{\\}}=2\pi\,\delta(k-q)$.:

$\displaystyle\,~{}\hat{U}^{\dagger\,}\hat{\chi}_{k}\,\hat{U}$
$\displaystyle~{}=~{}$
$\displaystyle\text{cos}(\zeta_{k})\,\hat{\chi}_{k}~{}-~{}\text{sin}(\zeta_{k})\,\hat{\psi}_{-k}^{\dagger}$
(156a) $\displaystyle\hat{U}^{\dagger\,}\hat{\psi}_{-k}^{\dagger}\,\hat{U}$
$\displaystyle~{}=~{}$
$\displaystyle\text{sin}(\zeta_{k})\,\hat{\chi}_{k}~{}+~{}\text{cos}(\zeta_{k})\,\hat{\psi}_{-k}^{\dagger}.$
(156b)

With (154) and (156), we obtain the entangler-transformed Hamiltonian,

$\displaystyle\hat{U}^{\dagger}\hat{H}\,\hat{U}~{}=~{}\int\frac{dk}{2\pi}$
$\displaystyle\Bigg{[}\Big{(}k\,\text{cos}(2\zeta_{k})-m_{\psi}\,\text{sin}(2\zeta_{k})\Big{)}\,\Big{(}{\hat{\chi}_{k}}^{\dagger}\,\hat{\psi}^{\dagger}_{-k}+\hat{\psi}_{-k}\,\hat{\chi}_{k}\Big{)}$
$\displaystyle~{}~{}+~{}\Big{(}k\,\text{sin}(2\zeta_{k})+m_{\psi}\,\text{cos}(2\zeta_{k})\Big{)}\,\Big{(}\hat{\chi}^{\dagger}_{k}\,\hat{\chi}_{k}-\hat{\psi}_{k}\,\hat{\psi}_{k}^{\dagger}\Big{)}\Bigg{]}.$

Using (5.2) and (148), we then obtain the vacuum expectation value of the
Hamiltonian,

$\displaystyle\big{\langle}\mathbf{0}\big{|}\hat{H}\big{|}\mathbf{0}\big{\rangle}~{}$
$\displaystyle=$
$\displaystyle~{}{\big{\langle}}\mathbf{P}\big{|}\hat{U}^{\dagger}\hat{H}\,\hat{U}\big{|}\mathbf{P}\big{\rangle}$
(158) $\displaystyle=$ $\displaystyle\,-\int
dk\,\,\delta(k-k)\,\Big{(}k\,\text{sin}(2\zeta_{k})+m_{\psi}\,\text{cos}(2\zeta_{k})\Big{)}.$

Minimizing $\langle\hat{H}\rangle$ with respect to entangler exponent
$\zeta_{k}$, we get:

$\displaystyle\frac{\partial}{\partial\zeta_{k}}\big{\langle}\mathbf{0}\big{|}\hat{H}\big{|}\mathbf{0}\big{\rangle}=0~{}~{}~{}~{}$
$\displaystyle\Rightarrow$
$\displaystyle~{}~{}~{}~{}k\,\text{cos}(2\zeta_{k})-m_{\psi}\,\text{sin}(2\zeta_{k})~{}=~{}0$
(159) $\displaystyle\Rightarrow$
$\displaystyle~{}~{}~{}~{}~{}~{}~{}\boxed{~{}\zeta_{k}\,=\,\frac{1}{2}\,\text{arctan}\bigg{(}\frac{k}{m_{\psi}}\bigg{)}\,}\,.$

Using (153), we finally obtain the fermionic entangler exponent in position
space,

$\boxed{~{}\zeta(x-y)~{}=~{}\frac{\,e^{-m_{\psi}|x-y|}\,}{4(x-y)}\;}\,.$ (160)

Finally, we can plug in the exact entangler expression into (5.2) to write the
entangler-transformed Hamiltonian as

$\displaystyle\hat{U}^{\dagger}\hat{H}\,\hat{U}~{}~{}$ $\displaystyle=$
$\displaystyle~{}\int\frac{dk}{2\pi}\,\omega_{k}\,\Big{(}\hat{\chi}^{\dagger}_{k}\,\hat{\chi}_{k}-\hat{\psi}_{k}\,\hat{\psi}_{k}^{\dagger}\Big{)},$
(161a) $\displaystyle=$ $\displaystyle~{}\int
dx\,dy~{}\,\Omega(x-y)\,\Big{(}\hat{\chi}^{\dagger}(x)\,\hat{\chi}(y)-\hat{\psi}(x)\,\hat{\psi}^{\dagger}(y)\Big{)}$
(161b) $\displaystyle=$ $\displaystyle~{}\int
dx\,dy~{}\,\hat{\overline{\Psi}}(x)\Omega(x-y)\,\hat{\Psi}(y).$ (161c)

where $\Omega(x-y)$ has already been defined in (116).

We now move on to the discretized theory to build the wMERA. Using the
Daubechies D6 wavelet basis131313In this section, we will assume that the
reference lattice spacing at $r=0$ is given by $a_{0}=1/m_{\psi}$. Therefore,
our definition of scaling functions $s^{r}_{i}(x)$ in (66) should be modified
by the replacement $m_{\phi}\to m_{\psi}$. to smear the fermion fields, we
define:

$\hat{\Psi}_{i}^{r}~{}=~{}\begin{pmatrix}\\!\hat{\chi}_{i}^{r}\\\
\,\hat{\psi}_{i}^{r\,\dagger}\end{pmatrix}~{}\equiv~{}\int
dx~{}s_{i}^{r}(x)\begin{pmatrix}\hat{\chi}(x)\\\
\hat{\psi}^{\dagger}(x)\end{pmatrix}\\!,$ (162)

so that the fermionic field representation at resolution $r$ is given
by141414Fermionic _wavelet_ modes can be analogously defined. However,
similarly to the scalar field example, we will variationally optimize each
wMERA after the IWT is performed, and therefore will only be working with
fermionic _scaling_ modes.:

$\hat{\Psi}^{r}(x)~{}=~{}\begin{pmatrix}~{}\hat{\chi}^{r}(x)\\\
\,\hat{\psi}^{r\,\dagger}(x)\end{pmatrix}~{}\equiv~{}\sum_{i}~{}s_{i}^{r}(x)\begin{pmatrix}\\!\hat{\chi}^{r}_{i}\\\
\,\hat{\psi}^{r\,\dagger}_{i}\end{pmatrix}\\!.$ (163)

Plugging (163) into (5.2), we obtain the discretized Hamiltonian at resolution
$r$:

$\hat{H}^{r}~{}=~{}m_{\psi}\sum_{i}\bigg{[}\hat{\chi}^{r\,\dagger}_{i}\hat{\chi}^{r}_{i}\,-\,\hat{\psi}^{r\,}_{i}\hat{\psi}^{r\,\dagger}_{i}-\,i\,2^{r}\sum_{j}\Big{(}\hat{\chi}^{r\,\dagger}_{i}\,\mathds{D}^{ss}_{ij}\,\hat{\psi}^{r\,\dagger}_{j}\,+\,\hat{\psi}^{r}_{i}\,\mathds{D}^{ss}_{ij}\,\hat{\chi}^{r}_{j}\Big{)}\bigg{]},$
(164)

where the kinetic term coefficients $\mathds{D}^{ss}_{ij}$ have been defined
in (61) and (62).

The ansatz for the unitary entangler that transforms the initial product state
into the ground state of the discretized Hamiltonian in (164) can be written
as:

$\hat{U}_{r}~{}=~{}\hat{R}_{r}(\theta)\;\exp\\!\left[-\,i~{}\sum_{i,j}\,\zeta_{ij}^{(r)}\;\Big{(}\hat{\chi}_{i}^{\,r\,\dagger}\hat{\psi}_{j}^{\,r\,\dagger}+\,\hat{\psi}_{j}^{\,r}\hat{\chi}_{i}^{\,r}\Big{)}\,\right],$
(165)

where entangling exponent ${\zeta}^{(r)}$ is an _anti_ -Hermitian matrix, and
the operator $\hat{R}_{r}(\theta)$, with a single additional variational
parameter $\theta$, is given by:

$\hat{R}_{r}(\theta)~{}=~{}\exp\\!\left[-\,\theta~{}\sum_{i}\,\Big{(}\hat{\chi}_{i}^{\,r\,\dagger}\hat{\psi}_{i}^{\,r\,\dagger}-\,\hat{\psi}_{i}^{\,r}\hat{\chi}_{i}^{\,r}\Big{)}\,\right].$
(166)

This operator mixes the two components of the Dirac spinor _locally_. It turns
out that, since we are working with the Dirac representation of
$\gamma$-matrices, $\theta=0$. Had we chosen to work with the Weyl
representation, we would find that $\theta=\pi/4$. For simplicity, in our
wMERA exercise we will set $\theta=0$ from the start, although including
$\theta$ as a ansatz parameter in our wMERA optimization would have been
straightforward and would not have affected the final results.

Setting $\hat{R}_{r}(\theta)=\mathds{1}$ in (165), the action of the entangler
$\hat{U}_{r}$ on the fermionic scaling modes is given by:

$\displaystyle~{}\hat{U}_{r}^{\dagger}\,\hat{\chi}_{i}^{\,r}\,\hat{U}_{r}$
$\displaystyle~{}=~{}$
$\displaystyle\big{(}\text{Cosh}\zeta^{(r)}\big{)}_{ij}~{}\hat{\chi}_{j}^{\,r}~{}\,-~{}i\,\big{(}\text{Sinh}\zeta^{(r)}\big{)}_{ij}\,\hat{\psi}_{j}^{\,r\,\dagger}$
(167a)
$\displaystyle\hat{U}_{r}^{\dagger}\,\hat{\psi}_{i}^{\,r\,\dagger}\,\hat{U}_{r}$
$\displaystyle~{}=~{}$
$\displaystyle\big{(}\text{Cosh}\zeta^{(r)}\big{)}_{ij}\;\hat{\psi}_{j}^{\,r\,\dagger}+~{}i\,\big{(}\text{Sinh}\zeta^{(r)}\big{)}_{ij}\;\hat{\chi}_{j}^{\,r}.$
(167b)

As with the scalar example in subsection 5.1, given a wMERA layer $r$, our
free variational parameters $\xi^{ww(r-1)}$ and $\xi^{sw(r-1)}$ will be
embedded in $\zeta^{(r)}$ via an IWT:

$\zeta^{(r)}_{ij}~{}=~{}\big{(}\mathbf{L}^{\mathsf{T}}\xi^{ss(r-1)}\mathbf{L}\big{)}_{ij}\,+\,\big{(}\mathbf{H}^{\mathsf{T}}\xi^{ww(r-1)}\mathbf{H}\big{)}_{ij}\,+\,\frac{1}{2}\big{(}\mathbf{L}^{\mathsf{T}}\xi^{sw(r-1)}\mathbf{H}\,-\,\mathbf{H}^{\mathsf{T}}\xi^{sw(r-1)\mathsf{T}}\mathbf{L}\big{)}_{ij},$
(168)

while $\xi^{ss(r-1)}$ will be fixed by the result of the optimization of the
previous layer through

$\xi^{ss(r-1)}_{ij}~{}=~{}\zeta^{(r-1)}_{ij}.$ (169)

Similarly to the scalar example, the entangling matrices $\zeta^{(r)}$,
$\xi^{ww(r)}$ and $\xi^{sw(r)}$ have their forms restricted due to discrete
translational invariance. However, unlike the scalar example, $\zeta^{(r)}$
and $\xi^{ww(r)}$ are _anti_ -Hermitian:

$\zeta_{ij}^{(r)}=-\zeta_{ji}^{(r)}=\,\zeta_{i-j}^{(r)},~{}~{}~{}~{}~{}~{}~{}~{}~{}\xi^{ww(r)}_{ij}=-\xi_{ji}^{ww(r)}=\,\xi_{i-j}^{ww(r)},~{}~{}~{}~{}~{}~{}~{}~{}~{}\xi^{sw(r)}_{ij}=\;\xi_{i-j}^{sw(r)}.$
(170)

We now proceed to obtain the entangler-transformed Hamiltonian using (167):

$\displaystyle\hat{U}_{r}^{\dagger}\,\hat{H}^{r}\,\hat{U}_{r}~{}=~{}m_{\psi}\,\sum_{i,j}$
$\displaystyle\\!\\!\\!\Bigg{[}\Big{(}\text{C}_{2\zeta}^{(r)}+2^{r}\,\mathds{D}^{ss}\,\text{S}_{2\zeta}^{(r)}\Big{)}_{\\!ij}\,\Big{(}{\hat{\chi}_{i}}^{\,r\,\dagger}\hat{\chi}_{j}^{\,r}-\hat{\psi}_{i}^{\,r}\,\hat{\psi}^{\,r\,\dagger}_{j}\Big{)}$
(171)
$\displaystyle~{}-i~{}\Big{(}\text{S}_{2\zeta}^{(r)}+2^{r}\,\mathds{D}^{ss}\,\text{C}_{2\zeta}^{(r)}\Big{)}_{\\!ij}\,\Big{(}{\hat{\chi}_{i}}^{\,r\,\dagger}\hat{\psi}_{j}^{\,r\,\dagger}+\hat{\psi}_{i}^{\,r}\,\hat{\chi}^{\,r}_{j}\Big{)}\Bigg{]},~{}~{}~{}~{}~{}~{}~{}~{}~{}$

where we have used the short-hand notation defined in (5.1). The vacuum
expectation value of the Hamiltonian is then given by:

$\displaystyle{\big{\langle}}\mathbf{0}\big{|}\hat{H}^{r}\big{|}{\mathbf{0}}\big{\rangle}_{\\!r}~{}$
$\displaystyle=$
$\displaystyle~{}{\big{\langle}}\mathbf{P}\big{|}\,\hat{U}_{r}^{\dagger}\,\hat{H}^{r}\,\hat{U}_{r}\,\big{|}{\mathbf{P}}\big{\rangle}_{\\!r}$
(172) $\displaystyle=$
$\displaystyle-\,m_{\psi}\;\text{Tr}\left[\text{C}_{2\zeta}^{(r)}+2^{r}\,\mathds{D}^{ss}\,\text{S}_{2\zeta}^{(r)}\right].$

Differentiating (172) with respect $\zeta^{(r)}$ and setting it to zero, we
obtain

$\text{S}_{2\zeta}^{(r)}+2^{r}\,\mathds{D}^{ss}\,\text{C}_{2\zeta}^{(r)}~{}=~{}0,$
(173)

which admits a formal solution given by

$\zeta^{(r)}_{ij}\,=\,-\frac{1}{2}\,\big{(}\text{Tanh}^{-1}2^{r}\,\mathds{D}^{ss}\big{)}_{\\!ij}\,.$
(174)

This solution is well-defined when the argument of the ArcTanh matrix function
falls within its domain of convergence. This, however, is not of our concern
here, since we will minimize (172) _numerically_ using the parametrization in
(168).

Further manipulating (172) using (61) and (62), we obtain:

${\big{\langle}}\mathbf{0}\big{|}\hat{H}^{r}\big{|}{\mathbf{0}}\big{\rangle}_{\\!r}~{}=~{}m_{\psi}\,\sum_{i=\infty}^{\infty}\Bigg{[}\\!-\big{(}\text{Cosh}\,2\zeta^{(r)}\big{)}_{\\!0}+\,2^{(r+1)}\sum_{j=1}^{4}D_{j}\,\big{(}\text{Sinh}\,2\zeta^{(r)}\big{)}_{\\!j}\Bigg{]}.$
(175)

We use this expression in our optimization algorithm, supplemented by the
truncations $i_{\text{max}}=7$, $j_{\text{max}}=9$, and $k_{\text{max}}=21$
for the $\xi^{ww}_{i}$, $\xi^{sw}_{j}$, and $\zeta_{k}$ arrays, respectively,
in all the layers of our calculation, which ranged from $r=$ 0 to 7.

With the numerical results for the entangling exponents, we can now
reconstruct the fermionic two-point function. Before proceeding, we first
write explicitly the exact continuum solution for the Feynman propagator in
position space:

$\displaystyle\big{\langle}\hat{\Psi}(x,t)\hat{\overline{\Psi}}(0,0)\big{\rangle}\,$
$\displaystyle=$
$\displaystyle\,\begin{pmatrix}\big{\langle}\hat{\chi}\,\hat{\chi}^{\dagger}\big{\rangle}&~{}\,-\big{\langle}\hat{\chi}\,\hat{\psi}\big{\rangle}\\\
\big{\langle}\hat{\psi}^{\dagger}\hat{\chi}^{\dagger}\big{\rangle}&~{}\,-\big{\langle}\hat{\psi}^{\dagger}\hat{\psi}\big{\rangle}\end{pmatrix}$
$\displaystyle=$ $\displaystyle\,m_{\psi}\,\mathds{1}_{2\times
2\\!}\left[\Theta\big{(}s^{2}\big{)}\,\frac{(-i)}{4}\,H_{0}^{(2)}\big{(}m_{\phi}\sqrt{s^{2}}\big{)}+\Theta\big{(}\\!\\!-\\!s^{2}\big{)}\,\frac{1}{2\pi}\,K_{0}\big{(}m_{\phi}\sqrt{-s^{2}}\big{)}\right]$
$\displaystyle~{}+m_{\psi}^{2}\,x_{\mu}\gamma^{\mu\\!}\left[\Theta\big{(}s^{2}\big{)}\,\frac{(-1)}{4}\,\frac{H_{1}^{(2)}\big{(}m_{\phi}\sqrt{s^{2}}\big{)}}{m_{\psi}\sqrt{s^{2}}}+\Theta\big{(}\\!\\!-\\!s^{2}\big{)}\,\frac{i}{2\pi}\,\frac{K_{1}\big{(}m_{\phi}\sqrt{-s^{2}}\big{)}}{m_{\phi}\sqrt{-s^{2}}}\right]\\!,$

where $t\geq 0$ and $s^{2}\equiv t^{2}-x^{2}$. From (5.2) we see that we need
to compute three expectation values (namely,
$\big{\langle}\hat{\chi}\,\hat{\chi}^{\dagger}\big{\rangle}$,
$\big{\langle}\hat{\psi}^{\dagger}\hat{\psi}\big{\rangle}$, and
$\big{\langle}\hat{\chi}\,\hat{\psi}\big{\rangle}$) in order to compare our
wMERA results with the continuum solution.

Plugging (173) into (171), we obtain the entangler-transformed Hamiltonian in
terms of the optimized entangler parameters:

$\hat{U}_{r}^{\dagger}\,\hat{H}^{r}\,\hat{U}_{r}~{}=~{}m_{\psi}\,\sum_{i,j}\,\bigg{(}\,{\hat{\chi}_{i}}^{\,r\,\dagger}\,\Omega^{(r)}_{ij}\,\hat{\chi}_{j}^{\,r}-\hat{\psi}_{i}^{\,r}\,\Omega^{(r)}_{ij}\,\hat{\psi}^{\,r\,\dagger}_{j}\bigg{)},$
(177)

with

$\Omega^{(r)}_{ij}~{}\equiv~{}\Big{(}\text{C}_{2\zeta}^{(r)}+2^{r}\,\mathds{D}^{ss}\,\text{S}_{2\zeta}^{(r)}\Big{)}_{\\!ij\,.}$
(178)

The entangler-transformed time-evolution operator is then given by:

$\displaystyle\hat{\mathds{U}}_{r}(t)~{}$ $\displaystyle\equiv$
$\displaystyle~{}\hat{U}_{r}^{\dagger}\Big{(}e^{-i\,t\,\hat{H}^{r}}\Big{)}\hat{U}_{r}$
(179) $\displaystyle=$
$\displaystyle~{}\text{Exp}\bigg{[}\\!\\!-i\,(m_{\psi}\,t)\sum_{i,j}\bigg{(}\,{\hat{\chi}_{i}}^{\,r\,\dagger}\,\Omega^{(r)}_{ij}\,\hat{\chi}_{j}^{\,r}-\hat{\psi}_{i}^{\,r}\,\Omega^{(r)}_{ij}\,\hat{\psi}^{\,r\,\dagger}_{j}\bigg{)}\bigg{]}.$

It transforms the fermionic scaling modes as :

$\displaystyle\hat{\mathds{U}}_{r}^{\dagger}(t)\;\hat{\chi}_{i}^{\,r}\;\hat{\mathds{U}}_{r}(t)$
$\displaystyle~{}=~{}$
$\displaystyle\big{(}e^{-i\,(m_{\phi}t)\,\Omega^{(r)}}\big{)}_{ij}\,\hat{\chi}_{j}^{\,r}$
(180a)
$\displaystyle\hat{\mathds{U}}_{r}^{\dagger}(t)\;\hat{\psi}_{i}^{\,r}\;\hat{\mathds{U}}_{r}(t)$
$\displaystyle~{}=~{}$
$\displaystyle\big{(}e^{-i\,(m_{\phi}t)\,\Omega^{(r)}}\big{)}_{ij}\,\hat{\psi}_{j}^{\,r}.$
(180b)

Figure 10: The reconstructed continuum fermionic two-point function at equal
times calculated using the numerical output of our wMERA optimization. The top
and bottom plots show, respectively,
$\text{Tr}\big{[}\big{\langle}\hat{\Psi}(x)\hat{\overline{\Psi}}(y)\big{\rangle}\big{]}=\big{\langle}\hat{\chi}(x)\hat{\chi}^{\dagger}(y)\big{\rangle}-\big{\langle}\hat{\psi}^{\dagger}(x)\hat{\psi}(y)\big{\rangle}$
and
$\text{Tr}\big{[}i\,\gamma^{1}\big{\langle}\hat{\Psi}(x)\hat{\overline{\Psi}}(y)\big{\rangle}\big{]}=i\,\big{\langle}\hat{\chi}(x)\hat{\psi}(y)\big{\rangle}+i\,\big{\langle}\hat{\psi}^{\dagger}(x)\hat{\chi}^{\dagger}(y)\big{\rangle}$.
The exact continuum solution is given by the solid black curves (see (5.2)).
In both plots, the functions were calculated for different ranges in
$m_{\psi}(x-y)$ using different wMERA layers, indicated by the resolution
index $r$. We recall that the lattice spacing is related to the resolution
index $r$ via $a_{r}=1/(2^{r}m_{\psi})$.

Using (167) and (180), we can obtain the fermionic two-point correlation
functions at any resolution $r$. For instance,

$\displaystyle\big{\langle}\hat{\chi}^{r}(x,t)\hat{\chi}^{r\,\dagger\\!}(y,0)\big{\rangle}$
$\displaystyle=$
$\displaystyle{\big{\langle}}\mathbf{0}\big{|}\big{(}e^{i\,t\,\hat{H}_{\sigma}^{r}}\big{)}\hat{\chi}^{r}(x)\big{(}e^{-i\,t\,\hat{H}_{\sigma}^{r}}\big{)}\hat{\chi}^{r\,\dagger\\!}(y)\big{|}{\mathbf{0}}\big{\rangle}_{\\!r}$
$\displaystyle=$
$\displaystyle\sum_{i,j}\,s^{r}_{i}(x)s^{r}_{j}(y)~{}{\big{\langle}}\mathbf{P}\big{|}\hat{\mathds{U}}_{r}^{\dagger}(t)\,\big{(}\hat{U}_{r\,}^{\dagger}\,\hat{\chi}_{i}^{r}\,\hat{U}_{r}\big{)}\,\hat{\mathds{U}}_{r}(t)\,\big{(}\hat{U}_{r\,}^{\dagger}\,\hat{\chi}_{j}^{r\,\dagger}\hat{U}_{r}\big{)}\big{|}{\mathbf{P}}\big{\rangle}_{\\!r}$
$\displaystyle=$ $\displaystyle
m_{\psi}\sum_{i,j}2^{r\,}s\big{(}2^{r}m_{\psi}\,x-i\big{)}\,s\big{(}2^{r}m_{\psi}\,y-j\big{)}\Big{(}\text{Cosh}^{2}\zeta^{(r)}e^{-i\,(m_{\phi}t)\,\Omega^{(r)}}\Big{)}_{\\!i-j\,.}$

Analogous steps can be followed to obtain

$\displaystyle\big{\langle}\hat{\psi}^{r\,\dagger\\!}(x,t)\hat{\psi}^{r}(y,0)\big{\rangle}$
$\displaystyle=$ $\displaystyle
m_{\psi}\sum_{i,j}2^{r\,}s\big{(}2^{r}m_{\psi}\,x-i\big{)}\,s\big{(}2^{r}m_{\psi}\,y-j\big{)}\Big{(}\\!\\!-\\!\text{Sinh}^{2}\zeta^{(r)}e^{-i\,(m_{\phi}t)\,\Omega^{(r)}}\Big{)}_{\\!i-j}~{}~{}~{}~{}~{}~{}~{}~{}~{}$
(182)

and

$\displaystyle\big{\langle}\hat{\psi}^{r\,\dagger\\!}(x,t)\hat{\chi}^{r\,\dagger\\!}(y,0)\big{\rangle}$
$\displaystyle=$
$\displaystyle-\,\big{\langle}\hat{\chi}^{r}(x,t)\hat{\psi}^{r}(y,0)\big{\rangle}$
$\displaystyle=$
$\displaystyle~{}\frac{i}{2}\,m_{\psi}\sum_{i,j}2^{r\,}s\big{(}2^{r}m_{\psi}\,x-i\big{)}\,s\big{(}2^{r}m_{\psi}\,y-j\big{)}\Big{(}\text{Sinh}2\zeta^{(r)}e^{-i\,(m_{\phi}t)\,\Omega^{(r)}}\Big{)}_{\\!i-j\,.}$

As usual, the expressions above are restricted to their regimes of validity,
$\sqrt{|t^{2}-(x-y)^{2}|}\sim\mathcal{O}(1-10)\times(2^{r}m_{\phi})^{-1}$.

Figure 11: The reconstructed continuum fermionic two-point function
$\big{\langle}\hat{\chi}(x,t)\hat{\chi}^{\dagger}(x,0)\big{\rangle}$ for
timelike separations calculated using the numerical output of our wMERA
optimization. The real and imaginary parts are shown in the top and bottom
plots, respectively. The oscillatory behavior of the wMERA results around the
exact continuum solution, given by solid black curves, stems from presence of
fermion doublers in the spectrum of the discretized theory—this effect becomes
significant at short time scales $m_{\psi\,}t<1$. In both plots, the functions
were calculated for different ranges in $m_{\psi\,}t$ using different wMERA
layers, indicated by the resolution index $r$. We recall that the lattice
spacing is related to the resolution index $r$ via $a_{r}=1/(2^{r}m_{\psi})$.

Piecing together the reconstructed two-point functions (5.2)-(5.2) at
different resolutions, we can compare them with the exact continuum solution
for $\big{\langle}\hat{\Psi}(\Delta x,\Delta
t)\hat{\overline{\Psi}}(0,0)\big{\rangle}$ in (5.2) at both spacelike and
timelike separations over several orders of magnitude in $\Delta x$ and
$\Delta t$, respectively. Our results are shown in figures 10, 11, and 12. We
find very good agreement between the wMERA and exact solutions for spacelike
separations (see figure 10). For timelike separations, there is good agreement
at a resolutions comparable to or lower than the system’s correlation length,
$m_{\psi}^{-1}$. However, at shorter time scales, there are two competing
effects worth mentioning: the first is a loss of accuracy and precision
stemming from the truncation choices in our numerical algorithm. We saw this
effect in the scalar case as well, and it can be easily mitigated with
modified truncations choices. The second effect is an oscillatory behavior of
the wMERA two-point function around the exact solution (see figures 11 and
12). This effect is actually _physical_ , and it reflects the fact that the
continuum and discretized theories have a different spectrum of states—an
effect commonly know as the “fermion doubler problem” in discretized field
theories with fermions. Indeed, our Daubechies D6 discretization is subject to
the Nielsen-Ninomiya theorem Nielsen:1981hk ; Friedan:1982nk . This means
that, in the massless limit, our D6-discretized theory of Dirac fermions has
indeed twice the number of zero-modes compared with the continuum theory
Fries:2018pcd . This explains the presence of oscillatory features in our
wMERA results: since the Hamiltonian mass term becomes less important as the
system flows to the UV, the system’s correlation function should approach the
massless solution at very high resolutions; therefore, it should not be
surprising that the effects of the doubler states become more and more
pronounced at shorter and shorter time scales.

Figure 12: Similar to figure 11, but for
$-\big{\langle}\hat{\psi}^{\dagger}(x,t)\hat{\psi}(x,0)\big{\rangle}$.

It is notable, however, that modulo these UV oscillations due to doublers, the
two-point function for D6-discretized theory _actually tracks_ the shape of
the continuum solution, a feature that is not shared by the naive lattice
field theory discretization based on Haar wavelets. Indeed, for the sake of
comparison, we also obtained the wMERA solution for the Haar-discretized
fermionic theory, and found that (i) the resulting two-point function had a
completely different shape than the continuum solution—this also follows from
solving the Haar-discretized theory using other methods, such as discrete
Fourier transforms; and (ii) the real-time Hamiltonian evolution of the Haar-
discretized theory suffered from much higher noise than the real-time
evolution of the D6-discretized theory. This noise was present at _all scales_
and was not mitigated by increasing the truncations orders in the numerical
algorithm. This seems to indicate that the D6 discretization might be a better
choice than the traditional, Haar discretization of fermionic theories. It is
not clear to this author, however, whether this preliminary conclusion will
hold once these theories are properly regularized to remove the doubler
states.

Given the ubiquity of the fermion doubler problem in lattice field theory,
many methods have been developed, in various contexts, to rid the spectrum of
fermionic doubler states (see, e.g., Rothe:1992nt ). At least some of these
methods should be adaptable to the Daubechies-based discretization.
Exploration of these possibilities, however, is deferred to future work.

##  6 Conclusions and Outlook

Traditionally, the concept of entanglement has been completely obscured in the
analytical and numerical treatments of QFTs. More recently, entanglement has
emerged as an important element in the context of the holographic principle,
black holes, and gravity Hayden:2007cs ; Maldacena:2013xja ; Swingle:2014uza ;
Pastawski:2015qua ; Brown:2015lvg , and as a useful measure to quantify
complexity and decoherence in quantum field theories Berges:2017hne ;
Carney:2017jut ; Agon:2014uxa ; Boyanovsky:2018fxl ; Balasubramanian:2011wt ;
Calabrese:2004eu . In this article, we hope to have conveyed another important
role that entanglement can play, namely, as a _computational tool_ 151515This
insight has emerged in condensed matter theory and has led to the emerging
field of tensor networks. to develop new numerical methods for quantum field
theories, which could potentially lead to novel insights.

In this study we introduced wMERA, an adaptation of entanglement
renormalization for quantum field theories that preserves several key features
of the original MERA, namely, (i) discretized degrees of freedom localized in
position space; (ii) a discrete, layered architecture of the ansatz
paralleling an inverse RG flow; and (iii) the quasi-locality of entanglers at
each layer. We considered two simple free field theory examples for which we
demonstrated the concrete implementation of the wMERA algorithm in detail. In
particular, we showed that wMERA allowed us to fully extract the two-point
correlators for spacetime separations spanning several orders of magnitude.
Notably, we were able to efficiently perform real-time evolution over many
decades in $\Delta t$ by working in the appropriate resolution for the ansatz
and for the system’s Hamiltonian. In other words, the use of the RG flow
played a key role in the efficiency of these computations.

In these specific wMERA examples, however, we did depart from the original
MERA implementation in one respect: instead of having the entanglers acting on
the Hilbert space (_à la_ Schrödinger picture), we chose to have the
entanglers transform the field operators (_à la_ Heisenberg picture). This,