1 Introduction↩︎

Tensor network states are a central tool in the quest to solve the quantum many-body problem classically. Tailored to the low energy states of local Hamiltonians, they already allow the simulation of large many-body systems numerically, with matrix product states (MPS) [1]–[6] in \(1\) space dimension, and projected entangled-pair states (PEPS) [5], [7]–[11] in \(2\) dimensions and more.

As their name suggests, tensor network methods are tied to a lattice. To study continuum models, like relativistic and non-relativistic quantum field theories (QFT), one has two strategies. First, one can discretize the model, solve it with standard tensor network techniques numerically, and then extrapolate the results to the continuum limit. In many cases, this is impressively accurate [12]–[15] as one can use the broad range of optimization methods developed for the lattice. However, one loses the guarantees of a variational method: it is difficult to know for sure if results on a finer lattice are actually better than those on a coarser one. Alternatively, one can take the continuum limit of the ansatz itself, and apply it to the continuum model directly. This yields variational results that need not be extrapolated. In \(1\) space dimension, continuous matrix product states (CMPS) [16]–[18], and their relativistic refinement [19], [20] work well, particularly in the thermodynamic limit. But they require new ad hoc optimization methods, and not all of the tensor network toolbox has been translated to the continuum yet. In \(2\) space dimensions and more, the situation is even more dire, and the continuous equivalent of PEPS [21], [22] is not yet practically usable for numerics.

In this article, we aim to develop a third way: use a discrete ansatz to solve a continuum theory variationally. To decompose states of quantum fields in position space, we consider an orthonormal family of sufficiently regular Daubechies wavelets [23], [24], which are a powerful tool in signal processing and numerical analysis. These functions are arbitrarily refinable into a discrete basis of \(L^2(\mathbb{R})\), differentiable, and compactly supported. Taking a subset of them, the Daubechies scaling functions at a fixed resolution \(r\), one can build a Fock space, which is a subspace \(\mathcal{H}^r\) of the continuum Hilbert space \(\mathcal{H}\). This Fock space has a natural tensor product structure, which respects spatial locality. It is thus natural to write candidate low energy states within \(\mathcal{H}^r\) as MPS, which is what we propose doing.

These wavelet-MPS (wMPS) have appealing properties: they are compatible with standard MPS algorithms (and thus applicable to interacting models), remain fully variational, and can be iteratively fine grained. One can start with a variational approximation of the ground state at a coarse scale, obtained very cheaply, and then refine it into a better approximation at a scale twice as fine. In practice, we easily chain more than a dozen such successive refinements without numerical instabilities. Conceptually, the construction is not in any way tied to \(1\) space dimension, and is transparently extended to \(d\geq 2\) space dimensions, albeit with likely tremendous numerical costs.

There has already been considerable discussion of wavelets in QFT and tensor networks, which it is important to contrast with what we are proposing.

First, there is a natural connection between wavelets and tensor network states, via the multi-scale entanglement renormalization ansatz (MERA) [25], [26]. More precisely, it has been shown that wavelet transforms can be written as a quantum circuit corresponding to a free fermionic or bosonic MERA [27]–[31]. Since the MERA itself is connected to the AdS/CFT correspondence, this has provided insights into holography [32]–[34]. However, despite important analytical insights, this MERA-wavelet connection has been used numerically only for free field theories.

Second, wavelets have long been considered as a candidate to discretize quantum fields in a controllable way, both for analytical questions [35] and for simulation on quantum computers [36]. Perhaps most closely related to our work is Alves [37], which expands quantum fields in position space using Daubechies scaling functions and then uses a wavelet MERA to approximate the ground state. As in previous work connecting wavelets and MERA, practical applications remain restricted to free field theories. Beyond free theories, Daubechies wavelets have recently been used to discretize a Hamiltonian formulation of 1+1 dimensional \(\phi^4\) theory in momentum space [38], [39], and solve it with Hamiltonian truncation (i.e. exact diagonalization on a subspace of the total Hilbert space). Because the Hamiltonian is non-local in momentum, this latter construction is not natural to use with tensor networks, and its cost remains exponential in system size.

Finally, there are recent discretizations of QFT that combine with tensor network states well, but that do not use wavelets. In [40], the authors cut space into small segments (which is akin to using the lowest order Daubechies wavelets, or Haar wavelets) on which they approximately diagonalize the local Hamiltonian. This provides a natural discrete tensor factorization of the Hilbert space, allowing one to find the ground state with the density matrix renormalization group (DMRG). The method is precise, but not strictly variational: connecting the segments continuously adds an infinite penalty term to the Hamiltonian, which needs to be approximated. Most recently, as we were putting the final touches to the present article, Shankar, Van Acoleyen, and Haegeman [41] have proposed a finite element discretization that is strictly variational, and preserves the locality structure that tensor network states need. Their construction is quite different from and complementary to ours, with some benefits – their family of functions has smaller support – and drawbacks – the finite elements are not orthonormal, which requires a modification of standard tensor network algorithms.

In Section 2, we present the mathematics of Daubechies scaling functions (the first half of Daubechies wavelet theory), before detailing our variational wavelet MPS approach. In Section 3, we use this fixed-resolution approach to study the ground state properties of the Lieb-Liniger model and assess convergence. In Section 4, we introduce Daubechies wavelet functions (the second half of Daubechies wavelet theory), which allow us to leverage the power of multi-resolution analysis, and iteratively refine our states until they are numerically free of finite scale artifacts.

2 Discrete basis for quantum fields↩︎

In this section, we present the continuous Fock space of non-relativistic QFT, using the Lieb-Liniger model as an example. We then give three criteria for a good discrete orthonormal family of functions, adapted to variational MPS, which ultimately motivate our choice of Daubechies scaling functions. We then review some of their key properties.

2.1 The (non-relativistic) QFT Hilbert space↩︎

To make our variational method concrete, we first give the prototypical example to which it applies: the Lieb-Liniger model. Our method may be used for other models with the same non-relativistic kinetic term and sufficiently local interactions, including models with more particle species, without difficulty.

The Lieb-Liniger model is a theory of non-relativistic bosons on an interval \(I\), interacting with a contact potential. Its Hamiltonian takes the following second-quantized form: \[\label{eq:LL} H=\int_I \partial_x \hat{\psi}^\dagger \partial_x \hat{\psi} + c\, \hat{\psi}^\dagger \hat{\psi}^\dagger \hat{\psi} \hat{\psi} - \mu \, \hat{\psi}^\dagger \hat{\psi},\tag{1}\] where \(c\) is the interaction strength, \(\mu\) the chemical potential, and \(\hat{\psi}(x)\) the bosonic annihilation operator, with commutation relations \([\hat{\psi}(x),\hat{\psi}(y)^\dagger]=\delta(x-y)\) and \([\hat{\psi}(x),\hat{\psi}(y)]=[\hat{\psi}(x)^\dagger,\hat{\psi}(y)^\dagger]=0\). This model is integrable, which makes it a convenient benchmark: it is possible to compute energy density errors, and to cleanly separate finite scale from finite bond dimension effects by comparing with CMPS values.

The Hilbert space of this model is the Fock space \(\mathcal{F}(L^2(I))\), obtained by creating excitations on top of the free Fock vacuum \(\ket{0}\) annihilated by all \(\hat{\psi}(x)\). A generic state \(\ket{\Psi}\) is a sum of states with fixed particle numbers \[\begin{align} \ket{\Psi}=& \displaystyle\sum\limits_{m\geq 0} \int_{I^m} \mathrm{d}x_1...\mathrm{d}x_m \; \varphi_m(x_1,...,x_m) \\ &\quad \quad\quad \quad \times\, \hat{\psi}^\dagger (x_1) ... \hat{\psi}^\dagger (x_m) \ket{0},\label{eq:secondq} \end{align}\tag{2}\] where \(\varphi_m\) is the \(m\)-particle wavefunction. Our goal is to write a discretely parameterized MPS that belongs to this Fock space, and that approximates the true ground state \(\ket{\text{GS}}\) of \(H\) well, especially in the thermodynamic limit where \(I\rightarrow \mathbb{R}\).

2.2 Requirements for a discrete orthonormal family↩︎

To write an MPS in the QFT Hilbert space, we first need a discrete orthonormal family of states, with a simple tensor product structure, spanning a subspace \(\mathcal{H}^r\) of \(\mathcal{F}(L^2(\mathbb{R}))\) [where \(r\) is a fixed scale, associated to a typical width \(\Delta x = 1/2^r\)]. A natural approach is to take first a discrete orthonormal family of functions \(f_n^r \in L^2(\mathbb{R})\), and then consider the tensor product of the Fock spaces associated to each mode function \(f_n^r\): \(\mathcal{H}^r = \bigotimes_n \mathcal{F}(f_n^r)\). This is the approach we follow in this article. The question now becomes: what makes a good orthonormal family of \(L^2(\mathbb{R})\) for our purposes?

For our method to be ultimately efficient, we require three properties:

1. Locality – The Hamiltonian \(H\) is local, and we would like it to remain local once written in our discrete tensor product subspace. This will ensure that \(H\) can be written exactly in MPO form, that its ground state has area law entanglement (up to log corrections) for this tensor product structure, and thus that MPS approximate its ground state well. We thus demand that the discrete orthonormal family of functions we consider is compactly supported.

2. Regularity – The kinetic term of a non-relativistic QFT Hamiltonian, as in 1 , involves the first derivative of the field operator. Thus, to obtain finite kinetic energies, we demand that the functions we consider are differentiable. This is crucial to keep our method truly variational.

3. Inclusion – Refining in scale, i.e. using functions of smaller width, should allow us to systematically get better variational approximations. To get this property, we demand that the vector space spanned by functions at resolution \(r\) is included in the vector space spanned by functions at resolution \(r+1\): \(\text{span}\{f_n^r,\; n\} \subset \text{span}\{f_n^{r+1},\; n\}\). Indeed, if this is the case, the inclusion can be promoted to Fock spaces, and implies that the Hilbert space considered grows with \(r\): \(\mathcal{H}^{r}\subset\mathcal{H}^{r+1}\).

Daubechies scaling functions with parameter \(N\geq 6\) satisfy all the properties listed above. However, as a start, it is helpful to see why these conditions are non-trivial, and why other natural choices of functions, which have been considered previously, violate at least one of our requirements.

The simplest choice is naive discretization, where discrete ladder operators at position \(n\) and resolution \(r\) are defined in terms of the bosonic field operators as \[\hat{a}_n^r = 2^{r/2} \int_{2^{-r}n}^{2^{-r}(n+1)} \mathrm{d}x\, \hat{\psi}(x) = \int_{\mathbb{R}} \mathrm{d}x\, s_n^r(x) \, \hat{\psi}(x),\] where we have introduced the so-called Haar scaling functions (which are simply windows) \[s_n^r(x) = \begin{cases} 2^{r/2}, & \text{if } 2^{-r}n \leq x < 2^{-r}(n+1), \\ 0, & \text{otherwise}. \end{cases}\] These functions indeed form an orthonormal family at fixed \(r\), have finite support, and the Haar scaling functions at scale \(r\) are simple sums of Haar scaling functions at scale \(r+1\). However, they are not differentiable. Thus, all the states in this Hilbert space but the trivial vacuum have infinite energy. The standard fix is to regulate the kinetic term by taking finite differences instead of derivatives, but this means changing \(H\), and thus losing the strict variational nature of the approach.

At least if we consider a finite interval \(I=[0,L]\), a natural choice would be a Fourier decomposition. As in [42], one could expand the bosonic operator \(\hat{\psi}(x)\) in quantized plane waves, \[\hat{\psi}(x) = \sum_{k \in \mathbb{Z}}\frac{1}{\sqrt{L}} e^{2\pi i kx/L} \hat{a}_k .\] This generates wildly non-local interaction terms. For instance, \[\int_{[0,L]} \hat{\psi}^\dagger \hat{\psi}^\dagger \hat{\psi} \hat{\psi} = \frac{1}{L} \sum_{k_1,k_2,k_3} \hat{a}_{k_1}^\dagger \hat{a}_{k_2}^\dagger \hat{a}_{k_3} \hat{a}_{k_1+k_2-k_3}.\] This non-locality is the sign that entanglement in momentum generically follows a volume law, and thus the corresponding tensor product structure is not adapted to MPS when \(L\rightarrow +\infty\)¹.

Finally, one could consider smooth localized bump functions, for example, \[f_n^r(x) = \begin{cases} \frac{2^{r/2}}{C} \exp\left(\frac{1}{4(2^rx -n)^2-1}\right), & \text{if } |x| < 1/2, \\ 0, & \text{otherwise}, \end{cases}\] with \(C=\sqrt{\int_{-1/2}^{1/2} \mathrm{d}x \exp\left(2/(4x^2-1)\right)}\). These functions are local, differentiable, and orthonormal at fixed \(r\). However, they violate our last inclusion condition: a bump function at resolution \(r\) cannot be written as a linear combination of bump functions at resolution \(r+1\).

2.3 Daubechies scaling functions↩︎

Daubechies scaling functions verify the three conditions we demand, on top of other favorable properties. They are presented thoroughly in standard references [24], [43]. In what follows, we introduce them in a slightly different way, to show why they are natural for our purposes.

The starting point is a single function \(s\), which we assume supported on \([0,N-1]\) with \(N\) even. We will specify it only later, once we understand what it should verify. To construct our family of functions from this single bump (or scaling) function, we introduce the dyadic dilation operator \(\mathcal{D}\) and translation operator \(\mathcal{T}\) defined in the following way: \[\begin{align} \mathcal{D}^r s (x)=& 2^{r/2} s(2^r x) \text{ for } r \in \mathbb{Z} \\ \mathcal{T}^n s(x) =& s(x-n) \text{ for } n \in \mathbb{Z}. \end{align}\] Combining these operations, we write the \(n\) translated scaling functions at resolution \(r\) as \[s_n^r(x) := \mathcal{D}^r \mathcal{T}^n s(x) = 2^{r/2}s(2^r x - n).\] We now demand that the functions \(s_n^r(x)\) are orthonormal at fixed resolution \(r\), and write \(\mathcal{S}_r= \text{span}\{s_n^r(x)\} |_{n \in \mathbb{Z}}\). Then \(\mathcal{S}_r\) is the subspace of \(L^2(\mathbb{R})\) that the \(s_n^r(x)\) span.

For our inclusion condition to be verified, we further ask that \(s_n^r\) can be written as a linear combination of \(s_m^{r+1}\). Since the functions are compactly supported, this has to be a finite sum, called a refinement relation: \[\label{eq:WTs} s_n^r(x) = \sum_{i=0}^{N-1}h_{i} s_{2n+i}^{r+1}(x),\tag{3}\] where \(h_i\) are the weights or filter coefficients. This relation naturally implies that \[\label{eq:sub95inc} \mathcal{S}_0 \subset \mathcal{S}_1 \subset \dots \subset \mathcal{S}_r \subset \mathcal{S}_{r+1} \subset \dots \subset L^2(\mathbb{R}) \, .\tag{4}\] For 3 to be true, the bump function \(s\) has to be very special. Indeed, using the definition of the dyadic dilation operator, it is equivalent to the fixed-point equation \[\label{eq:RGs} s(x) = \sum_{i=0}^{N-1} h_i s_i^1(x) = \underset{\mathfrak{h}\cdot s(x)}{\underbrace{\sqrt{2} \sum_{i=0}^{N-1} h_i s(2x-i)}}\, .\tag{5}\] This is an infinite dimensional linear equation of the form \(s= \mathfrak{h} \cdot s\), and thus it has a non-trivial solution only if \(\mathfrak{h}\) has an eigenvalue \(1\).

Demanding that the fixed-point equation has a non-trivial solution, that the resulting solution is regular, and that the integer-translated solutions are orthogonal is not quite sufficient to fix the weights \(h_i\). Daubechies scaling functions are obtained uniquely by demanding on top that one can represent polynomials of degree \(N/2-1\) exactly², i.e. that there exist coefficients \(c_n^{(k)}\) such that \[\label{eq:polys} x^k = \sum_{n=-\infty}^{\infty} c^{(k)}_n s_n(x),\text{ for }0 \leq k \leq N/2-1 \, .\tag{6}\] This property implies that the regularity of \(s\) increases with \(N\), and for \(N \geq 6\) the function \(s\) has a continuous first derivative [44], [45]. The weights \(h_i\) are obtained by solving the non-linear algebraic equations implied by the orthonormality constraint, the fixed-point equation 5 , and the representation equation 6 . These equations can be solved exactly with radicals up to \(N=10\) (and otherwise computed numerically).

Once the weights are known, one can solve eq. 5 numerically to find \(s\). One simple strategy is to write eq. 5 for \(x\) set to integer \(n\), which gives a finite linear system since \(s\) is compactly supported. Once \(s(n)\) is known, one can plug the values in the right-hand side of 5 to get \(s\) at half integer values. Recursively, one gets exact values for all dyadic rationals, which are dense in \(\mathbb{R}\). As this fractal construction may suggest, the resulting function \(s\) takes a peculiar form, at least for \(N>2\) (\(N=2\) corresponds to Haar scaling functions). At \(N=6\), we can visualize the refinement relation 3

which is quite non-trivial! In practice, we need to compute integrals of products of \(s_n^r\), which can be obtained directly from the fixed-point equation. This is all we need to present our variational method at a fixed resolution.

3 MPS at fixed resolution↩︎

We now make our variational MPS manifold explicit, illustrate the minimization on the Lieb-Liniger model, and compute its energy density and correlation functions.

3.1 The variational manifold↩︎

The Hilbert space \(\mathcal{H}^r\) at fixed scale \(r\) is defined as the tensor product of the Fock spaces associated to the modes \[\label{eq:a95def} \hat{a}_n^r:= \int \mathrm{d}x \, s_n^r(x) \, \hat{\psi}(x).\,\tag{7}\] At this stage, each Fock space is still infinite dimensional. To use standard MPS algorithms we thus consider the truncated Fock space \(\mathcal{H}^r_d\subset \mathcal{H}^r\) where each mode is populated with at most \(d-1\) excitations. Since the particle density is finite in the ground state, we expect that when \(r\) is large, it is not necessary to increase \(d\) beyond \(3\) to get an arbitrarily good variational approximation.

For example, for \(I=[0, L]\) (\(L\) integer) and periodic boundary conditions, one can fit \(\ell = 2^r L\) such discrete truncated modes, and a basis of \(\mathcal{H}_d^r\) is \[\label{eq:numstate} \ket{m_0,...,m_{\ell-1}}_r := \frac{(\hat{a}_0^{r\dagger})^{m_0} (\hat{a}_1^{r\dagger})^{m_1} \dots (\hat{a}_{\ell-1}^{r\dagger})^{m_{\ell-1}}}{\sqrt{m_0! m_1! \dots m_{\ell-1}!}}\ket{0},\tag{8}\] where \(m_0,\dots, m_{\ell-1}\) are each in \(\{0,1,2,\dots,d-1\}\), and we have periodized the scaling functions.

Our final step is to consider an MPS manifold \(\mathcal{M}_\text{MPS}^{r,d,\chi}\) at finite bond dimension \(\chi\) within \(\mathcal{H}_d^{r}\). This manifold is made of the states of the form: \[\label{eq:MPS95def} \begin{align} \ket{A, r} &:= \!\!\sum_{m_0,\cdots, m_{\ell-1}=0}^{d-1} \!\!\!\!\!\!\!\!\!\!\text{Tr}(A^{m_0}_{0} ... A^{m_{\ell-1}}_{\ell-1}) \ket{m_0,... , m_{\ell-1}}_r \\ &= \quad \vcenter{\includestandalone[width=0.54\textwidth]{figures/MPSpbc}}\;, \end{align}\tag{9}\] where \(A_k^{m}\) are \(\ell\times d\) matrices of size \(\chi\times \chi\) and we have used the standard tensor network notation [46] on the second line.

In the translation invariant case, one can drop the \(k\) dependency, and there are just \(d\) matrices. They contain the variational parameters that one optimizes over to find the ground state. Further taking the thermodynamic limit, we obtain infinite MPSs \[\label{eq:iMPS95def} \ket{A,r} = \quad \vcenter{\includestandalone[width=0.54\textwidth]{figures/MPSgs} } \, ,\tag{10}\] which are a submanifold of the original Hilbert space \(\mathcal{H}\) as illustrated in Fig. 1. We call these states \(\ket{A,r}\) wavelet-MPS or wMPS³.

3.2 Evaluating the Hamiltonian↩︎

We now aim to evaluate the expectation value of operators, and in particular the Hamiltonian (density) on the states of our manifold. To this end, we first note that for any normal-ordered polynomial \(P\) in \(\hat{\psi}^\dagger\) and \(\hat{\psi}\), and a state \(\ket{\xi} \in \mathcal{H}^r\), we have \[\label{eq:normal95ordered95at95r} \bra{\xi} :\! P(\hat{\psi}^\dagger,\hat{\psi})\!: \ket{\xi} = \bra{\xi} :\!P(\hat{\psi}^\dagger_r,\hat{\psi}_r)\!: \ket{\xi},\tag{11}\] with \(\hat{\psi}_r\) the continuous mode operator projected at the resolution \(r\): \[\hat{\psi}_r(x) := \sum_{n\in\mathbb{Z}} s_n^r(x)\, \hat{a}_n^r \, . \label{eq:test}\tag{12}\] This is intuitive, and follows formally from the definition of the modes at scale \(r\) 7 and the canonical commutation relation \([\hat{\psi}(x),\hat{\psi}^\dagger(y)]=\delta(x-y)\).

This implies that we have \(\bra{\xi} H \ket{\xi} = \bra{\xi} H^r \ket{\xi}\) for \(\ket{\xi} \in \mathcal{H}^r\) with the Hamiltonian at resolution \(r\): \[\begin{align} H^r &=\int_I \mathrm{d}x \sum_{n,m}\hat{a}_n^{r\dagger} \hat{a}_m^r \partial_x s_n^r(x) \partial_x s_m^r(x) \\ &+ c \int_I \mathrm{d}x\sum_{n,m,l,k}\hat{a}_n^{r\dagger} \hat{a}_m^{r\dagger} \hat{a}_l^r \hat{a}_k^r s_n^r(x) s_m^r(x) s_l^r(x) s_k^r(x) \\ &- \mu \int_I \mathrm{d}x\sum_{n,m}\hat{a}_n^{r\dagger} \hat{a}_m^r s_n^r(x) s_m^r(x) \label{eq:fullH}\,. \end{align}\tag{13}\] Inverting the sum and the integral over \(x\) gives \[\begin{align} H^r&= 2^{2r}\sum_{n,m} K_{m-n} \hat{a}_n^{r\dagger} \hat{a}_m^r \\ &+ 2^r c \sum_{n,m,l,k} \Gamma^4_{m-n, l-n, k-n} \hat{a}_n^{r\dagger} \hat{a}_m^{r\dagger} \hat{a}_l^r \hat{a}_k^r \\ &- \mu \sum_n \hat{a}_n^{r\dagger} \hat{a}_n^r, \end{align}\] where the kinetic coefficients \(K\) are given by the integrals over the derivative of scaling functions \[\begin{align} \label{eq:K} \int \mathrm{d}x\, \partial_x s_n^r(x) \partial_x s_m^r(x) &= 2^{2r} \int \mathrm{d}x\, \partial_x s(x) \partial_x s_{m-n}(x) \nonumber \\ &:= 2^{2r} K_{m-n}, \end{align}\tag{14}\] and the four-point interaction coefficients \(\Gamma^4\) are defined from the integral of \(4\) scaling functions \[\label{eq:Gamma4} \begin{align} &\int \mathrm{d}x \; s_n^r(x) s_m^r(x) s_l^r(x) s_k^r(x) \\ &\quad =2^{r} \int \mathrm{d}x\; s(x) s_{m-n}(x) s_{l-n}(x) s_{k-n}(x)\\ &\quad := 2^r \Gamma^4_{m-n, l-n, k-n}. \end{align}\tag{15}\] These coefficients can be evaluated using the fixed-point equation 5, as detailed in the Appendix. Because \(s\) is compactly supported, both the kinetic and interaction coefficients are zero for site separations larger than \(N-1\) (for the interaction term this means the separation between any pair of sites must be less than \(N-1\) for the term to be non-zero). This implies that \(H^r\) is local, albeit with interactions beyond nearest neighbor.

In what follows, we consider Daubechies wavelets of order \(N=6,8\) and compute the coefficients \(K\) and \(\Gamma^4\) for each order once and for all. Since it is local, we can then express \(H^r\) as a Matrix Product Operator (MPO), or an infinite MPO [47] in the thermodynamic limit: \[\label{eq:Hr95mpo} H^r = \quad \vcenter{\includestandalone[width=0.54\textwidth]{figures/MPO}}\,.\tag{16}\] This MPO encoding is exact, which is crucial to keep our method variational. However, because the interactions are not nearest neighbor, especially in the 4-point interaction term, the MPO bond dimension is quite large at \(\chi_\text{MPO}=702\) for \(N=6\) and \(\chi_\text{MPO}=3306\) for \(N=8\), and would grow even larger for larger \(N\). Despite this large MPO bond dimension, calculations remain numerically tractable.

3.3 Energy density↩︎

We are now ready to minimize the energy density of the Lieb-Liniger model over our variational manifold at fixed resolution \(r\). This step is now completely standard, and we can use any off-the-shelf algorithm. In practice, we use MPSKit.jl [48] to run the variational uniform matrix product states algorithm (VUMPS) [49], [50], before switching to gradient descent on the Grassmann manifold [51] as the state nears convergence. This switch method is standard, and typically motivated by the fact that it is faster for certain benchmark models [51]. In our case, the primary motivation is that VUMPS, though fast, becomes unstable at high resolutions when nearing convergence. For \(r\geq 6\), energy minimization becomes difficult no matter the method, which motivates the introduction of iterative refinement in Sec. 4.

We compare the relative error in the ground state energy density we obtain with the exact Bethe-Ansatz solution [52], [53] in Fig. 2. We first consider a bond dimension \(\chi=16\), which is large enough that the error for \(r\leq 6\) is dominated by the finite lengthscale introduced. We observe that the relative error then decays exponentially in \(r\), i.e. linearly in the typical lengthscale \(2^{-r}\). Although this is the same decay one could find with more naive discretization methods, our method has the advantage of being truly variational. One may suspect that one could get a faster decay with smoother scaling functions. The \(D6\) scaling functions have Hölder exponent \(\alpha=1.0878\) [24], so they are just above the threshold of differentiability. The \(D8\) scaling functions have \(\alpha=1.6179\), and are thus substantially more regular (although not twice differentiable). As we can see in Fig. 2, the convergence rate seems to be unrelated to this regularity.

It might be worth trying even more regular functions, but we suspect that it is not possible to get faster convergence without abandoning our approach consisting in first defining an orthonormal family of \(L^2(\mathbb{R})\), and then building local Fock spaces from it. For example, in [40], the authors obtain a better error scaling by defining the \(n\)-body wavefunctions on small intervals separately (namely, states with two excitations are not made from tensor products of the one-excitation states present in the basis already). In that case, the Hilbert space they consider, analogous to our \(\mathcal{H}^r_d\), is not a truncated Fock space. This requires a subspace construction tailored to the Hamiltonian considered (and in principle, rebuilt for every coupling). Further, in the case of [40], it makes the method not variational since one needs to add a formally infinite penalty term to connect intervals.

Our solution to the slow convergence problem is different. Since the error decreases exponentially in \(r\), we can reach any precision desired provided we can increase \(r\) arbitrarily, through iterative refinements, without re-optimizing from a random state every time. This is the method we present in Sec. 4, and that leverages the full power of the multi-resolution analysis and makes use of wavelets.

We now compare our results with CMPS values obtained with CMPSKit.jl [54]. To this end, we fix \(\chi\) and push to much larger \(r\) (using the refinement method described in the next section). As expected the relative error in the energy density reaches a plateau, where finite entanglement is the bottleneck. Interestingly, this plateau is precisely the CMPS value. This leads us to conjecture that, in fact, our wMPS at \(r\rightarrow +\infty\) are exactly CMPS.

3.4 Correlation functions↩︎

Once the approximate ground state \(\ket{A,r}\) has been obtained, computing the expectation value of observables is straightforward. We insist that, just like with the Hamiltonian, the observables are not discretized or approximated in any way. In particular, unlike with a naive lattice discretization, we can compute correlation functions at arbitrary real separation \(x\), even if the scale \(r\) is fixed. To this end, we use eq. 11 to express the exact expectation value of \(\ket{A,r}\) as an expectation value involving only the scaling modes \(\hat{a}_n^r\). For example, for the two-point function, we have

Crucially, this is a finite sum, because \(s\) has compact support. The density-density correlator \(\bra{A, r}\hat{\psi}^\dagger(x)\hat{\psi}^\dagger(0)\hat{\psi}(x)\hat{\psi}(0)\ket{A, r}\) can be computed in an analogous way.

The two-point and density-density correlation functions are shown for different resolutions in Fig. 3. We deliberately chose small \(r\) here so that the finite \(r\) artifacts are visible. We get the expected behavior at sufficiently large \(r\). In particular, for the density-density correlator, our results become visually indistinguishable from the exact Tonks-Girardeau (TG) limit \(c\rightarrow +\infty\), with the Friedel oscillations visible already at smaller \(r\).

Interestingly, for the smallest \(r\), we see that the correlators are not just coarser versions of the higher \(r\) ones (with the rough pattern of the \(D6\) scaling function imprinted on top), but are actually globally shifted. This is expected of a true variational method: we are not just coarse graining the exact solution, we are minimizing over coarser states, and thus reaching slightly different approximate ground states for each \(r\).

4 Multi-resolution analysis↩︎

So far we have only considered wMPS at fixed resolution: each time the resolution is increased, we start from a random initial state, and run the full ground-state optimization routine. However, by construction \(\mathcal{H}^{r}\subset \mathcal{H}^{r+1}\), and thus, at least for some \(\tilde{\chi}\geq \chi\), we have \(\mathcal{M}_\text{MPS}^{r,d,\chi}\subset \mathcal{M}_\text{MPS}^{r+1,d,\tilde{\chi}}\). Thus, in analogy with bond-space expansion techniques [50], [55]–[58], we should be able to start from an approximate ground state at resolution \(r\), and promote it to an approximate ground state at resolution \(r+1\) with at most the same energy density. This refined approximate ground state could then be taken as a starting point for energy minimization, which would converge far quicker than when starting from a random state.

This refinement can be done using the full strength of multi-resolution analysis, which is what we do in this section. We first introduce the second half of the theory of Daubechies wavelets that we had postponed, then present our strategy to implement the refinement, explain in detail the different steps, and finally test it. Quite unexpectedly, we find that refining without optimizing further already substantially lowers the energy density.

4.1 Wavelet functions↩︎

Using the refinement relation 3 , one can write modes at a coarse scale in terms of modes at a finer scale, and thus naturally write a state in \(\mathcal{H}^r\) as a state in \(\mathcal{H}^{r+1}\). We will later need to implement this embedding locally, with a unitary circuit. When doing so, it is helpful to understand the orthogonal complement of \(\mathcal{H}^r\) within \(\mathcal{H}^{r+1}\), and this is when wavelet functions become helpful.

The wavelet space \(\mathcal{W}_r\) is simply defined as the orthogonal complement of \(\mathcal{S}_r\): \(\mathcal{S}_{r+1} = \mathcal{S}_r \oplus \mathcal{W}_r\). Applying this recursively, one covers \(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 \dots \; .\] The subspace \(\mathcal{W}_r\) is spanned by a new set of functions, \(\{w_n^r(x)\} |_{n\in \mathbb{Z}}\), the wavelet functions. They are orthogonal to the scaling functions, satisfy \(w_n^r(x) = 2^{r/2}w(2^r x - n)\), and can be written explicitly in terms of the scaling functions: \[\label{eq:WTw} w_n^r(x) = \sum_{i=0}^{N-1}g_{i} s_{2n+i}^{r+1}(x),\tag{17}\] where \(g_i := (-1)^i h_{N-1-i}\) are the wavelet function weights.

At resolution \(r+1\), the orthonormal families \(\{s_n^{r+1}(x)\} |_{n \in \mathbb{Z}}\) and \(\{s_n^r(x)\} |_{n \in \mathbb{Z}} \bigcup \{w_n^r(x)\} |_{n \in \mathbb{Z}}\) span the same subspace. The explicit relation is given by the inverse wavelet transform: \[\label{eq:IWT} \begin{align} s_{2n}^{r+1}(x) =& \sum_{i=0}^{\frac{N}{2}-1} h_{2i} s_{n-i}^r(x) + g_{2i} w_{n-i}^r(x) \\ s_{2n+1}^{r+1}(x) =& \sum_{i=0}^{\frac{N}{2}-1} h_{2i+1} s_{n-i}^r(x) + g_{2i+1} w_{n-i}^r(x). \end{align}\tag{18}\] The requirement that scaling functions can locally represent polynomials implies that the moments of the wavelet functions up to order \(N/2-1\) are zero: \[\int_\mathbb{R}\mathrm{d}x\, x^k w(x) = 0,\text{ for }0 \leq k \leq N/2-1.\] This can be easily seen from eq. 6 and the fact that \(\int \mathrm{d}x s_n(x)w(x) = 0\) for all \(n \in \mathbb{Z}\). Just like the scaling function \(s\), the wavelet function \(w\) is compactly supported on \([0,N-1]\), becomes more regular as \(N\) increases, and has continuous first derivative for \(N \geq 6\).

Moving from \(L^2(\mathbb{R})\) onto Fock spaces, one can introduce the wavelet modes \[\hat{b}_n^r= \int_\mathbb{R}\mathrm{d}x\, w_n^r(x) \hat{\psi}(x) \, ,\] which commute with \(\hat{a}_n^r\) and verify \[[\hat{b}^r_n, \hat{b}^{q\dagger}_m] = \delta_{rq}\delta_{nm},\,\,\,\,\,\, [\hat{b}^r_n, \hat{b}^q_m] = [\hat{b}^{r\dagger}_n, \hat{b}^{q\dagger}_m] = 0.\] As was done in [36], one can now express the continuous annihilation operator \(\hat{\psi}(x)\) exactly in terms of scaling and wavelet modes: \[\label{eq:DWT} \hat{\psi}(x) = \sum_n \left(\hat{a}_n^r s_n^r(x) + \sum_{j=0}^{\infty} \hat{b}_n^{r+j} w_n^{r+j}(x) \right)\, .\tag{19}\]

4.2 Strategy↩︎

We now give the general strategy of our refinement method, to motivate the subsequent technical subsections. For simplicity, we restrict to \(N=6\) Daubechies wavelets from now on.

We start from a matrix product state \(\ket{A, r}\) optimized at resolution \(r\). Our refinement algorithm is presented in Fig. 4 and involves three steps:

1. Applying the IWT 18 , written as a quantum circuit, to \(\ket{A, r}\) with the wavelet modes in the vacuum,

2. Contracting the resulting tensor network with \(3\) steps of the infinite time-evolving block decimation (ITEBD) algorithm [59],

3. Restoring translation invariance, to get back a wMPS state \(\ket{C, r+1}\).

After step (ii), the MPS has a two-site unit cell, but should be the same state as the original one, up to numerical truncation errors. However, the third step does change the state, and is not a pure change of gauge. Intuitively, since the ground state is translation invariant, forcing translation invariance at a finer scale could be helpful, and this is indeed what we observe.

4.3 The inverse wavelet transform as a quantum circuit↩︎

Writing the inverse wavelet transform as a quantum circuit is quite a standard step, that can be found in various forms in the literature (see e.g. [29]). We present it in detail here both for completeness, and also because we need a more explicit form. Namely, we need not only a formal expression at the mode level (which suffices for free theories), but the coefficients in Fock space.

The IWT is given in eq. 18 in terms of scaling and wavelet functions. It implies the following Bogoliubov transform at the level of modes \[\begin{align} \label{eq:IWTop} \hat{a}_{2n}^{r+1} =& h_0 \hat{a}_{n}^r& + h_2 \hat{a}_{n-1}^r + h_4 \hat{a}_{n-2}^r \nonumber \\ &&+ g_0 \hat{b}_{n}^r + g_2 \hat{b}_{n-1}^r + g_4 \hat{b}_{n-2}^r \nonumber \\ \hat{a}_{2n+1}^{r+1} =& h_1 \hat{a}_{n}^r& + h_3 \hat{a}_{n-1}^r + h_5 \hat{a}_{n-2}^r \nonumber \\ &&+ g_1 \hat{b}_{n}^r + g_3 \hat{b}_{n-1}^r + g_5 \hat{b}_{n-2}^r\, , \end{align}\tag{20}\] or in abstract form \[\left(\begin{array}{c} \cdots \\ \hat{a}_{2n}^{r+1}\\ \hat{a}_{2n+1}^{r+1}\\ \cdots \end{array}\right) = \mathcal{U} \left(\begin{array}{c} \cdots\\ \hat{a}_{n}^{r}\\ \hat{b}_{n}^{r}\\ \cdots \end{array}\right)\,.\] The first step is to write this unitary transformation \(\mathcal{U}\) in mode space as a local brick-wall circuit. For \(N=6\), it has depth \(3\): \[\label{eq:bogo95as95brickwall} \mathcal{U}= \vcenter{\includestandalone[width=0.65\textwidth]{figures/brickwall}}\, .\tag{21}\] The three unitary two-by-two matrices—\(u^{(i)}\) for \(i=1,2,3\)—acting in operator space, can be determined from eq. 20 and one gets: \[\label{eq:ai} \begin{align} u^{(1)} =& \frac{1}{\sqrt{h_1^2+h_4^2}}\begin{pmatrix} h_1 & -h_4 \\ h_4 & h_1 \end{pmatrix}, \\ u^{(2)} =& -\begin{pmatrix} \sqrt{h_2^2+h_3^2} & -\sqrt{1-h_2^2-h_3^2} \\ \sqrt{1-h_2^2-h_3^2} & \sqrt{h_2^2+h_3^2} \end{pmatrix}, \\ u^{(3)} =& \frac{1}{\sqrt{h_4^2+h_5^2}}\begin{pmatrix} h_4 & h_5 \\ h_5 & -h_4 \end{pmatrix}. \end{align}\tag{22}\] So far the circuit 21 is a map acting on mode operators (or equivalently, on single particle wavefunctions). We need to promote it to a circuit acting on the whole (truncated) Fock space. This is conceptually straightforward even if practically tedious. To this end, consider a number state of two modes \(a,b\) and act with a \(2\times 2\) unitary \(u\) at the mode level: \[\begin{align} \ket{n_a,n_b} &:= \frac{\hat{a}^{\dagger n_a} \hat{b}^{\dagger n_b}}{\sqrt{n_a!n_b!}}\ket{0}\\ &\underset{u^{(i)}}{\longrightarrow} \frac{(u_{11} a + u_{21} b)^{n_a} (u_{12} a + u_{22} b)^{n_b} }{\sqrt{n_a!n_b!}}\ket{0} \label{eq:mode95transformed95factorized}\\ &:= U^{(i)} \ket{n_a,n_b}\, . \end{align}\tag{23}\] where \(u_{nm}\) are the matrix entries of \(u\). This defines a unitary operator \(U^{(i)}\) acting directly on Fock space, and that has the same effect as applying \(u^{(i)}\) at the mode level. By expanding the right-hand side of 23 , we find the rows of \(U^{(i)}\) one by one:

We then get the brick-wall quantum circuit advertised by promoting \(u^{(i)}\) to \(U^{(i)}\) in 21 .

At the end of this procedure, we have simply rewritten the state \(\ket{A,r}\in\mathcal{H}^r\) as a tensor network adapted to the tensor product structure of \(\mathcal{H}^{r+1}\): \[\label{eq:iwt95equality} \vcenter{\includestandalone[width=0.9\textwidth]{figures/iwt_equality}}\, ,\tag{24}\] where we use the sign \(\cong\) to mean that these different tensor networks, with legs living in different local Hilbert spaces, represent the same state in the global Hilbert space \(\mathcal{H}^r\subset\mathcal{H}^{r+1}\).

As discussed in depth in [27]–[31], a circuit that introduces wavelet modes and then performs a wavelet transform via two-site unitaries is an entanglement renormalization circuit. However, we do not stack several such layers at once, and thus depart from standard discussions connecting MERA and wavelets.

4.4 Refinement map contraction via ITEBD↩︎

We now need to contract the right-hand side of 24 to bring it back to MPS form. We do so with \(3\) steps of the ITEBD algorithm, contracting the brick-wall layers one at a time as summarized in Fig. 5.

For the first layer, at each site we contract the MPS tensor \(A\), the wavelet state, and \(U^{(1)}\). We carry a singular value decomposition (SVD) on the resulting tensor

to get an infinite two-site MPS with tensors \(U_1\) and \(s_1V_1\). The second layer is similarly contracted and decomposed as

and the process for the final layer follows by analogy. For the last two layers, the contraction includes the singular values \(s_{i-1}\), and the SVD yields \(\tilde{U}_is_i\tilde{V}_i\) with the \(s_{i-1}\) singular values absorbed into \(\tilde{U}_i,\tilde{V}_i\). To obtain the final diagram in the figure above, one takes inverses of the singular values such that \(U_i = s_{i-1}^{-1}\tilde{U}_i\) and \(V_i = \tilde{V}_i s_{i-1}^{-1}\).

The resulting state is an MPS with a 2-site unit cell at resolution \(r+1\), with tensors \(A\) and \(B\) (for example, we take \(A= s_3 U_3\) and \(B=s_2 V_3\)), and much larger bond dimension. In practice, we truncate the bond dimension back to \(\chi\) at the end, and observe that this only minimally impacts the energy density of the state.

4.5 Approximate projection onto translation invariant MPS↩︎

If we iterate the steps (i) and (ii) that we presented before \(M\) times, the unit cell will have size \(2^M\). Hence, before refining further, we need to get back to an MPS with a \(1\)-site unit cell. We should do so in a way that does not degrade the energy density. Fortunately, we expect the ground state to be translation invariant, and thus we expect the two-site MPS we obtain to be approximately translation invariant.

If the state were exactly translation invariant, then it would imply that there exists some unitary \(U\) for which \(U^\dagger A = BU = C\), explicitly restoring translation invariance. Note that translation invariance would automatically be restored by subsequent optimization (e.g. with VUMPS). However, by approximately projecting onto translation invariant MPS directly, we get two benefits: first, VUMPS is cheaper when run on single-site MPS, and second, we see that we have a discrete freedom in the projection that allows us to reach a substantially lower energy.

To find \(U\), we first put \(A,B\) in right-orthonormal gauge, to get \(A_R,B_R\), respectively. We then construct the following transfer matrix [60]:

\(C_R=U^\dagger A_R = B_R U\) which is possible only if the state is exactly translation invariant (in which case, making translation invariance explicit, and finding \(C\), is a pure gauge choice). Importantly, \(U\) is a fixed point of this transfer matrix with eigenvalue 1:

This allows one to find \(U\), then \(C_R\), and thus approximate the refined two-site MPS as a single-site state \(\ket{C_R, r+1}\). We say approximate, because in practice, \(U^\dagger A_R\) and \(B_R U\) are not exactly equal, so we are free to choose whichever one gives a state with lower energy density. Quite surprisingly, we observe that one choice systematically lowers the energy compared to the original two-site state, almost bringing it to the energy of the fully optimized state at resolution \(r+1\)!

This is the final step of the refinement algorithm. We now have a process which starts at an initial single-site MPS \(\ket{A, r}\) at resolution \(r\) with energy density \(e_0(r)\), and outputs a single-site state \(\ket{\tilde{A}, r+1}\) with lower or equal energy density.

4.6 Numerical results↩︎

We first look at the behavior of the energy density at each of the scale refinement steps in Fig. 6. As expected, the energy stays constant as we contract the brick-wall circuit and our resolution \(r\) MPS to get a two-site MPS at resolution \(r+1\). Then, quite surprisingly, approximately projecting onto a one-site MPS drastically lowers the energy density – it even gets closer and closer to its post-optimization value at \(r+1\) when \(r\) increases. Adding a few steps of an optimization algorithm at the end of the refinement gets us to the best possible energy at \(r+1\).

This iterative refinement approach is in fact the only way we manage to find approximate ground states in a reasonable amount of time for \(r\geq 6\).

We consider that we have sufficiently refined in scale when the energy density error is dominated by the finite \(\chi\) error, and not the scale error. On this toy model, we can distinguish the two cleanly by computing the CMPS energy density: when \(e_0^\text{wMPS}(r) - e_0^\text{CMPS} \ll e_0^\text{CMPS} - e_0^\text{exact}\) we are at a scale that is small enough. In Fig. 7, we show that without iterative refinement, optimizing at \(r=6\) directly, using VUMPS (on the edge of stability) followed by Grassmann gradient descent, the optimization does not converge (staying far from what iterative refinement reaches at the same \(r\)) and stays quite far from the threshold where the finite bond dimension error dominates. In contrast, when we alternate refinements and Grassmann gradient descent, we can reach much higher \(r\) and go to a regime where scale error is negligible.

5 Discussion↩︎

We have presented a new method to solve interacting continuum models with Daubechies scaling functions and MPS. The main advantage of the approach is that it is strictly variational: it gives certified upper bounds on the energy, and a clean way to compare the quality of different ground-state approximations. Using an iterative refinement scheme, which exploits multi-resolution analysis, it is possible to get accurate results on the example of the Lieb-Liniger model. In that case, for moderate bond dimension, the residual error is dominated by finite entanglement effects, and not the way one deals with the continuum.

Importantly, this method is compatible with the entirety of the discrete tensor network toolbox, as it is already implemented in standard libraries. While we demonstrated the power of the method on a translation invariant ground state problem with uniform MPS – because it makes comparison with exact results easier – it can be used on problems with open boundary conditions, for simulations of real-time dynamics, or to compute excitation spectra, simply by plugging in the corresponding MPS routine. This flexibility is the main advantage of wMPS over dedicated continuous algorithms like CMPS. The latter may however keep an edge in the thermodynamic limit when high accuracy and large bond dimensions are required.

The link between wMPS and CMPS remains to be fully elucidated. We did observe that a wMPS at \(r\rightarrow + \infty\) gives the same energy density as a CMPS of identical bond dimension. This naturally raises the question of how well one can discretize a CMPS into wMPS, which would allow one to conveniently change representation when needed. Since the theory of CMPS can be mapped to the theory of continuous measurement [61], this problem is related to optimal ways of binning time-continuous measurements into discrete records [62], [63], which is particularly relevant for the manipulation of superconducting qubits.

A natural next step would be to consider broader classes of wavelets beyond the Daubechies family. For example, relaxing orthonormality, with biorthogonal wavelets (as used e.g. in [31]), would allow the use of more symmetric and regular wavelets with the same support. More generally, it would be interesting to understand what criterion one should optimize to get the “best” function with fixed support for typical variational problems.

Another extension would be to consider relativistic QFT, where extra divergences (related to short-distance scale invariance) bring additional difficulties. In this context, relativistic CMPS (RCMPS) solve the problem, and are truly variational, but at the cost of extra non-locality, which translates into a higher numerical cost and the need for an ad hoc optimizer. If it is possible to use standard MPS algorithms while remaining variational in the relativistic case, this would be particularly useful.

Finally, there is in principle no obstruction to generalizing our wMPS to higher dimensions, with wavelet Projected Entangled Pair States (wPEPS). In practice, the bond dimension of the MPO in one dimension is already large, and so we anticipate it will be expensive to perform a numerical analysis with a wPEPS. In this context, however, there is less competition with continuous methods. In particular, despite some recent progress with continuous tensor network states [21], [22], [64]–[66], which can now be contracted in some non-trivial cases [67], there is no functional tensor network algorithm to treat the 2d Bose gas directly in the continuum.

6 Evaluating scaling function integrals↩︎

We detail the calculation of \(K_{a}\) in eq. 14 and \(\Gamma^4_{a,b,c}\) in eq. 15 , for any \(N\). As described in [68], the key to evaluating these coefficients without performing any integration is the fixed-point equation 5 . Starting with the kinetic coefficient, with \(-N+2 \leq a \leq N-2\), we use the fixed-point equation to find \[\begin{align} K_a &=& 2 \sum_{m,n=0}^{N-1} h_m h_n \int \mathrm{d}x\,\partial_x s(2x -m) \partial_x s(2x-2a -n) \nonumber \\ &=& 4 \sum_{m,n=0}^{N-1} h_m h_n \int \mathrm{d}x\,\partial_x s(x) \partial_x s(x-2a -n + m) \nonumber \\ &=& 4 \sum_{m,n=0}^{N-1} h_m h_n K_{2a+n-m} \nonumber \\ &=& 4 \sum_{n'=-N+2}^{N-2} \sum_{m=0}^{N-1} h_m h_{n'+m-2a} K_{n'} \nonumber \\ &:=& \sum_{n'=-N+2}^{N-2} \mathcal{A}_{an'} K_{n'} \end{align}\] where we defined \(\mathcal{A}_{an'} = 4 \sum_{m=0}^{N-1} h_m h_{n'+m-2a}\). This then makes it clear that \(K_a\) is an eigenvector of \(\mathcal{A}\) with eigenvalue 1. It is normalized by the condition \[\sum_{a=-N+2}^{N-2} a^2K_a=-2\] which is derived in [68] from the property that \(s_n(x)\) can locally represent polynomials up to degree \(N/2-1\).

The four-point coefficient in eq. 15 can be calculated in a similar fashion. Via the fixed-point equation 5 , \(\Gamma^4\) becomes

where we defined \[\mathcal{B}_{n'l'k'}^{abc}= 2 \sum_{m=0}^{N-1} h_m h_{n'+m-2a} h_{l'+m-2b} h_{k'+m-2c}.\] Then \(\Gamma^4_{a,b,c}\) is a fixed point of the tensor \(\mathcal{B}\), with eigenvalue 1. Again, by the polynomial representation property, we normalize \(\Gamma^4\) using \[\sum_{a,b=-N+2}^{N-2} \Gamma^4_{a,b,c} =\delta_{c0}.\] In this way, we find the kinetic and interaction coefficients for \(N=6,8\).

References↩︎