Gaussian Continuous Tensor Network States for Simple Bosonic Field Theories

Tensor networks states allow to find the low energy states of local lattice Hamiltonians through variational optimization. Recently, a construction of such states in the continuum was put forward, providing a first step towards the goal of solving quantum field theories (QFTs) variationally. However, the proposed manifold of continuous tensor network states (CTNSs) is difficult to study in full generality, because the expectation values of local observables cannot be computed analytically. In this paper, we study a tractable subclass of CTNSs, the Gaussian CTNSs (GCTNSs), and benchmark them on simple quadratic and quartic bosonic QFT Hamiltonians. We show that GCTNSs provide arbitrarily accurate approximations to the ground states of quadratic Hamiltonians, and decent estimates for quartic ones at weak coupling. Since they capture the short distance behavior of the theories we consider exactly, GCTNSs even allow to renormalize away simple divergences variationally. In the end, our study makes it plausible that CTNSs are indeed a good manifold to approximate the low energy states of QFTs.


I. INTRODUCTION
Quantum field theories (QFTs) are difficult to solve out of the perturbative regime with deterministic techniques. Apart from lattice Monte-Carlo algorithms [1][2][3][4], an option would be to solve strongly coupled QFTs variationally. In a nutshell, this would mean guessing a "good" manifold M of states |ψ ν described by a manageable number of parameters ν, minimize the energy ψ ν |H|ψ ν over this class M, and hope that the answer is close enough to the real ground state |0 . As was noted by Feynman already [5], finding such a good manifold for typical QFTs is a highly non-trivial task. In particular, apart from simple Gaussian states such as free ground states, it seems impossible to have a sparsely parameterized state with easily computable local observables ψ ν |O(x 1 ) · · · O(x n )|ψ ν while keeping an extensive ansatz -the latter requirement excluding e.g. simple expansions in the particle number basis.
On the lattice, the situation has proved more favorable in the last two decades. Tensor network states (TNSs) have essentially provided what one was looking for: a sparse and extensive parameterization of many physically relevant many-body quantum states [6][7][8]. In this approach, the quantum state is obtained from low-rank tensors, contracted along the links of a network. In the translation-invariant case, all tensors are identical, making parameter economy and extensivity manifest. Tensor networks have proved successful numerically in d = 1 space dimension, with Matrix Product States (MPSs) [9], which are at the root of the earlier density matrix renormalization group (DMRG) [10][11][12], and more recently in d ≥ 2 with projected entangled pair states (PEPSs) [13][14][15][16]. As is often the case, a computationally successful method brings theoretical insights, and tensor network states have * antoine.tilloy@mpq.mpg.de allowed a succinct classification of symmetry protected [17][18][19][20] and topological phases of matter [21,22].
Given their undeniable success on the lattice, it is tempting to try to bring tensor networks to the continuum. This can in principle be done in two ways: by discretizing continuum theories to the lattice, where powerful techniques can be applied more or less out of the box, or by bringing the tensor network toolbox itself to the continuum. While the efficiency of the first approach is so far unmatched [23][24][25], we here wish to explore the second, longer term option. It turns out that bringing tensor network states to the continuum can be done rather straightforwardly in d = 1 space dimension, with the so called continuous matrix product states (CMPSs) [26] which have been applied successfully to a few QFTs [27][28][29]. Going to d ≥ 2 space dimensions has proved more difficult. Recently, a candidate higher dimensional continuous tensor network state (CTNS) was presented [30]. It is obtained as a continuum limit of a lattice tensor network state, and many of the properties of the discrete follow through to the continuum. However, the efficiency of this candidate at approximating low energy states of QFTs has not been demonstrated yet, primarily because carrying computations in the general case is substantially more difficult than with CMPSs.
Our objective here is to assess the soundness of CTNSs for the approximation of ground states of simple QFTs. To this end, we will restrict ourselves to an easily manageable subclass, the Gaussian CTNSs (GCTNSs). Naturally, this class is very restrictive, and not dense in the space of low energy states of interacting QFTs. It can approximate with arbitrary precision only the ground states of Hamiltonians quadratic in creation and annihilation operators. But understanding this class is a necessary sanity check: if GCTNSs do not even work in this simple setup, they should probably be abandoned right away. Apart from the assessment of the promise of CTNSs in general, the study of GCTNSs can also provide an economical approximation of physically relevant Gaussian states, which in theory require an infinite number of parameters (a full continuous two-point correlation function) to be defined in the context of QFTs.
We start our exploration by recalling the definition of CTNSs, characterize their Gaussian submanifold, and give basic computational tools in Sec. II. We then apply GCTNSs in Sec. III to the task of finding the ground state of a simple quadratic Hamiltonian in d = 1 and then d = 2 space dimensions. The higher dimensional setting comes with subtleties related to the infinite energy density of the ground states. Finally, we briefly study in Sec. IV a true interacting model in d = 1 space dimension, to show how the GCTNS can be efficient in some regimes even if the model under investigation is not quadratic.

II. CONTINUOUS TENSOR NETWORK STATES
In this section we define continuous tensor network states (CTNSs), explain how local observables can be computed with them, and introduce a Gaussian subset, the Gaussian CTNSs (GCTNSs) which are analytically tractable.

A. Definition
A CTNS is a quantum state |V, α , formally belonging to the Fock space F[L 2 (R d )], and defined by the functional integral [30]: where |vac is the "physical" Fock vacuum state, (ψ † ,ψ) are the canonical bosonic creation-annihilation operators on this Fock space, [ψ(x),ψ † (y)] = δ d (x − y). The "auxiliary" field φ integrated over has D components, , and ∇φ 2 := k ∇φ k · ∇φ k . This number D is the field bond dimension or simply bond dimension and is the continuous analog of the bond dimension for discrete tensor network states. We have restricted ourselves to the translation invariant case and taken the thermodynamic limit, which spares us the discussion of what happens at the boundaries. Our objective in this paper is to use this quantum state (1) as an ansatz for the ground state of a QFT Hamiltonian of interest.
Some quick comments are in order. The state is not normalized, and not all choices of functions V and α even yield a state at all (for example if V [φ] = −φ 2 ). We just assume that we choose functions such that the functional integral in (1) at least formally makes sense.
The state is parameterized by two (complex) functions, which suggests that there is an infinite number of parameters even for a number D of auxiliary fields fixed. In practice, one could expand both functions as polynomials in the fields: The maximum degrees κ V , κ α of these two expansions, together with D, then give a measure of the expressiveness of the class of states considered. Formally, the coefficients in the expansion are also tensors, and so we recover the simple idea that a tensor network state should associate a quantum state to a few elementary low-rank tensors.
Finally, we may try to give some intuition of the connection between this CTNS ansatz (1) and discrete TNSs, for the reader already familiar with the latter. A tensor network state is obtained by taking a product of elementary tensors and contracting a fraction of their indices (the bond indices) along the edges of a lattice. For CTNSs, the equivalent of the product of tensors is the exponential of the integral, the equivalent of the contraction of discrete indices is a product of integrals over auxiliary fields, which becomes a functional integral in the limit [30]. The gradient square term in (1) comes from the fact that the tensors are connected to their nearest neighbors. In this paper, understanding the derivation of CTNSs as the continuum limit of TNSs is not needed, since we will directly test the validity of CTNSs in the continuum, without relying on a discretization.

B. Generating functional
To compute expectation values of local observables on a CTNS, the most straightforward method is to introduce the generating functional Z j ,j for the normal ordered correlation functions: For example, it can be used to compute the simple twopoint function Using the Baker-Campbell-Hausdorff formula to commute the two exponentials in (2) and then using the formula for the overlap of unnormalized field coherent states, one obtains [30]: It is important to note that powers of the field in the expansion of α come multiplied and connect together the two auxiliary fields coming from bra and ket, as in a Schwinger-Keldysh functional integral. In general, if arbitrary powers of the field appear, the functional integral (4) might be diverging. Assuming that the divergences can be properly substracted, then actually computing correlation functions remains difficult non-perturbatively. Apart from Monte-Carlo techniques, a boundary CMPS method was suggested in [30] but we will not explore this further here.

C. Gaussian subset
The functional integral in (4) can be computed exactly if the expansions of V and α are truncated to quadratic and linear order respectively [31]: We call states defined by such restricted V and α Gaussian CTNS. Note that these states are Gaussian in the usual sense. Indeed, one can carry the Gaussian integral directly in the state definition (1) to get: where This expression allows to spot a lot of redundancy in the parameterization. The first and simplest observation is that V (0) does not appear because it merely changes the state normalization which we do not keep track of. Second, we notice that the second term in (7) can be incorporated into α (0) , and thus we may fix V (1) = 0 without lack of generality. This is quite intuitive: giving the auxiliary field a non-zero expectation value can be compensated by a constant source. Finally, under the mild assumption that V (2) is diagonalizable V (2) = U −1 (M/2)U , we have a straightforward rewriting: This expression could be obtained directly by taking Thus, without lack of generality, we can now assume that we have diagonal "mass" matrix M := diag(m 1 , . . . , m D ) for the auxiliary field. In the end, a GCTNS is simply parameterized by two complex vectors α (1) and m, and a scalar α (0) , that is 2D + 1 complex parameters.
We may now go back to the computation of the generating functional (4). Carrying out the Gaussian integral yields: where the operator K fulfills , and it is convenient to go to Fourier space: With this, we can compute various expectation values of the state, for example the two-point functions using (3).

D. Variational optimization
We now summarize the strategy to variationally optimize GCTNSs in practice. In what follows, we will study models specified by a local bosonic Hamiltonian: where h(ψ † ,ψ)(x) contains product of the operatorsψ,ψ † and its derivatives. For a GCTNS |V, α we introduce the associated energy density Our objective is to minimize it to find an approximation to the ground state |0 and an upper bound to the ground energy density e 0 : To carry out the minimization, one needs to be able to compute h V,α , which reduces to the computation of a sum of correlation functions ofψ,ψ † , which we know how to compute in general from the generating functional (see Appendix B). Whether we use simple gradient descent or more advanced optimization algorithm like Broy-denâĂŞFletcherâĂŞGoldfarbâĂŞShanno (BFGS) [32], we also need the gradient of h V,α with respect to the 2D + 1 complex coefficients parameterizing the state. Since we have explicit expressions for all correlation functions, this presents no fundamental difficulty and is done in appendix B.

III. A QUADRATIC MODEL IN 1 AND 2 SPACE DIMENSIONS
In this section, we present a simple quadratic, thus exactly solvable, Hamiltonian and approximate its ground state with our GCTNS ansatz.

A. The model
We first consider a model with a Hamiltonian quadratic in creation and annihilation operators and that thus has a Gaussian ground state. In fact, for a single species of spinless bosons and the usual nonrelativistic kinetic term, it is essentially the most general one can write. Such a Hamiltonian can typically be obtained as the mean field approximation of a weakly interacting Bose gas, but we take it as an exact starting point here. Another instructive way to interpret this Hamiltonian is to see it as the regularized Hamiltonian of the relativistic free boson [28]: whereπ,φ are the traditional canonically conjugate fields , and Λ is a non-relativistic momentum cutoff. This Hamiltonian H Λ fb reduces to (15) with the field mappinĝ and the parameters Closing the gap, which happens when λ/µ → f c = 1/2, is equivalent to lifting the non-relativistic regulator (m Λ). This model is exactly solvable and one finds (see appendix A) that its ground state energy density is (21) which is infinite when d ≥ 2. Consequently, in d = 1, we will be able to directly optimize the energy, whereas in d ≥ 2, we will have to renormalize away the divergent part. The corresponding two-point functions can also be computed exactly and we have e.g.

B. Variational optimization in 1 space dimension
To compute the ground state energy with our ansatz, we simply compute the energy density, its gradient with respect to the parameters, and use a standard BFGS solver to find the point yielding the minimal energy. The results are shown in Fig. 1. We observe that for parameter values of order 1 away from the gap closing (say f = λ/µ = 0.25 = f c /2), the convergence to the exact value is extremely fast in Dto the point that it is difficult to probe large values of D because of machine precision issues. As we get closer to the gap closing, the convergence becomes slower, but moderate values of D still give accurate values, even for λ/µ = 0.99f c . This is compatible with the TNS folklore that gapped systems can be precisely approximated with low bond dimension, and that larger values have to be used as we get closer to a critical point.
In QFT, one might worry that optimizing the energy does not give a fast convergence of the state itself (summarized by its two-point functions in the Gaussian case). Here, because the theory is regular (or equivalently nonrelativistic), this is not the case, and we observe a fast uniform convergence of the two-point function, at least away from the gap closing (see Fig. 2).

C. A theoretical aparté
Before we move on to the trickier d = 2 space dimensions case, it is helpful to understand better the structure of GCTNS correlation functions and compare them to the exact one (22). Using the expression for the generating functional, it is straightforward to see that C V,α (p), the Fourier transform of ψ † (x)ψ(0) V,α , is of the form where w k are the complex eigenvalues of W defined in (10) and a k are complex coefficients (see appendix B for more detail). Putting all the fractions in (23) on the same denominator shows that C V,α (p) is an even rational function of degree at most 4D. Clearly, this means that there is no chance to capture C 0 (p) exactly for a finite D, since it contains a square root. However, using the identity we have that with uniform convergence for all p as long as λ/µ < 1/2. This is the same structure as a GCTNS correlation function, except that the expansion in rational functions is truncated at order 4D in p for GCTNSs.
At short distances, p → +∞, only the first term in the expansion matters. It can be reproduced exactly already by a GCTNS with D = 1, which means the UV behavior of the QFT can be captured by the simplest non-trivial GCTNS. At long distances, p 0, the series (25) is still absolutely convergent with an error decreasing exponentially with the number of terms. Hence, for a GTCNS the error should be dominated by the infrared and at most O([2λ/µ] 2D ).
Naturally, we perform a variational optimization of the energy and not a perturbative term by term optimization of the two-point function, and as a result the error obtained in practice could scale differently. And indeed, we observe in Fig. 1, at least for small D, that the error decreases faster than naively expected.

D. Variational optimization in 2 space dimensions and renormalization
In d = 2 space dimensions, several two point functions of interest diverge when taken at equal points. In particular, the kinetic energy ∇ψ † ∇ψ andψψ +ψ †ψ † terms diverge when evaluated on GCTNSs. This can be traced back to the fact that the corresponding momentum integrals (see appendix B) diverge logarithmically. This divergence can be renormalized in a way we now explain. First, we introduce a hard momentum cutoff Λ such that correlation functions are finite. We then observe that the energy density evaluated on a GCTNS reads: such that the energy can be split into a regular and log divergent part. For the Hamiltonian (15) we consider, the log divergent part can be evaluated exactly and we find: where we used the simplified notation α (1) j = α j . Importantly, h div can be made negative and minimized exactly, yielding the condition: This condition defines a submanifold of "maximally divergent energy" GCTNS on which the parameters can be numerically tuned to minimize the remaining finite part h r . We obtained the condition (28) in a variational way, only asking that the energy be minimal and taking the cutoff to infinity. As a welcome surprise, it provides the same divergence of the energy density as the exact solution! Indeed, as can be seen from (21), the latter diverges as −λ 2 ln Λ 2 /(π), exactly as for GCTNS on the submanifold defined by (28). So not only can a GCTNS capture the UV behavior of the exact ground state, it captures it exactly upon optimization.
In what follows and for comparison, we consider only the renormalized part of the energy density h r := lim Λ→∞ h − λ 2 ln(Λ 2 )/(4π). For the exact solution it gives the "renormalized" energy density (29) which is of course finite and which we can compare to h r since the counter terms used in both cases (the divergent parts) are identical. Again, we insist that this optimization procedure and the associated renormalization of the energy density do not require knowing the exact solution. Results are shown in Fig. 3 and we observe that the convergence for the renormalized energy density and two point function is qualitatively as good as in the d = 1 case.

IV. A QUARTIC MODEL IN 1 SPACE DIMENSION
In this section we study a simple quartic model, the Lieb-Liniger model, that has a non-Gaussian ground state. Consequently, there is no hope to approximate it with arbitrary precision with GCTNS, but we may still capture qualitative features.

A. Lieb-Liniger model
The Lieb-Liniger model is about the simplest model of interacting bosons in d = 1 space dimension and is given by the Hamiltonian where c is the strength of the coupling. The number of particles is conserved and another important parameter in the model is the particle density ρ = ψ †ψ . The physics of the model depends on the adimensional coupling γ = c/ρ. This model is integrable and with the Bethe Ansatz it is possible to write an exact equation for the energy density in the ground state, which can be solved numerically to essentially arbitrary precision or expanded in a power series at weak and strong coupling [33,34].
The ground state of this model is not a Gaussian state, and as a result a GCTNS cannot approximate it with arbitrarily good precision even for large D. However, it is possible for a GCTNS to give reasonable approximation in some regime, which is what we aim to explore here. To this end, we will compare with two other simple approximation techniques: classical solution and mean field. For us, the classical solution is simply what we obtain by minimizing the energy in the space of coherent states, or equivalently GCTNS with D = 0. The meanfield approximation corresponds to the ground state of a different Hamiltonian, namely the mean field quadratic Hamiltonian of the same model. In appendix B 2, we explain how to deal with the quartic terms and how to obtain the mean field Hamiltonian.
Our analysis can be seen as the continuum analog of the one carried recently for the Bose-Hubbard model [35], where a generic Gaussian state approximation was compared with standard classical and mean field solutions. In our case, aside from dealing with the continuum, we have the refinement that we do not use the most general Gaussian states in the first place (which would anyway require infinitely many parameters), but a tower of more and more expressive submanifolds indexed by D.

B. Results
In practice, we simply minimize the energy density of the model over GCTNSs of fixed D keeping ρ = 1 fixed with gradient descent. As GCTNSs are Gaussian states, the expectation value of the quartic term is simply computed with Wick's theorem (see appendix B), and thus the energy density and its gradient are easily evaluated.
In Fig. 5, we can see that the upper bound provided by GCTNS approaches the exact ground energy as the coupling γ gets smaller. This is expected: the ground state of a weakly interacting Bose gas becomes Gaussian when the coupling goes to zero.
What is remarkable is that the simplest GCTNS ansatz for D = 1 is already sufficient to get all the expressive power of Gaussian states in this case. Almost all the improvement from the classical solution D = 0 is reached for D = 1. The refinements obtained with larger D are not necessary in the sense that they bring improvements in the energy density much smaller than the distance between the best Gaussian energy density and the true energy density. This is rather intuitive: if a Gaussian state is anyway not the exact solution, we do not gain much by getting the absolute best Gaussian state, and a crude approximation of the best Gaussian state can do qualitatively as well.

V. DISCUSSION
Let us briefly summarize our results. We were mainly interested in knowing if CTNSs had the right properties to be good trial wave functions for quantum field theories, mirroring the efficiency of their discrete lattice counterparts. To this end, we focused on a subclass of analytically tractable CTNSs, the Gaussian CTNSs, which are also a submanifold of general Gaussian states. Optimizing GCTNSs on a simple non-relativistic quadratic Hamiltonian, we obtained a very good match with the exact solution, both for the energy density and the state itself (parameterized by its two point function). Importantly GCTNSs have the right UV behavior, even for the minimal number of auxiliary fields D = 1. This allows to exactly renormalize the divergent part of the Hamiltonian density in 2 space dimensions, all variationally, without requiring any knowledge of the exact solution.
With GCTNSs of moderate field bond dimension D, it is even possible to go near the gap closing, corresponding to the relativistic limit of the model we considered. Hence, GCTNSs have exactly the right UV properties to approximate non-relativistic QFTs and they can also accommodate relativistic theories provided the cutoff scale is only reasonably far from the physical mass scale (at least up to Λ/m ∼ 10 2 ) Naturally, interesting Hamiltonians are not quadratic and thus do not have a Gaussian ground state. In this context, a more general CTNS would be required, but it is worthwhile to see if GCTNSs can help already. For the Lieb-Liniger model, we showed GCTNSs allowed to obtain good approximations of the energy density, at least in the weak coupling regime. Importantly, the lowest bond field dimension D = 1 already allows to capture essentially everything general Gaussian states (with infinitely many parameters) can capture. This means that GCTNSs allow a drastic compression of Gaussian states for quantum field Hamiltonians, yielding a potentially large gains for methods that build upon them. All these results are encouraging, and demonstrate that CTNSs indeed have the right properties expected from their discrete tensor network analogs, and, as a result, deserve to be studied further.
Promising extensions of this work are already possible while staying in the relatively easy realm of Gaussian states. The states we considered could be extended to deal with Fermions, where richer physics already appears for quadratic Hamiltonians. Dealing with multiple species of Bosons / Fermions could also enable the exploration of topological phases and see if their characterization for GCTNS matches what can be seen on the lattice. Further, in general, one can obtain much more from Gaussian states than mere ground states, and one could obtain the spectrum and real time dynamics with GCTNS extending the geometric methods developed in [35,36]. Finally, the success of tensor network methods has been well understood from their entanglement properties, and it would be useful to see if such an analysis can be done as well in the context of GCTNSs. In particular it is still unknown if the bond field dimension can upper bound the prefactor in the area law scaling of entanglement entropy, as it does in the discrete.
Going beyond the Gaussian setting to deal with genuine interacting theories could be done in different ways. A first step could be to stay with GCTNSs but considering a sum of them, which is no longer Gaussian. In this context, the fact that low field bond dimension and thus very few parameters give already good approximations of the best Gaussian states would allow to consider large sums. This would be prohibitively expensive in the more brutal approach of considering a sum of generic Gaussian states. GCTNSs could also be used to construct a better basis of states for Hamiltonian truncation methods. In this approach (see e.g. [37]), one diagonalizes an interacting Hamiltonian in a truncated basis made from the low energy sector of the free Fock space. With GCTNSs, this free Fock space could be replaced by the Fock space built from the excitations above a GCTNSs optimized on the interacting Hamiltonian.
Another, more radical option is to use genuinely non-Gaussian CTNSs. There, the difficulty is that it is not possible to compute correlation functions exactly, and in particular to compute the energy density one typically wants to minimize. As a first step, such correlations could be evaluated with Monte-Carlo or perturbatively. Note in the latter case we would still have an overall nonperturbative method: even at the lowest order of Taylor expansion, the Gaussian part would already contain nonperturbative effects in the model coupling constant. The most appealing option and the one also most in the spirit of tensor networks would be to evaluate correlation functions of a CTNS in 2 dimensions using the transfer matrix method in 1 dimension as proposed in [30]. In 1 space dimension, one can use CMPSs to efficiently find the largest eigenstate of an operator, here the transfer matrix. This would reduce the problem of computing CTNS correlation functions to that of optimizing a CMPS. This would likely require an improvement of the efficiency of existing CMPS algorithms, but does not seem out of reach. Ultimately, although a lot remains to be done to make CTNSs practically useful in the context of interacting QFTs, we hope that the present work offers evidence that this is a path worth pursuing.
Note: While finishing the present paper, we got aware of work conducted in parallel at the University of Ghent by Aelbrecht, Mortier, and Haegeman [38].
Note that in the models we are considering one can choose the gauge where a (0) ∈ R. We set because it corresponds to the C V,α when α (0) = 0. Let us now compute the real space correlation function at equal points, which is needed to get the energy density. As the momentum integral of the zero mode is trivial, we focus on the contribution of the contribution from C V,α (p) α (0) =0 . First we diagonalize matrix K(p) with an unitary 2D × 2D matrix U , such that W = U −1 LU and L is a 2D × 2D diagonal matrix with eigenvalues λ 1 , λ 2 , ..., λ 2D . Note that the matrix W needs to be positive definite for the state to be physical, thus Re[λ i ] > 0. We get and hence This allows to find the equal point 2-point function with the integral This integral is convergent in d = 1 and logarithmically divergent in d = 2. However, the divergences cancel each other in the sum (B7) as they do not depend on λ i and thus the particle density is finite in d = 2 space dimensions (and in fact even in d = 3). We can proceed in the same way to compute the other correlation functions ψ (x)ψ(x) and ψ † (x)ψ † (x) . For these, the logarithmic divergences do not cancel each other in d = 2 and contribute to the divergence of the energy density we explain how to renormalize in the main text.
To compute the kinetic energy density, we simply take the derivative of the two point function lim x→y ∂ x ∂ y ψ † (x)ψ(y) . Ultimately, this yields the same formula as before with the replacement of I 1 by This latter integral is linearly divergent in d = 1, but again this divergent part is independent of λ i and thus cancels in the expression for the kinetic energy. In d = 2, the leading divergence is quadratic and cancels in the sum but a subleading logarithmic divergence subsists in the expression of the kinetic energy density, as well as in ψ (x)ψ(x) and ψ † (x)ψ † (x) , contributing to the overall logarithmic divergence of the energy density.

Four-point function
The Lieb-Liniger Hamiltonian contains quartic terms and, as a result, evaluating its energy density requires computing a 4-point function. As our states are Gaussian, we can use Wick's theorem or the expression for the generating functional Z j ,j to get: ψ †ψ †ψψ =| ψ | 4 + 4| ψ | 2 ψ †ψ | α (0) =0 where all the operators are taken in the same point x which we omitted since the problem is translation invariant. We have also split the 2-point functions into a part that does not depend on α (0) and the zero mode contribution: The latter corresponds to the condensed fraction in the Lieb-Liniger model. Taking the mean field approximation is equivalent to neglecting the last two lines in (B8) as one assumes that the zero mode ψ dominates.