{'fraction_non_alphanumeric': 0.03378538334066083, 'fraction_numerical': 0.009537015117183537, 'mean_word_length': 4.608877086494689, 'pattern_counts': {'":': 0, '<': 3, '<?xml version=': 0, '>': 5, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 7, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Non-Gaussian Variational Wavefunctions for Interacting Bosons on the Lattice

A variational method for studying the ground state of strongly interacting quantum many-body bosonic systems is presented. Our approach constructs a class of extensive variational non-Gaussian wavefunctions which extend Gaussian states by means of nonlinear canonical transformations (NLCT) on the fields of the theory under consideration. We illustrate this method with the one dimensional Bose-Hubbard model for which the proposal presented here, provides a family of approximate ground states at arbitrarily large values of the interaction strength. We find that, for different values of the interaction, the non-Gaussian NLCT-trial states sensibly improve the ground state energy estimation when the system is in the Mott phase.


I. INTRODUCTION
One of the major problems in many body physics is to understand the phenomena associated to strongly coupled systems. This includes a huge variety of effects ranging from quark confinement to strongly correlated electron systems in condensed matter physics. To do so, nonperturbative methods are required. Many useful techniques such as path integral approaches, large N expansions or numerical methods (such as density matrix renormalization group or tensor networks) have been developed to address these problems. While these methods have successfully been applied to a broad range of problems, the complexity of the theoretical proposals have obscured the understanding of the physical phenomena. In this sense, it is acknowledged that the use of variational methods allow to tackle these problems to some extent in a relatively simple way by means of variational wavefunctions. When it comes to nonperturbative aspects, there are situations in which the use of wavefunctionals exhibits clear advantages. For example, path integral methods are especially suited to compute quantities that have no perturbative contributions and can be addressed through a saddle point approximation. Nevertheless, in case the observables of interest can receive both perturbative and nonperturbative contributions, the path integral approach becomes more difficult [1].
Choosing appropriate variational wavefunctions for a strongly interacting many body system is a challenging issue [2]. 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. Due to the enormous size of the Hilbert space in a many body system, it is very difficult to identify by mere intuition the relevant features that have to be grasped by the ansatz. Thus, a systematic method to build them would be desirable. Secondly, one must face the problem of calculability. Even possessing a reasonable and flexible ansatz for the wavefunctional, one wants to evaluate expectation values of operators/observables of interest in this state which, in general, will be challenging. Given the very limited ability to evaluate expectation values with non-Gaussian wavefunctionals, the calculability requirement on the trial wave functional has constrained the form of the trial wavefunctionals to Gaus-sian states.
Gaussian states are given by the exponentials of quadratic functionals of creation and annihilation operators of the fundamental fields of the theory under consideration. The expectation values of physical observables can be efficiently computed for these states [3], which obey Wick's theorem, thus allowing to express expectation values of arbitrary products of mode operators in terms of products of pairs [4]. Gaussian approximations as the Hartree-Fock-Bogoliubov one are used to approximate the dynamics of interacting bosons [5,6]. While Gaussian states represent the exact ground state in noninteracting systems, they have some important limitations in capturing phenomena related to interacting systems. Therefore, it is interesting to extend those ansatze to non-Gaussian generalizations. Given, however, the very limited ability to evaluate expectation values with non-Gaussian wavefunctionals, the calculability requirement on the trial wave functional is typically restricted to the form of the trial wavefunctionals to Gaussian states.
Finally, one must address the problem of the ultraviolet modes. The main objective of a variational calculation in a strongly interacting system 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 system 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.
In this work, following the pioneering works [7][8][9] we develop a class of variational non-Gaussian wavefunctions which extend Gaussian states by means of nonlinear canonical transformations (NLCT) on the fields of the theory under consideration. The presented scheme is self-consistent and requires no other assumptions than the choice of a NLCTvariational manifold. These variational wavefunctions, has been shown to possess strong entanglement between the microscopic degrees of freedom [10], while retaining most of the calculability of Gaussian wavefunctions. There exist several related methods, like the cluster expansion [11,12], the t-expansion [13], and the δ -expansion [14], which, even going (slightly) beyond the Gaussian ansatz, they yield expectation values that cannot be exactly calculated and have to be arXiv:2211.04320v1 [cond-mat.str-el] 8 Nov 2022 approximated by additional series expansions.
We illustrate the NLCT method with the Bose-Hubbard model in one dimension. Our approach provides a family of approximate ground states at arbitrarily large values of the interaction strength. Quantum interacting bosons in 1D provide an exciting area where condensed matter, low-temperature, and ultracold atomic converge. The Bose-Hubbard (BH) model describes a system of interacting spinless bosons on the lattice. This model provides a theoretical description of interacting cold atoms in optical lattices [15,16]. Besides, some numerical solutions exist for the BH model, that allows us to benchmark our non-Gaussian ansatze [15].
In line with our approach, in Ref. [17], authors built non-Gaussian wavefunctions for many body systems in the lattice using a set of unitary transformations on Gaussian states. The nature of the transformation depends on the theory under consideration. Its general form was inspired by canonical transformations in condensed matter physics, such as the polaron transformations in electron-phonon systems and flux attachment in FQHE. While the construction in [17] may work for mixed boson-fermion systems, it offers a very limited class on purely non-Gaussian bosonic states. This is due to the challenge posed by what we call the truncation problem for bosons that will be discussed later in this paper.
In the context of high energy physics, NLCT-wavefunctions have been recently used to build a nonperturbative version of continuous entanglement renormalization tensor networks in order to explore connections between tensor networks and the AdS/CFT holographic correspondence [18][19][20]. The tensor network circuit there, implements a series of scale-dependent NLCT. It was shown that the leading contribution to the entanglement entropy, comes from the Gaussian part of the ansatz and are always related to the leading area term in the holographic calculation. On the other hand, the subleading contributions are given by the non-Gaussian part of the ansatz and are related with quantum corrections to the holographic entanglement entropy.
The paper is structured as follows: A description of the NLCT transformation method and some technical aspects, as the calculation of expectation values, are presented in Section II. In Section III we introduce the 1D Bose-Hubbard model and briefly review the variational approach to its ground state using Gaussian states (Gaussian Variational Approach, GVA). In Section IV we detail the calculation of the energy expectation value for the BH Hamiltonian in a concrete NLCT-Non-Gaussian state. We discuss on the suitability of this choice of NLCT-state for the problem at hand and provide expressions for the energy functional that will be subsequently optimized. Section V details the optimization of the ground state energy functional and the results where, our approximation to the ground state of the 1D Bose-Hubbard model for arbitrarily large values of the interaction strength, shows a sensible improvement on the estimation of the ground state energy with respect to the Gaussian case, especially when the system is in the Mott phase. We finish with a discussion on the results and an outlook in Section VI.

II. NON GAUSSIAN VARIATIONAL ANSATZ
In a bosonic QFT on a lattice, let us consider a normalized Gaussian variational wavefunctional where the creation and annihilation operators for the bosonic modes are denoted by b † k and b k respectively such that [b k , b † k ] = δ k,k , and Λ kl is a matrix depending on the variational parameters λ . Then, extensive non-Gaussian trial wavefunctions can be non-perturbatively built as with U(h) = exp(B(h)), and B(h) † = −B(h) an anti-Hermitian non-quadratic operator (in terms of b † k and b k ) that non-perturbatively adds new variational parameters labelled by h to those λ defining |Ψ(λ ) G .
The exponential nature of U ensures the correct extensive volume dependence of observables and specifically the total energy of the system. As U is unitary, the normalization of the state is not affected. It is straightforward to see that the expectation value of any operator O with |Ψ(λ , h) NG ≡ |Ψ NG amounts to the calculation of a Gaussian expectation value for the transformed operatorÕ = U † O U, In principle, for bosonic systems, both in finite and infinite lattices, any non-quadratic choice for B(h), while leading to a non-Gaussian trial state, it induces an infinite commutator expansion via the Hadamard's lemmã This spoils any possibility on having an ansatz with finite calculability properties as any computation of an expectation value, amounts to the evaluation of an infinite series of Gaussian expectation values. This is what we call the truncation problem. In fermionic systems in finite lattices, this basic problem is alleviated by the anticonmuting nature of its operator algebra. Namely, for a non-quadratic fermionic operator has only a finite number of terms in its formal power expansion series or, equivalently, [F, O] n = 0 at least for n ≥ N, with N the number of sites in the lattice. This of course does not ensure that any further structure is needed in order to find fermionic ansatze with nice calculability properties (see [19] for examples in fermionic field theory), but it is worth to mention that the truncation problem is especially ill-posed in bosonic systems.
Remarkably, the method presented here, provides an ansatz for bosonic systems that automatically implements a controllable truncation in Eq. (5). This reduces the calculation of expectation values of operators to a finite number of Gaussian expectation values. The operator B consists of a product of bosonic operators π's and φ 's, which is given by We will denote these operators from here in advance symbolically as B ≡ π φ m . Here, s is a variational parameter that 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 (5) terminates after the first nontrivial term. Namely, the constraints in (9) are the responsible for this truncation when the Hadamard's lemma is applied. The action of U on the canonical field operators φ k and π k is given bỹ where Φ k , Π k are defined as the nonlinear field functions, Being U unitary, the canonical commutation relations (CCR) still hold under the nonlinear transformation of the fields, giving [φ p ,π q ] = iδ p,−q . For this reason, the above transformations are known as nonlinear canonical transformations (NLCT).
Regarding observables, it is of particular interest to consider n-point correlation functions φ k 1 · · · φ k n NG . To evaluate this, we use (9) and (10), to obtain where · · · refers to a an expectation value taken w.r.t. the Gaussian state. That is to say, the calculability of the ansatz allows us to compute the expectation value of observables such as correlation functions in terms of a finite number of Gaussian expectation values. In particular, the terms proportional to s j in the non-Gaussian n-point correlation function correspond to (n + m( j − 1))-point Gaussian correlators, where j = 0, . . . , n and m is the power associated to the operator B = πφ m .

III. GAUSSIAN VARIATIONAL APPROACH TO THE BH MODEL
The Bose-Hubbard model describes a system of interacting spinless bosons on the lattice. It is a theoretical description in a wide variety of contexts such as interacting ultracold atoms in optical lattices, 4 He in various confined geometries and granular superconductors [21]. Its Hamiltonian in 1D is given by: In the first term, the kinetic energy, t is the hopping amplitude. The second term stands for the on-site interaction with repulsion strength U (> 0). For convenience, we include the chemical potential, µ.
The dynamics of the BH model is given by the interplay between boson tunneling (with amplitude t), and the repulsion between two bosons on the same site (of energy U). At fixed density, for small U/t, the bosons can be considered as nearly free, so at low temperatures T they condense into a superfluid with macroscopic occupation of the zero momentum singleparticle state. In the opposite limit of large U/t, the repulsion between the bosons localizes them into a Mott insulator. The Mott insulator is a state adiabatically connected to the product state with one boson in each potential minimum of the lattice. Interestingly, there is a second order quantum phase transition between these states which occurs at a critical value of U/t. This transition is driven by phase fluctuations and belongs to the XY universality class. In addition, there is another universality class for the Mott transition. A transition that occurs by changing the chemical potential (density), which is driven by density fluctuations and belongs to the mean-field universality class [15,22].
Here we discuss the Gaussian varational ansatze (GVA) for the BH model in D = 1. For that, it is convenient to write the Hamiltonian in momentum space: refers to the (non-interacting) dispersion relation with N the number of lattice sites. The generalization to D > 1 is trivial.
The GVA is based on building a Gaussian variational trial state given by with |Ω being the trivial vacuum b k |Ω = 0 and Here S a squeezed operator and D a displacement operator (which accounts for boson condensation in the broken symmetry phase) given by We note that S(λ ) is particular case of Eq. (1) where Λ k l is diagonal in the mode basis with Λ k l ≡ λ k δ k+l, 0 . The unitary transformation in (15) yields a canonical linear transformation (Bogoliubov transformation) on the field operators given by where u k ≡ cosh(λ k ) and v k ≡ sinh(λ k ) The GVA amounts to finding the optimal variational parameters that minimize the ground state energy, that is with upon which a numerical minimization must be carried out. It is, however, instructive to consider some previous analysis. By setting λ k = 0, we obtain a coherent state ansantz (S(λ ) = 1). In this case, the energy is minimized for |β 0 | equal to leading to the expectation value While the Gaussian ansatz (15) nicely captures some features of the BH ground state in the superfluid phase. the (15) wavefunction fails to capture the essential features of the Mott insulating phase. At zero order in t/U in the Gutzwiller approximation [23], the Mott-phase wavefunction is given by which it is obviously a non-Gaussian state (n 0 is the number of bosons on each site).

IV. NLCT-WAVEFUNCTIONS FOR THE THE BOSE-HUBBARD MODEL
Given the discussion above, here we apply the NLCT denoted by B = π φ 2 to the BH model. Before entering into many details, we discuss on the suitability of this choice for the problem at hand.
The low energy physics in the vicinity of the critical point between the superfluid and Mott insulator driven by the ratio t/U is described by a quantum field theory with an emergent Lorentz invariant structure [24]. In this low energy effective theory, the speed of 'light', c, is given by the speed of sound in the superfluid phase. In terms of the long-wavelength boson annihilation operator ψ, the Euclidean time action for the field theory is given by the quartic self-interacting scalar theory S = dτdx c 2 |∂ x ψ| 2 +|∂ τ ψ| 2 +w|ψ| 2 +u|ψ| 4 , where u is the coupling which tunes the system across the quantum phase transition at some u = u c . For u > u c , which corresponds to the Mott insulator, the field theory has a mass gap and no symmetry is broken. The gapped particle and antiparticle states associated with the field operator ψ correspond to the 'particle' and 'hole' excitations of the Mott insulator. These can be used as a starting point for a quasiparticle theory of the dynamics of the Mott insulator. The other phase with u < u c corresponds to the superfluid where the global U(1) symmetry of S is broken and a quasiclassical Gaussian theory of the superfluid phase is possible. This phase transition in the ψ 4 theory has been investigated through NLCT. More concretely, the B = π φ 2 has first been used in [7][8][9],and lately in the context of continuous tensor networks [18,19] and holography [20]. It is thus sensible to apply this transformation to our problem, by remarking that the choice for different transformations must include a justification for the regimes of the theory one is interested to study. Furthermore, let us elaborate on the effect of the NLCT transformation on wavefunctionals. In [8,19] it is shown that U generates a translation of the argument in the configuration space of the theory that symbolically reads with Φ = h φ 2 . In compact notation, writing Ψ G [φ ] = exp − 1 4 (φ · G −1 · φ ) , with G −1 the Gaussian kernel defining the correlation matrix of Gaussian states, the above identity can be cast as from which one immediately infers that the new wavefunctional has been enhanced with variational skewness and kurtosis terms (the term proportional to s and s 2 respectively) that cannot be captured by a Gaussian ansatz. For esxample, in [25], these terms have been characterized under different regimes of the Bose-Hubbard model. The bosonic field transformation (9) is explicitly given by where a summation index convention is assumed. In order to proceed it is convenient to provide an ansatz for the variational parameters h p,q 1 ,q 2 that fulfills the truncation constraints in (9). This can be easily achieved 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 where ∆ 1 and ∆ 2 , are variationally optimized, coupling dependent momentum cutoffs and Γ(x) = Θ(1 − |x|), with Θ the Heaviside step function. With this, our objective is to compute To facilitate the task we describe, below, the calculation of the energy with the non-Gaussian ansatz on a term-by-term basis.
A. Kinetic Term.
For the kinetic term we obtain where E kin G = ε 0 β 2 0 + ∑ k ε k v 2 k and After a lengthy albeit straightforward calculation the result for the kinetic term can be written as where, defining G k = 1 2 exp (−2λ k ), χ i are given by For the interaction term, we obtain As it will be justified below, our calculations will be carried out in such a way that s 2 s 4 notwithstanding that a large value of the interaction strength U can be taken. Due to this, the E int NG will take the form where The explicit expressions for the Σ's in terms of the variational parameters are rather lengthy so we refer the reader to Appendix to find them.

C. Final expression
At this point, and following [8] it is convenient to write the non-Gaussian energy density in terms of a different set of variables. To this end, we note that after the transformation π φ 2 , b k NG = β 0 + s χ 1 depends on several parameters of the ansatz. We define a new parameter ψ c ≡ sχ 1 , and, thus, one of the original variables can be eliminated. The resulting energy density is in general different from the Gaussian case and can be written as From this last expression it is straightforward to find an optimal s for a fixed value of ψ c and optimized values of ∆ 1 and ∆ 2 by, Analogously, it is possible to find a set of equations for finding the optimal values of the variational parameters ∆ 1 and ∆ 2 by In the case of the λ φ 4 theory these are a kind of Fredholm integral equations [9] that can be solved numerically. In our case the equations are rather involved so we use an alternative numerical procedure to find the optimal values. This will be described in the next section.

V. RESULTS
We carry out the optimization of the ground state energy functional in (40) following the lines exposed in [7,8,19]. This consists in obtaining the variational parameters of the ansatz in two consecutive steps. First, we optimize the Gaussian ansatz to obtain the optimalβ 0 ,λ parameters. Once these parameters are obtained, they are fixed in order to carry out the optimization of the non-Gaussian parameters in a separate fashion, that is, To further simplify the process, we fix ∆ 2 to its maximum allowed value of 2π(N − 1)/N leaving the parameter ∆ 1 (which can take values from 0 to 2π(N − 1)/N) as the only one parameter needed to determine the ground state energy.
After carrying out the optimization procedure described above, we note that the χ's and Σ's in (40), directly depend on the total number M of non-linearly modified modes φ k . Intuitively, this means that h k,p,q weighs how much a low energy mode φ k (with k ≤ ∆ 1 ) is modified by high energy modes φ p and φ q (with ∆ 1 ≤ p, q ≤ ∆ 2 ). According to our ansatz, a given value, for instance, the optimal value of the variational parameter ∆ 1 , defines the total amount of the non-linearly modified modes as M = (∆ 1 /2π) N < N.
Following [26], such dependence can be written as where a 2 , c 0 , c 1 and c 2 are slowly varying functions of t,U and µ. Given the definition of the independent parameter ψ c , that is, a non Gaussian correction to the Gaussian condensate, in terms of χ 1 , it is sensible to assume that ψ c = b M 1/2 . That is to say, assuming that b is an O(1) constant, the non-Gaussian correction to the condensate value is proportional to M 1/2 in such a way that its correction to the energy density is ∝ ψ 2 c ∼ M. Consequently, we obtain where ε G is the energy density of the Gaussian ansatz. In order to get further insights, we use that, after optimizing the energy functional (40), a numerical analysis yields α 2 = γ 1 = 0, γ 0 = γ 2 = 1, a 2 , c 1 , c 2 < 0 and c 0 > 0 which leaves Thus, defining z = ∆ 1 /2π, in the limit where z 1. i.e., for M N ,the optimal s can be written as This result implies that our truncated estimate of the ground state energy at order s 2 is justified and can be written as More explicitly, for µ > 0, the leading contributions to the variation of the energy density estimation w.r.t. the Gaussian case reads The last expression allows us to make some specific quantitative predictions for the behaviour of the ansatz at different points of the phase diagram. For instance, for a fixed t/U 1, one might expect that With Eq. (40), we compute the predictions of our ansatz for the ground state energy density. The performance of our method is tested by comparing this quantity with the ones obtained through a Coherent state, a Gaussian state and DMRG [27] (see Fig. 1). Both the energy densities obtained through a coherent and a Gaussian state are higher than the DMRG one, especially as one deeps into the MI phase. While, the Gaussian variational family provides a consistent class to approximate the ground state of the Bose-Hubbard model in the superfluid phase, our non-Gaussian approach is capable of both, i) providing a consistent non-perturbative (U-independent) truncation of the bosonic opèrator growth under a non-Gaussian transformation and ii) predicting a much better approximation to the ground state energy than previous methods. This is quite manifest in the MI phase, where the energy estimation improvement over the Gaussian ansatz is qualitatively explained by (52). Remarkably, the non-Gaussian ansatz performance on the superfluid phase equals the one obtained by the Gaussian ansatz, as far as the fraction of non-linearly modified modes z ∼ 0 in this phase.
We also show in Fig. 2 the energy density computed through the NLCT-ansatz, for varying t/U with µ fixed. The numerical results show agreement with the qualitative predictions yielded by (51) and with results in [25] where authors employed a high-order perturbative expansion to characterize the ground state of the Mott phase of the 1D Bose-Hubbard model.
Results in Fig (1) and Fig (2) are shown for a total lattice sites N = 32. We have numerically checked that the results are similar for larger values of N (N = 128, N = 512). This might be expected in terms of Eq. (51). There, the improvement over the Gaussian prediction is given in terms of the ratio Here, we have used that b ∼ O(1) parameter. Thus, numerical results based on a construction of the wavefunction ansatz that is fully non-perturbative, strongly suggests that our method accurately approximates ground state energy per lattice site across the critical point.

VI. DISCUSSION AND OUTLOOK
In this work, we have used a class of variational non-Gaussian wavefunctions which extend Gaussian states by means of NLCT. This technique was proposed as a variational non-perturbative method to study phase transitions in QFT [7][8][9]26] and has been recently developed for building controlled settings to study the AdS/CFT holographic duality through tensor networks [10,[18][19][20]. Here, the NLCT-wavefunction method has been applied to the 1D-Bose-Hubbard model. We obtain a family of approximate ground states for arbitrarily large values of the interaction strength. Our results show that a particular class of these states is able to sensibly improve the ground state energy estimation when the system is in the Mott phase [27].
We have used the ground state energy as a test-bed to benchmark the validity of the NLCT approach. It has been left as a future problem to address, to what extent the full nonperturbative structure of the ground state of the BH model is being captured by our concrete ansatz. This connects with the question of choosing the most appropiate non-Gaussian operator B = πφ n . The study of connected correlation functions, which can be systematically addressed through Eq. (11) (see [19,20]) might shed light on this. This amounts to quantify the skewness and kurtosis [25] of the system in different regimes and compare them with our ansatz. An additional benchmark to test the NLCT-ansatz is the prediction of the critical point (t/U)| crit for the superfluid-Mott insulator transition. For this, it is necessary to perform the explicit evaluation of the s 4 -term in the ground state energy functional, Eq. (40). This procedure has been used in the λ φ 4 theory [8] (see reference [28] for a comparison with other methods including tensor networks). For the BH model in particular, a comparison with [29] and [30]) might further help to elucidate to which extent our ansatz covers the non-perturbative structure of the BH ground state.
Another interesting possibility offered by our method is given by the translational effect on the wavefunctional argument exhibited in (25) and which has been partially explored in [10,20,31]. Given the calculability bonanza shown by the variational ansatze used in this paper, it is possible to com-pute the entanglement entropy of arbitrary regions in the BH ground state by using the prescription for the ground state of free theories. In free theories, the entanglement entropy is fully determined by the two-point correlation functions. For the interacting case, it is shown that non-Gaussian contributions to the entanglement entropy can be obtained through a Gaussian prescription by replacing the Gaussian two point functions by their non-Gaussian counterparts obtained by the NLCT method [10].
Finally, it is worth to mention possible extensions of the NLCT formalism to study out-of-equilibrium interacting bosons. Specifically, it is interesting to understand how those systems respond to quenches. Extending the NLCT-procedure would allow us to study the spreading of correlations and entanglement in the strongly interacting regime of the Bose-Hubbard model [32,33]. As the GVA possesses a well defined time-dependent extension, we expect that the structure of wavefunctions built from NLCT could help to establish a well defined extension of the non-Gaussian ansatze treated here in order to address these problems in the future.
where we have introduced the shorthand notation