Non-Gaussian Entanglement Renormalization for Quantum Fields

In this work, a non-Gaussian cMERA tensor network for interacting quantum field theories (icMERA) is presented. This consists of a continuous tensor network circuit in which the generator of the entanglement renormalization of the wavefunction is nonperturbatively extended with nonquadratic variational terms. The icMERA circuit nonperturbatively implements a set of scale dependent nonlinear transformations on the fields of the theory, which suppose a generalization of the scale dependent linear transformations induced by the Gaussian cMERA circuit. Here we present these transformations for the case of self-interacting scalar and fermionic field theories. Finally, the icMERA tensor network is fully optimized for the $\lambda \phi^4$ theory in $(1+1)$ dimensions. This allows us to evaluate, nonperturbatively, the connected parts of the two- and four-point correlation functions. Our results show that icMERA wavefunctionals encode proper non-Gaussian correlations of the theory, thus providing a new variational tool to study phenomena related with strongly interacting field theories.


Introduction and Summary
The multiscale entanglement renormalization ansatz (MERA) [1,2], which was originally proposed as a variational method to obtain the ground state of spin chains systems, consists of a real space renormalization group technique that, iteratively, removes the quantum correlations between small adjacent regions of space at each length scale. A continuous version of MERA (cMERA) was proposed for free field theories [3,4]. Motivated, among others, by the conjecture that cMERA is a realization of the AdS/CFT correspondence [5][6][7][8][9][10], a rigorous and (non)perturbative formalism for interacting theories turns out to be essential to advance in this program.
Precisely, one of the major problems in quantum field theory (QFT) is the understanding of the phenomena associated to strongly coupled systems. To do so, nonperturbative methods are required. While these are difficult problems to solve exactly, it is acknowledged that the use of nonperturbative variational methods allow to tackle these problems to some extent. Despite field theory was initially formulated within the Hamiltonian framework, these methods lost relevance against path integral techniques for many years.
Nevertheless, when it comes to nonperturbative aspects, there are situations in which the use of wavefunctionals on configuration space exhibits clear advantages. Namely, nonperturbative path integral methods are especially suited to compute quantities that have no perturbative contributions and can be addressed through a saddle point approximation.
However, when the observables of interest receive both perturbative and nonperturbative contributions, dealing with the path integral becomes more difficult. In addition, with the raising appeal in understanding QFT from a quantum information point of view, the Hamiltonian framework involving wavefunctionals seems more suitable than other approaches.
In recent years, tensor networks, a new class of variational states, have proven to be very useful in the study of a huge variety of interacting many body systems. Initially devised for lattice systems, through the Rayleigh-Ritz variational principle, a tensor network representation of the wavefunction provides an efficient approximation to the ground state of an interacting many body system by systematically identifying the relevant degrees of freedom for the physics at low energies. As an instance, the MERA tensor network [1,2], implements a variational real space renormalization group on the wavefunction that represents the ground state of the system at different length scales. In spite of their success in analyzing 1D systems on the lattice, several difficulties arise when trying to generalize tensor networks to higher dimensions and/or to interacting field theories. In this sense, the continuous generalization of the the matrix product state (cMPS), while proving efficient to describe the low energy physics of nonrelativistic systems in 1+1 dimensions, suffers from regularization ambiguities when dealing with relativistic systems in 1+1 di-mensions [11].
A continuous tensor network circuit designed to work in arbitrary spatial dimensions for (non)relativistic field theories is the continuous version of MERA, cMERA [3,4]. This tensor network builds a multi-layered representation of the ground state wavefunctional through a variationally optimized pattern of entanglement between the relevant degrees of freedom at any length scale. cMERA amounts to a real space renormalization group (RG) of the wavefunctionals in a Hamiltonian framework such that, each layer of the network corresponds to a step of the RG flow. Regrettably, cMERA has only been explicitly formulated for free theories of bosonic, fermionic and gauge fields [3,4,12]. In these formulations, the cMERA renormalization in scale is generated by a quadratic operator, and thus, the resulting state is given by a Gaussian wavefunctional. Obviously, while this fact dramatically limits the interest of the Gaussian cMERA ansatz for interacting QFTs, such trial states result useful to reproduce correlation functions and entanglement entropy in free field theories [13]. However, to clarify whether cMERA is a possible realization of the holographic duality [5][6][7][8][9][10], a more general formulation to study nonperturbative interacting field theories becomes crucial.
Our aim in this work is to develop a truly non-Gaussian cMERA tensor network circuit able to nonperturbatively capture relevant phenomena of interacting QFTs. In this respect, let us comment on some aspects that stem from applying a variational method (e.g., tensor networks) to QFT: generality, calculability and ultraviolet modes.
Firstly, one has the problem of the generality of the trial state. Namely, the trial state should be general enough to capture the most salient physical features of the phenomena under study through the variation of its parameters. Due to the enormous size of the Hilbert space in a QFT, it is very difficult to identify by mere intuition the relevant parameters that have to be probed. Thus, a systematic method to build ansatze is desirable.
Secondly, there is the problem of calculability [14]. That is to say, even when possessing a reasonable and flexible ansatz for the vacuum wavefunctional, one needs to evaluate expectation values of operators/observables of interest in this state, This amounts to the evaluation of an Euclidean functional integral in which the square of the wavefunctional acts as the partition function. Thus, in QFT, given the very limited ability to evaluate non-Gaussian path integrals, the calculability requirement on the trial wave functional is certainly severe. Indeed, it is so severe that it has constrained the form of the trial wavefunctionals to Gaussian states. Despite this, Gaussian trial states have been successfully applied to self-interacting relativistic scalar and spinor field theories [15], where a great amount of nontrivial exact results in the large N limit has been obtained (among others, a proposal to build a Gaussian approach to cMERA for interacting field theories [16]). Namely, the Gaussian ansatz (that is the exact ground state in noninteracting QFTs) works very well for settings in which the relevant nonperturbative physics of the system is dominated by a single condensate. However, it is well known that the connected part of N -point correlation functions distinguishes the ground states of interacting theories from those of noninteracting ones. While for Gaussian states, the connected correlation functions of order higher than two vanish, those of interacting systems are generally nonzero. For this reason, despite being successful in capturing some nonperturbative effects, with the aim of going beyond the Gaussian approach, it would be desirable to develop in a systematic way, ansatze which variationally borne in, some kind of "generalized" condensates while keeping the calculability of the Gaussian ansatz intact.
Finally, but not least important, one is faced to the problem of the ultraviolet modes.
The main objective of a variational calculation in a strongly interacting field theory is to obtain the correct configuration for the low momentum modes of the field in the vacuum wavefunctional. Due to the interaction between the high and low momentum modes in an interacting QFT, it is thus desirable to have a method that yields variational parameters that optimally integrate out the effects of high energy modes into the low energy physics.
Bearing in mind these features, various proposals have been recently raised to go beyond a purely Gaussian wavefunctional. In [17,18], authors have developed techniques to carry out systematic perturbative calculations of cMERA circuits restricted to a weakly interacting regime. Another recent approach, [19], proposes a particular realization of a Gaussian cMERA with an UV structure analogous to that of the cMPS. It is expected that this connection could yield cMERA wavefunctionals that are able to capture nonperturbative physics beyond the Gaussian approximation. In [20], authors presented a truly nonperturbative method to build non-Gaussian cMERA wavefunctionals for interacting QFTs. This approach relies on nonlinear canonical transformations (NLCT) [21][22][23][24] to build a set of scale dependent extensive wavefunctionals which are certainly non-Gaussian.
With this prescription, we showed that observables such as disconnected correlation functions can be analytically calculated in a closed form.
Let us remark the last proposal. The scale dependent NLCT in [20], which were constructed from the product of two unitary operators, obviously admit a realization through a unitary operator. This is what could constitute a more satisfactory and natural continuum generalization of the lattice MERA algorithms. In the present work we refine our method and solve this issue using the firm ground provided by the conceptual framework presented in [20]. In addition, we carry out an exhaustive optimization process and show the efficiency and predictability of our method.
Let us comment on the structure of the present work and summarize the main results of each section. In Section 2 we have reviewed and generalized the cMERA formalism for any entangler beyond the Gaussian one. Using a quantum mechanical toy model, one can prove the following result (see (23), [17,18]): [order 2, order k] ≤ order k .
That is to say, only when commuting with order-2 operators, an order-k operator does not increase the order. On the other hand, entanglers containing order-k operators, k ≥ 3, are necessary to generate non-Gaussian states. However, the Hadamard's lemma and the above result show that they are very problematic when acting on the field operators of the theory. Thus, if the approach must be nonperturbative, an alternative solution appears to be necessary to solve this puzzle.
In Section 3 we have reviewed the quadratic Gaussian entangler introduced in [3], which performs scale dependent Bogoliubov transformations on the fields. Based on (2), we can grasp why both free and interacting theories can be treated with Gaussian trial wavefunctions: when applying the Gaussian entangler, the order of the renormalized operators does not grow. Afterward, upon studying the free and the λφ 4 scalar theory, we give evidence of the limitations of the Gaussian entangler to study interacting theories.
For example, we introduce the connected N -point correlators and justify why the 4-point correlators or the kurtosis automatically vanish for Gaussian trial wavefunctions.
In Section 4 we have firstly studied the nonquadratic operator (50) Naturally, F will play the role of a fermionic entangler when scale dependence is included.
In Section 5, we have considered a circuit formed by the exponential of the entangler where u parametrizes the scale of the renormalization and P is the u-ordering operator. Here L is the dilatation operator and the generating operator K(u) is the so-called "entangler". The renormalization scale parameter u is usually taken to be in the inter- is the scale at the UV cutoff , and the corresponding momentum space UV cutoff is Λ = 1/ . u IR is the scale in the IR limit.
The state |Ψ Λ ≡ |Ψ U V is the state in the UV limit and it may be the ground state of a quantum field theory. The state |Ω is such that there is no entanglement between spatial regions upon which the cMERA flow builds correlations at successively smaller distance scales. We will impose |Ω to be scale invariant with respect to spatial dilatations, so that e −iLu |Ω = |Ω or, equivalently, In this work, we will assume |Ω to satisfy the following condition 1 for all momenta p, where ω Λ = √ Λ 2 + m 2 with m the mass of the particles in the free theory and χ 0 = Ψ Λ |φ(x)|Ψ Λ . This state satisfies The nonrelativistic dilatation operator L does not depend on the scale u and it is only governed by the scaling dimensions of the fields. It is taken as the "free" piece of the cMERA Hamiltonian and is given by The transformation properties of the field under the scale transformation L are given by: The entangler K(u), which contains all the variational parameters to be optimized, creates entanglement between field modes with momenta |p|< Λ, where Λ is a generic cutoff. Actually, as we will see in Section 5, various cutoff parameters can co-exist simultaneously, in such a way that they regulate the strength of each of the variational parameters at different regions of momentum space. As a consequence, the entangler is considered the "interacting" part of the cMERA Hamiltonian. In our approach, only the choice of K(u) will fully determine whether |Ω transforms into a Gaussian or a non- can be understood as a Hamiltonian evolution with K(u) + L along the scaling parameter u. As such, it is useful to define cMERA in the "interaction picture" through the unitary transformation of states In this picture, the entangler is given bŷ and the u-evolution is determined by the unitary operator Thus, one may write the wavefunctional in the interaction picture as

Non-Gaussian States: Beyond Quadratic Entanglers
As the ground states of interacting theories are generically non-Gaussian, our approach to construct cMERA for interacting theories consists of considering scale dependent non-Gaussian trial wavefunctionals. Such states will be generated by the action of an entangler that contains nonquadratic operators on a Gaussian state. Precisely, this formulation is in correspondence with the fact that in QFT, trial states created by introducing polynomial corrections to a Gaussian state represent a finite number of particles and those are suppressed in the thermodynamic limit. As a consequence, in going beyond the Gaussian ansatz, we are forced to use a class of variational extensive states for which the energy density does not depend on the volume. Regrettably, commutators of nonquadratic operators yield operators with increasingly larger products of φ's and π's [17] and hence, do not close as an algebra. This is a significant obstacle to systematically define unitaries that build non-Gaussian states from a Gaussian or another non-Gaussian state, which is precisely the crucial step to define cMERA flows for interacting field theories.
Without loss of generality, we will translate this problem to (0+1)-dimensional quantum mechanics. Our conclusions can be straightforwardly translated to (d + 1)-dimensional field theory. In this case, we introduce the Bender-Dunne algebra of Hermitian operators Let us consider the subalgebra {T m,n } +∞ m,n=0 and study the unitaries of the form The order of an operator is defined as the largest value of m + n, for which c m,n = 0.
Motivated by cMERA, we are interested in calculating quantities like ψ|U † T m,n U |ψ , where T m,n would be the Hamiltonian terms. Then we will have ψ|U † T m,n U |ψ = ψ|exp −i ad ∞ p,q=0 cp,qTp,q T m,n |ψ , where we have used the Hadamard's lemma Interestingly, one can show that the commutator of an order-2 and an order-k operator gives rise to, at least, an order-k operator: [order 2, order k] ≤ order k .
In particular, for k = 2, operators form the Lie algebra Heis(3, R) ⊕ sl(2, R), and generate unitaries that map Gaussian states to Gaussian states. Namely, if Q, Q are order ≤ 2, then the commutator [Q, Q ] will also have order ≤ 2. As a result, if |Ψ is a Gaussian state, then e −iQ |Ψ is also a Gaussian state. In this sense, we note that the entangler operator in Gaussian cMERA [1] can be seen as the operator K 2 with a scale dependent variationally optimized coefficient c (1) 2 (p 1 , p 2 ; u). Actually, (23) ensures that even interacting terms, as for example λφ 4 , remain under control when transformed with unitaries containing order-2 entanglers, as it was done in [16]. However, because trial states are still Gaussian, non-Gaussian effects will be absent. For example, quantities as kurtosis or connected 4-point correlators automatically vanish [20].
When entanglers containing higher order operators are considered, the situation is less trivial: If the order of the entangler is higher than 2, then the quantity U † T m,n U is an operator of infinite order. Namely, according to (23), the nested commutators induce higher-order operators at every step which propagate out of control. This is why order-2 entanglers are the only ones that can be straightforwardly used on any free or interacting theory.
To get over this issue, various methods have been proposed. For example, a method based on the following unitary where is assumed to be a small parameter, proves to be useful for a perturbative treatment. In [17] it is shown how to perturbately obtain the ground state of the quantum mechanical toy model of the anharmonic oscillator H = p 2 /2m + mx 2 /2 + λx 4 from a Gaussian state. In particular, upon Taylor expanding the unitary (24) above = 0, identifying = λ and choosing a set of T p,q operators with p + q > 2, the non-Gaussian ground state |ψ anharm is given by where |ψ harm is the ground state of the harmonic oscillator. Based on these principles, several entanglers have been proposed in field theory at a perturbative level [17,27].
In [20], we have proposed a different approach to truncate the infinite series of higher order operators. Based on [21,22], we consider a cMERA entangler in terms of the following field theory anti-Hermitian operator: where m ≥ 2 is a parameter to be fixed. 2 Here, s is the variational parameter whereas g(p, q 1 , . . . , q m ; u) consists of a combination of scale dependent variational cutoff parameters. The explicit dependence on the variational cutoff parameters will certainly be essential for the truncation of the series of nested commutators and for the optimization.
In contrast to quadratic entanglers, this operator clearly goes beyond Bogoliubov transformations and induces nonlinear canonical transformations on the fields φ and π. It is important to emphasize that the variational parameter s is not necessarily perturbative 3 . In addition to this, in Section 5 we will show that, when applied, scale dependent non-Gaussian effects are unambiguously captured: the connected part of multi-point correlators, kurtosis or skewness clearly result scale dependent nonvanishing quantities [20].
In summary, upon reviewing the cMERA formalism in full generality, we have reduced the problem of having (non-)Gaussian trial states to a particular choice of the operators entering the entangler. In the following sections we will review the Gaussian cMERA formalism and explain the construction of a scale dependent non-Gaussian cMERA circuit.
Finally, we will perform the optimization on the λφ 4 theory and study various truly non-Gaussian observables.

Quadratic Entangler
For free scalar theories in (d + 1) dimensions with d, the spatial dimensions of the theory,

K(u) is given by the quadratic operator [3, 4]
where (27) is the only variational parameter to optimize in the cMERA circuit.
This function factorizes as where Γ(x) ≡ Θ(1 − |x|) and Θ(x) is the Heaviside step function; g(u) is a real-valued function known as density of entanglers and Γ(p/Λ) implements a high frequency cutoff such that p ≡ Λ p [3,4]. The sharp cutoff function ensures that K(u) acts locally in a region of size Λ −1 . It is also possible to define the entangler K(u) through localized smooth smearing functions instead of sharp cutoff functions [19].
In the interaction picture the entangler operator reads as, where the integral in the second line will be suppressed by the cutoff for |k|≤ Λe u .
With this, the operator U G (0, u) ≡ Pe −i u 0 du (K(u )+L) defines the cMERA evolution in terms of the scale-dependent linear transformation of the fields We have used the subscript G because it is straightforward to show that the cMERA wavefunctional (6) with the quadratic entangler (27) and a reference state such as (8) can be written as the Gaussian wavefunctional given by where N is a normalization constant and the scale dependent Gaussian kernel F (p; u) is defined through the variational cMERA parameter g(p; u) by [16] with F (p; 0) = (p 2 + m 2 ) −1/2 . We note that this wavefunctional is built as where the operator that shifts the argument of any functional (and specifically the Gaussian wavefunctional) by a constant χ 0 , is given by This is an standard type of transformation which leaves invariant the measure of integration in the functional path-integral formalism invariant, i.e., Dφ = Dφ where φ = φ − χ 0 , the shift of the variable of integration by a fixed background field configuration. A more general possibility, which will be analyzed in the next section, was introduced in [21 -23] and amounts to shifting part of the field modes φ(k) by a nonlinear polynomial functional of other field modes. In geometrical terms this is tantamount to nonlinearly deform the infinite-dimensional configuration vectors {|φ(x) } of the field theory. For the case of free theories, generating the set of "fixed-background" shifted scale dependent cMERA Gaussian wavefunctionals (31) is enough from a variational point of view.

Free Scalar Theory
In order to "solve" the cMERA ansatz for the relativistic free massive scalar theory, one applies the variational principle to minimize the energy functional given by with ω k = √ k 2 + m 2 [3,4]. The total energy E is given by Then, we impose which yields Using that one obtains and

Self-Interacting Scalar λφ 4 Theory
As commented in Section 1 and based on the discussion of Section 2.1, the Gaussian ansatz has been widely used in the context of interacting field theories. In this sense, in [16], a cMERA circuit based on the quadratic entangler (27) with respect to the ansatz wavefunctionals Ψ[φ, u] of the form (31). As shown in [16], to solve the variational problem we compute . Then we impose δE/δF (p) = 0, which yields an optimized Gaussian kernel F (p), where M is the modified mass of the propagating Gaussian quasi-particles Finally, through (32), we obtain the optimized Gaussian cMERA variational parameter Remarkably, the wavefunctional optimization over an infinite-dimensional space of kernels F has reduced to solving a single nonlinear equation for M . The Gaussian cMERA wavefunctional thus obtained is a vacuum state for a free theory with mass given by (44).
The optimized ansatz captures all 1-loop 2-point correlation functions, and additionally captures the resummation of all cactus-like diagrams [28]. Nevertheless, the ansatz is unable to capture the effects of the interaction on higher order correlation functions. Let us extend a bit on this point.
Quantum field theories are characterized by correlation functions of the form where O is a field operator, N is the order of the correlation and, in general, the expectation value is taken with respect to the ground state of the system. In absence of interactions, the relevant information is encoded only in the 2-point correlation function G (2) , with higher-order correlations G (N >2) decomposing as a sum of products of only . This is the case when the ground state of the system is Gaussian. In order to characterize the influence of interactions, it is useful to decompose the N -th point dis refers to the disconnected part of the correlation function and it is determined by lower order correlation functions G (N <N ) . In this sense, the G (N ) dis does not possess any new information of the system at order N . On the other hand, G (N ) con , the connected part of the correlation function, has access to proper and characteristic information about the system at order N . As a result, a complete factorization of higher-order correlation functions as in the Gaussian case, is equivalent to have G (N ≥2) con = 0. In other words, the connected part of N -point correlation functions distinguishes the ground states of interacting theories from those of noninteracting ones: i.e., while for Gaussian states the connected correlation functions of order higher than two vanish, those of interacting systems are generally nonzero. This is the reason for which, despite being successful in capturing some nonperturbative effects through the gap equation, it would desirable, in order to study the phenomena of strongly interacting systems, to have powerful tools able to quantify the effect of perturbative and especially, nonperturbative quantum processes that contribute to the connected part of N th-point correlators.

Non-Gaussian cMERA
In this section we present a systematic and model independent 4 formalism to go beyond the Gaussian approach in interacting field theories. This formalism was generalized to nonperturbatively build non-Gaussian cMERA wavefunctionals in [20]. The method relies in performing a set of nonlinear transformations on the fields of the theory that extends the linear transformations defining the Gaussian approach for free theories.

Nonlinear Canonical Transformations
According to [21][22][23], extensive non-Gaussian trial states wavefunctionals can be built as Hadamard's lemma 5 , This reduces the calculation of expectation values of functionals to a finite number of Gaussian expectation values. The exponential nature ofŨ ensures the correct extensive volume dependence of observables and specifically the total energy of the system. Furthermore, asŨ is unitary, the normalization of the state is preserved. The operator B consists of a product of π's and φ's, which is given by with m ∈ N. We will denote these operators from here in advance symbolically as B ≡ π φ m . Here, s is a variational parameter that, as it will be shown later, tracks the deviation of any observable from the Gaussian case. h(p, q 1 , . . . , q m ) is a variational function that must be optimized upon energy minimization. It is symmetric w.r.t. exchange of q i 's and is constrained to satisfy: These conditions ensure that the commutator series terminates after the first nontrivial terms. Namely, the constraints in (51) are the responsible for this truncation when the Hadamard's lemma (22) is applied. The action ofŨ on the canonical field operators φ(p) and π(p) is given byφ where the quantities with a bar are defined as the nonlinear field functionals, Thus, we are considering a class of field transformations, The constraints (51) on the non-Gaussian variational parameter h(p, q 1 , . . . , q m ) can be accomplished by taking the decomposition where it is imposed that η(p) · ζ(p) = 0, i.e., the domains of momenta, where η and ζ are different from zero have to be disjoint, up to sets of measure zero. A suitable ansatz for η and ζ is given by [21,22] where ∆ 0 and ∆ i (i = 1, . . . , m), are variationally optimized, coupling dependent momentum cutoffs such that |∆ 0 ,i |≤ Λ. Γ(x) refers to the sharp cutoff function in (28).
the half-mean width of the functional is proportional to (k 2 + µ 2 ) −1/4 . Thus, as the variational mass µ increases in the Gaussian kernel F (k), a stronger suppression for nonclassical configurations than in the free case occurs. This is what happens when the Gaussian ansatz is used, for instance, in the λφ 4 theory [16]. That said, we now illustrate the action ofŨ on a Gaussian wavefunctional by choosing the transformation B = π φ 2 for clarity: Hence, while for the Gaussian case the decaying slope of the wavefunctional is now the decaying slope is corrected by new terms as Observables: the calculability of the ansatz In this sense, it is of particular interest to consider n-point correlation functions. In general, we will have where the subscript N G refers to an expectation value taken w.r.t. the non-Gaussian state (48). To evaluate this, we use (52) and (53), which yields That is to say, the calculability of the ansatz allows to compute the expectation value of ob-

Non-Gaussian cMERA Formalism
Following [20], now we use the method depicted above to generate cMERA non-Gaussian trial statesΨ[φ, u] implementing a renormalization group flow of the wavefunction for interacting field theories. This can be cast as where U G (u) = U G (u, u IR ) = P e −i u u IR (K(u )+L) du is the cMERA unitary operator for the free theory defined through the variational parameter (28) andŨ = exp (B). To go beyond Gaussian approach, we have to consider operators B that at least are cubic in the products of π and φ fields. Here, we are focusing on the simplest case B = πφ 2 . Recalling (30), the transformation U acts on the fields as follows: where, for compactness, we have dropped δ(p − q 1 − q 2 ) appearing inside the integrals and we have definedh From a cMERA point of view,h(p, q 1 , q 2 ; u) can be interpreted as a variational, scale and coupling dependent momentum cutoff function. While this happens automatically in the nonperturbative cMERA circuit, let us elaborate more on this now. In the Gaussian entangler the variational parameter g(p; u) was decomposed in (28) as i.e., a scale dependent function g(u) times a sharp momentum cutoff function Γ(p/Λ).
Similarly, we would expect the nonlinear transformations on the field modes (66) (and thus, any observable built from products of these field modes) to have a similar structure, h(p, q, r; u) ≡ g(p, q, r; u) · Γ B (p, q 1 , q 2 ) , where Γ B is a generalization of the Γ cutoff function in (28). However, from (67) and (57) we identify Then,h(p, q, r; u), which is a variational coupling dependent momentum cutoff function, implies the optimization of both the Gaussian parameter f (p; u) and the cutoff momenta Λ 0,i . As we will see in Section 5, this variational scheme captures nonperturbative and non-Gaussian interaction effects, which turn out to be essential at the regime at which the Gaussian quasi-particle picture is no longer valid. 6 With this, in [20], using the set of scale dependent nonlinear transformations given in (66), nonperturbative cMERA states for interacting field theories were built where, as before, it is straightforward to show that the non-Gaussian scale-dependent wavefunctional Ψ[φ; u] can be written asΨ Nevertheless, from a circuit/tensor network viewpoint, this methodology seems incomplete. Namely, the sequence of scale dependent non-Gaussian wavefunctionalsΨ[φ; u] is not generated by a cMERA circuit with an entangler containing the whole set of transformations, including the nonquadratic ones. This would be a satisfactory and more natural continuum generalization of the lattice MERA algorithms. However, this is trivially guaranteed, as the product of two unitaries turns out to be another unitary. In this sense, an important result of this work is to explicitly obtain such unitary operator using the firm ground provided by the conceptual framework presented in this section.

Fermionic non-Gaussian cMERA
Before going into the proposal of a cMERA circuit that implements the scale dependent nonlinear canonical transformations for an interacting scalar field, let us comment on a formulation of nonlinear canonical transformations for fermionic fields. As advanced in [20], the method to generate non-Gaussian cMERA wavefunctionals applies to fermionic field theories. Here we present some transformations.
Despite they can be generalized to any dimension and/or fermion, in this work, we will consider non-Gaussian transformations for two-dimensional Dirac fermions. Being model independent, our proposal is specially well suited to analyze the Gross-Neveu model (GN) [29]. This is a renormalizable, asymptotically free two-dimensional theory which displays chiral symmetry breaking and dynamical mass generation. The model describes where ψ andψ ≡ ψ † γ 0 are two-dimensional Dirac spinor fields, g 0 denotes the coupling constant and the γ-matrices are given by with σ i being the Pauli matrices. It is useful in order to deal with the GN model, to decompose the Dirac spinor ψ, into its chiral projections, ψ = ψ L + ψ R , In terms of the chiral projections, the GN model action can be written as As in the bosonic case, we build extensive non-Gaussian fermionic trial states as wavefunctionals of the form where Ψ G [ψ L , ψ R ] is a normalized Gaussian wavefunctional andŨ F = exp(F), with F, an Hermitian operator that nonperturbatively adds the new variational parameters s and g(q 1 , · · · , q 4 ) to those in the Gaussian wavefunctional. Here, we consider the transformation F = s Using the Fierz identities, we can prove that momenta q 2 and q 4 can be exchanged by including an extra minus [30,31]. From here in advance we will drop the momentum conservation delta. The Hermiticity of F imposes Upon this constraint, the effect of the transformation on the chiral spinor fields is given by: where, for the spinor index α, we have Due to the Hermiticity of the operator F, the canonical anticommutation relations are preserved for the transformed Ψ L,R (p) fields. With this, as in the bosonic case, it is possible to build a set of scale dependent non-Gaussian fermionic wavefunctionals using (78) and the Gaussian cMERA for free Dirac fermions given in [3,4].
Being formally guaranteed, we will discuss the scale dependent version of these non-Gaussian transformations in a future work.

icMERA: non-Gaussian cMERA Circuits
In this section we first introduce a nonquadratic cMERA entangler to generate non-

Nonquadratic Entangler
Based on the concept of scale dependent nonlinear canonical transformation and non-Gaussian cMERA wavefunctionals exposed above, in this section we propose, as a definition of a nonperturbatively built non-Gaussian cMERA circuit (icMERA), the Hamiltonian evolution produced by the entangler where K 0 (u) is the quadratic Gaussian entangler given in (27) and B(u) is This is a scale dependent operator that nonperturbatively incorporates nonquadratic interaction terms to the cMERA evolution through the variational function g(q 1 , q 2 , q 3 ; u), with s being a variational parameter, as in the previous section. g(q 1 , q 2 , q 3 ; u) is symmetric under exchange of momenta q 2 , q 3 and is parametrized in terms of variational cutoff functions.
Thus, as a generalization of the Gaussian cMERA, we define the icMERA evolution operator in the interaction picture as Having a real wavefunctional requires time reversal invariance of the icMERA evolution which in turn amounts to having an odd number of π operators in B(u). In this case B(u) is formally equivalent to B = π φ 2 , i.e., it incorporates cubic interactions into the cMERA evolution. Furthermore, as s in (83) is a variational parameter related to the coupling strength of the theory, the standard Gaussian cMERA evolution (15) is recovered when With this, one may write the icMERA wavefunctional in the Schrödinger picture as The action of the icMERA operator (84) on the field modes φ(p) is given by, Defining the following scale integrals as we obtain (88) with g(p; u) given by (28) and having imposed the constraints (51) on f (p, q 1 , q 2 ; u). That is to say, in contrast to the function h in (50), we impose constraints on variational functions once the integration over u has been done. For the function f (p, q 1 , q 2 ; u), which satisfies g(pe −u , q 1 e −u , q 2 e −u ; u) = ∂f (p, q 1 , q 2 ; u) ∂u , we will assume the following ansatz where g B (u) modulates the strength of the scale u and Γ B (p, q, r), which is a generalization of the Gaussian Γ(x) function cutoff, satisfies the conditions (51). While in this work, we fix g B (u) = 1, the case for a generic smooth function will be addressed in a future.
Following the discussion in Section ??, we pick with Once again, we denote Γ(x) ≡ Θ(1 − |x|) with Θ(x) being the Heaviside step function while ∆ 0 and ∆ 1 can be understood as two variational and coupling dependent momentum cutoffs, with ∆ 0 , ∆ 1 ≤ Λ. The optimal function Γ B (p, q, r), which will depend on the scale through the momenta, must be found self-consistently by determining both cutoffs.
Different from the Gaussian set-up, this scheme illustrates how the strength of the interaction variationally determines the region in momentum space that will be relevant in the optimization procedure. This connection turns out to be essential for strongly-coupled systems, which exhibit some regimes at which the Gaussian quasi-particle picture is no longer valid. In Section 5.3, we give evidence of some nonperturbative effects captured by this method.
With this, the action of the icMERA operator K(u) + L on φ(p) is given by where for convenience we have dropped the momentum conservation delta δ(q 1 + q 2 − p) in the integrand. In the same manner, the action on π(p) can be written as (95) e iK π(p) e −iK = e f (p,u) π(p) + 2s This yields a full action of the icMERA operator on π(p) given by (96) That is to say, changing the notation to f (p, u) ≡ f u (p) and f (p, q 1 , q 2 , u) ≡ f u (p, q 1 , q 2 ), the scale dependent nonlinear transformations on the fields generated by the icMERA circuit (84) arẽ π(e −u q 1 )φ(e −u q 2 ) . (97)

Generic non-Gaussian transformation
When the non-Gaussian entangler B(u) contains n scalar fields (in analogy to B = π φ n introduced in the previous section), (where we have dropped the momentum conservation δ(k+q 1 +· · ·+q n )), then the action of the icMERA operator is given bỹ whereh u andȟ u are defined as and At this point, it is worth to discuss (100) in more detail. We have built the icMERA circuit in terms of an entangler K(u) = K 0 (u) + B(u) that includes nonquadratic π φ n interaction terms in a nonperturbative way by means of B(u). The crosstalking between the quadratic part K 0 (u) and the purely non-Gaussian term B(u) is obvious from (100), where we appreciate that the coefficients of the nonlinear transformations on the fields induced by icMERA, are the product of the variational parameter f u (p, q 1 , · · · , q n ) related with B(u) multiplied by a term that is a function of the variational parameter f u (p) that defines the Gaussian term of the entangler. In this sense, the quadratic Gaussian part analysis, let us assume that in an optimized icMERA state, the variational parameter of the quadratic term, f u (p), has the simple, Gaussian form where µ Λ is a variational mass parameter that depends on the interaction strength.
As it will be shown when discussing the optimization of the circuit, this assumption is correct for a variety of situations of physical interest. With this, we fix our attention into the term in (100). Namely, the quadratic Gaussian contribution of the icMERA entangler K(u) suppresses the contribution of the non-Gaussian term f u (p, q 1 , · · · , q n ) as one probes the system with very small momenta, i.e., close to the IR. This suppression smoothly grows and takes a value O(1) as one reaches the UV, i.e., p ∼ Λ. This is the regime at which the non-Gaussian contribution to building the state is maximum. In this sense, we interpret the quadratic term K 0 of the icMERA ansatz as a regulator which drives the IR state to be purely Gaussian.
From the point of view of the entanglement flow along an icMERA circuit, as posed in [16], the quadratic term is responsible for generating pairwise entanglement between modes as a function of scale [9], whereas the B(u) term generates n-tuplet quantum correlations as a function of scale. How to characterize these higher order quantum correlations is an interesting topic that we leave for future investigations.

icMERA Correlation Functions
In other words, the icMERA circuit goes beyond the Gaussian approximation and captures scale dependent nonperturbative contributions for the N -th order correlator, which are arranged in powers of the variational parameter s. As commented above, in order to quantify to which extent the icMERA ansatz is nonperturbatively characterizing the interactions of the theory under consideration, we have to explicitly calculate the expressions for the connected part of these correlators. In this respect, the connected 2-point and 4-point correlation functions in real space (see also [20]) given by the icMERA circuit with a cubic non-Gaussian entangler B(u), i.e., with n = 2 in (105), are given by G (2) c (ab; u) = D(ab, u) + s 2 χ 2 (ab, u) , where D(ab; u) is and bracketed quantities, which are given in Appendix A, correspond to a series of per-mutations of the loop integrals 7 χ 2 (ab, u) = 1 2 p,q Υ 2 (p, q, u)F (p u ) F (q u ) e i(p+q)·x ab χ 5 (ab, cd, ef , u) = p,q,r Υ 5 (p, q, r, u) F (p u ) F (q u ) F (r u ) e i(p·x ab +q·x cd +r·x ef ) χ 6 (ab, cd, ef , gh, u) = p,q,r,s × e i(p·x ab +q·x cd +r·x ef +s·x gh ) . (108) We have introduced the variational "vertices" Υ 5 (p, q, r, u) = c(p, q; u) c(p, r; u) , Υ 6 (p, q, r, s, u) = c(p, q; u) c(p, r; u) c(q, s; u) c(r, s; u).
and have used the compact notation a ≡ x a , ab ≡ x ab ≡ x a − x b with p u ≡ p e −u . The scale dependent variational functions are encoded in c(p, q; u), which is given by It is worth to mention that icMERA, which involves scale dependent nonperturbatively generated wavefunctionals, allows us to study regimes that interpolate between weak and strong coupling. For example, the expressions for the connected parts of the 2-and 4-point functions in (106) are fully determined and, depending on the optimal value of s, they can receive both perturbative and nonperturbative contributions. As it will be shown in the next section, in the self-interacting scalar λ φ 4 model, the parameter s ∝ λ. As a result, one may infer from (106) that in the perturbative regime, G c ∼ O(λ 2 ) while the nonperturbative effects are captured by the term s 4 .

icMERA Circuit for the Scalar λ φ 4 Theory
In order to fully solve the icMERA tensor network and evaluate the previous expressions for a concrete theory described by a Hamiltonian H, one must obtain the optimal values for the variational parameters f (p, u), f u (p, q 1 , · · · , q n ) and s. This is addressed by minimizing the expectation value of the energy density w.r.t the icMERA ansatz for a fixed length scale u, i.e., H u = Ψ u | H |Ψ u . Our aim here is to obtain the optimized parameters for an icMERA tensor network circuit representing the ground state of the self-interacting λφ 4 scalar theory in (1+1) dimensions. The Hamiltonian density for this model reads where m and λ are the bare mass and the bare coupling of the theory respectively. The λφ 4 scalar field theory in two dimensions provides an example of a nontrivial interacting field theory. According to [34], this model experiences a second order phase transition at which the vacuum changes continuously from a symmetric to a nonsymmetric state. However, the rigorous proof [35] of this fact does not allow to compute the critical coupling. An estimate was obtained by the variational Gaussian approximation [36], but this yields a wrong critical behavior as it predicts a first order phase transition.
The momentum integrals χ's, which are given in AppendixA, here are evaluated at the same spatial point x and The variational terms in the integrals χ i can be understood as a kind of generalized condensates yielded by the icMERA ansatz. To see this, we note that, for small s, the term ∼ φ c χ 3 (u) in (113) is the major contribution to the improvement of the energy value compared to the Gaussian estimate (which is obtained when taking the limit s → 0 in (113)) [20][21][22]. This term is formally equivalent to the one coming from the interaction with a background field. In a Gaussian ansatz, that background field is given by χ 0 in (31).
Indeed, it has been shown that the optimal χ 3 (u) contains an infinite series of contributions to the two-point function that correspond to the "cactus"-diagrams resummation obtained from a pure Gaussian ansatz [24]. This is in agreement with understanding χ's in (113) as a series of generalized non-Gaussian condensates.
In this sense, in some limits the Gaussian wavefunctional (31) turns out to be a special case of the πφ 2 -kind of wavefunctionals yielded by icMERA [22]. The idea is to define P and Q as the domain supports for the transformed p-modes of the field and the support for the shifting q-modes respectively (see (97)). P is a sphere with volume V P and center at the origin and Q is a spherical shell which surrounds P with volume V Q . Taking the limit of small ε ≡ V P /V Q it is possible to show that the only "condensate" independent of ε is χ 1 while for the remaining χ's one can easily obtain upper bounds which depend on ε → 0 [22]. By considering φ c as a new parameter, the final energy expectation value may be obtained directly from the Gaussian result by substituting χ 0 → φ c . Again, this strongly suggests that the parameters χ act as a kind of "higher order" non-Gaussian condensates that expand the ability of the ansatz to improve the variational estimation of ground state energy.

Optimization
The effective potential computed with the icMERA state Ψ[φ; u] is defined as, The optimal values of the variational parameters f (p; u), s and f (p, q, r; u) have to be found by setting Despite this can be done in full generality, φ c has to be fixed, in order for the trial wavefunctions to be consistent with the Rayleigh-Ritz method [32,33]. In other words, one has to fix φ c , and minimize with respect to the rest of the variational parameters. In this form, this yields a set of nonlinear coupled equations that must be solved numerically and self-consistently.
From the above equations it is clear that V G (φ c ; u) only depends on the variational function f (p; u). It is thus reasonable to proceed as follows: we first optimize the variational parameter f (p; u) by choosing an appropriate scale, in our case the UV, by taking This yields exactly the same result as Eq. (43), i.e., F (p; 0) ≡ F (p) reduces to F (p) = 1 with µ a mass variational parameter given by the gap equation, With this, one may fix the variational parameter f (p; u) of the icMERA circuit through (32), that is to say, f (p; u) reduces to the simple Gaussian form given in Section 5: with ω Λ = Λ 2 + µ 2 .
Once f (p; u) has been fixed through µ, one is left to optimize the variational function f (p, q, r; u) defined in (92) and (93). According to those equations, determining that function amounts to finding optimal values for the variational momentum cutoff parameters ∆ 0 and ∆ 1 (renamed ∆ M and ∆ N respectively from here in advance). To this end, for the fixed value of µ obtained through the gap equation, we consider the non-Gaussian corrections to the effective potential ∆V (φ c , µ; u) as a function of s, ∆ M and ∆ N . Then, by imposing δ∆V (φ c , µ, ; u) δf (p, q, r; u) = 0 and δ∆V (φ c , µ, ; u) we must find the optimal values of s, ∆ M and ∆ N . The solution for the optimal s is given by s = − 2η 2 χ 3 (u) λ φ c 2η 2 (2m 2 χ 2 (u) + λχ 5 (u) + 4χ 7 (u)) + λχ 2 (u) log (m 2 /µ 2 ) , with η ≡ 2π. When plugging the expression for the optimal s into ∆V (φ c , µ; u), we obtain an expression that only depends on the χ-integrals. Thus, finding the optimal values of the variational cutoffs ∆ M and ∆ N , immediately leads to establishing the optimal value of s.
To this end, and in order to have a fully operational equation to optimize these parameters, we note that the χ-integrals defining ∆V (φ c , µ, ; u), assuming that µ can be explicitly written in terms of ∆ M and ∆ N (see expressions given in Appendix A).
This allows us to express ∆V (φ c , µ; u) only in terms of the variational cutoffs and thus carry out a direct numerical minimization over them.
To finish the procedure it is necessary to choose a convenient scale u to carry out the numerical optimization. A significant further simplification of the equations occurs by choosing the scale u IR → −∞. With this, the leading behavior of the optimal s is given bys where the barred quantities refer to their values when evaluated at ∆ M , ∆ N and u IR .
With the optimal values of the parameters µ, ∆ M and ∆ N at hand, the icMERA circuit could be understood as a device for probing, at different length scales, the labyrinth of variational mass scales created by these parameters.
The optimized icMERA circuit allows now, for any given λ and φ c , to calculate expectation values of observables at different length scales. To illustrate this, in Figure 1, we show the behavior of ∆V (φ c ) for two different values of λ and various renormalization scales. We remark that the results which are consistent with our approximations are those close to φ c ∼ 0. With this, we note that icMERA detects the second order phase transition commented above. While a simple Gaussian approximation yields a first order phase transition, remarkably, this phase transition is turned to second order by the influence of the πφ 2 kind of transformation of the optimized icMERA state [21,22]. The critical behavior, inherent in the effective potential, can be analyzed through the behavior of ∆V (φ c ≈ 0; u). Indeed, we obtain that for λ crit ≈ 41.4. For coupling constants bigger that λ crit , e.g., for λ = 61.2, the derivative above is strictly negative, thus signing that a second order phase transition has occurred.
Noteworthy, a pure Gaussian ansatz would have predicted a first order transition at this point. These results agree with those posed in [21,22] and thus are taken as a validity test for our optimization procedure.

Correlation Functions
As mentioned above, once the optimal variational parameters of the icMERA-πφ 2 circuit, f (p, u) and f (p, q, r; u), are obtained for the λ φ 4 theory, then higher-order correlation functions can be computed through equations (106). As commented previously, the knowledge of higher order correlation functions is necessary in order to distinguish the ground states of an interacting from of a noninteracting system. With the aim to illustrate the performance of an optimized icMERA circuit, we have carried out computations of the connected part of two and four point correlation functions in (106) for the λ φ 4 theory.
To this end we have optimized the icMERA circuit for this model under the prescriptions    approaches the deep IR regime. Our results also show that icMERA is able to capture proper and characteristic scale dependent information about the system at larger orders than those provided by the Gaussian ansatz. Such scale dependence occurs because of the variationally optimized "vertices" Υ 2 , Υ 5 and Υ 6 in Eq. (109).

Discussion
The non-Gaussian circuit icMERA presented in this work, introduces a new variational tool to address strongly interacting field theories by means of a systematic building of non-Gaussian wavefunctionals. In this sense, icMERA provides a tool to study how the structure of interactions are encoded in the correlations and the entanglement patterns of the wavefunctionals of interacting field theories. In our proposal, we have shown that an icMERA circuit can be viewed as a device for probing at different length scales, a variational labyrinth of momentum scales that give account for the structure of the interactions in a theory. Regarding the aforementioned entanglement patterns, those are well understood for the case of free theories in terms of the RG flow implemented by the Gaussian cMERA. In the case of interacting theories, it is expected that non-Gaussian correlations establish more complex patterns of entanglement at different length scales.
Thus, it would be very interesting to carry out a systematic study of these quantum correlations in future works. The multiscale approach provided by an icMERA circuit, may be useful to address recent experimental data on higher order correlation functions in many body systems [37,38]. As a fact, despite the most fundamental laws of physics are usually explored in experiments probing the smallest distances, it has been recently shown that models designed to give account of the observations in these high energy experiments, may emerge as effective descriptions of many-body systems at lower energies, e.g. in condensed matter physics and quantum simulation experiments. It is worth to investigate how an icMERA circuit could give account of these data by fixing the laboratory cutoff and the scale energy at which the experiment is performed.
Finally, having such a robust prescription to address the entanglement renormalization of interacting theories at a nonperturbative level, we expect to unveil in a near future the concrete holographic realizations that icMERA is able to exhibit. For example, it would be interesting to probe the definition of complexity in QFT from a cMERA point of view [8,27,39,40]

A χ integrals
The loop integrals χ i which are related to the circuit π φ 2 depend on both positions and the renormalization scale u. Once the optimal variational parameters f (p; u) and f (p, q 1 , q 2 ; u) are obtained for a particular theory, then higher order correlation functions can be computed through them. Denoting c(p, q; u) ≡h u (|p + q|e −u , pe −u , qe −u ) , their explicit expressions can be written as χ 2 (ab; u) = 1 2 pq c(p, q; u) 2 F (pe −u )F (qe −u ) e i(p+q)·x ab , × F (pe −u )F (qe −u )F (re −u )F (se −u )e i(p·x ab +q·x cd +r·x ef +s·x gh ) , + c(|p + q|, −q; u) 2 F (qe −u ) F (pe −u ) −1 .
With this, the quantities in brackets appearing in (106) are given by where the explicit dependence on u has been dropped for clarity.

A.1 Asymptotic form of χ integrals
Once we pick the ansatz for the variational function f (p, q, r; u) in the form given by where η ≡ (2π).