The Inhomogeneous Gaussian Free Field, with application to ground state correlations of trapped 1d Bose gases

Motivated by the calculation of correlation functions in inhomogeneous one-dimensional (1d) quantum systems, the 2d Inhomogeneous Gaussian Free Field (IGFF) is studied and solved. The IGFF is defined in a domain $\Omega \subset \mathbb{R}^2$ equipped with a conformal class of metrics $[{\rm g}]$ and with a real positive coupling constant $K: \Omega \rightarrow \mathbb{R}_{>0}$ by the action $\mathcal{S}[h] = \frac{1}{8\pi } \int_\Omega \frac{\sqrt{{\rm g}} d^2 {\rm x}}{K({\rm x})} \, {\rm g}^{i j} (\partial_i h)(\partial_j h)$. All correlations functions of the IGFF are expressible in terms of the Green's functions of generalized Poisson operators that are familiar from 2d electrostatics in media with spatially varying dielectric constants. This formalism is then applied to the study of ground state correlations of the Lieb-Liniger gas trapped in an external potential $V(x)$. Relations with previous works on inhomogeneous Luttinger liquids are discussed. The main innovation here is in the identification of local observables $\hat{O} (x)$ in the microscopic model with their field theory counterparts $\partial_x h, e^{i h(x)}, e^{-i h(x)}$, etc., which involve non-universal coefficients that themselves depend on position --- a fact that, to the best of our knowledge, was overlooked in previous works on correlation functions of inhomogeneous Luttinger liquids ---, and that can be calculated thanks to Bethe Ansatz form factors formulae available for the homogeneous Lieb-Liniger model. Combining those position-dependent coefficients with the correlation functions of the IGFF, ground state correlation functions of the trapped gas are obtained. Numerical checks from DMRG are provided for density-density correlations and for the one-particle density matrix, showing excellent agreement.


Introduction
Most gapless 1d quantum systems fall into the Luttinger liquid universality class, an effective field theory approach that accounts for their low-energy (or large distance) excitations [1][2][3][4][5]. This paradigm is well known for being intimately related to certain 2d conformal field theories (CFT) [6] with central charge c = 1, namely free massless boson theories at different compactification radii, that are themselves at the heart of the Coulomb gas picture of 2d statistical models developed in the 1970s and 1980s [7][8][9][10]. Nowadays, those free theories are playing a fundamental role in modern mathematics, especially at the intersection of probability theory and conformal geometry, where they are known as the "Gaussian Free Field" (GFF) 1 [11].
In this paper, our goal is to explore the case where the assumption of a constant parameter K is relaxed. This is natural in many physically relevant situations. Perhaps the most notable example is that of a 1d gas of bosons, modeled by the Lieb-Liniger model [34], trapped in an external potential V (x), where x is the spatial coordinate. In this model, the Luttinger parameter K acquires a spatial dependence, As we will explain shortly, contrary to the case of constant Luttinger parameter K, the underlying field theory is no longer a GFF. Instead, it is an "inhomogeneous" generalization of the GFF, with a spatially varying coupling constant, which we will dub "Inhomogeneous GFF" (IGFF). Because the IGFF is a free (or Gaussian) theory, calculating correlation functions in the IGFF boils down to solving some boundary value problem by calculating its Green's function. This will be discussed in full detail in Sec. 2. In Sec. 3, we will apply that formalism to calculate ground state correlation functions in the trapped Lieb-Liniger model.
In the rest of this introduction, we explain how exactly this work differs from previous ones on inhomogeneous Luttinger liquids, and then motivate the introduction of the IGFF, defined by the action (1.8) below.

Relation with previous works on inhomogeneous Luttinger liquids
Over the past twenty years, some of the results we will derive or use in this paper have been partially reported in the literature. Here, we give a brief account of the existing works that aimed at the same direction, to the best of our knowledge. In 1995, Maslov and Stone [35] and (independently) Safi and Schulz [36] investigated the Landauer conductance of an interacting electron wire. Both ends of the wire are connected to a lead, represented by free electrons. In that setup, the Luttinger parameter jumps from K = 1 in the leads to some value fixed by the interactions in the wire. So does the velocity v of gapless excitations, jumping from the Fermi velocity in the leads to some other value in the wire. Thus, the problem of calculating reflection and transmission coefficients reduces to studying the Luttinger liquid Hamiltonian with K(x) and v(x) that are step functions. To our knowledge, this is the first occurence of an "inhomogeneous Luttinger liquid" with non-constant Luttinger parameter K(x). It turns out that, in this particularly simple setup, the Green's functions can be expressed analytically. Maslov and Stone [35] used a Lagrangian formulation and therefore wrote the action of the IGFF (1.8) -see Eq.
(3) in their paper -; to our knowledge, this is the first time that action appeared in the literature. Maslov and Stone also derived a differential equation for the propagator (Eq. (6) in their paper) that is similar to the generalized Poisson equation from Sec. 2 below. The same model was studied by Fazio, Hekking and Khmelnitskii [37] in the context of thermal transport. However, the physical quantities studied in Refs. [35][36][37] were simply defined in terms of integrals of the propagator, so the authors did not have to push further the calculation of more general correlation functions.
About a decade later, in 2003, Gangardt and Shlyapnikov [38] had similar insights, and wrote the Hamiltonian of the inhomogeneous Luttinger liquid (see the equation above Eq. (12) in their paper), this time with the purpose of computing correlation functions of a 1d Bose gas trapped in an harmonic potential. They took the Luttinger liquid Hamiltonian [2], assumed that K and v were both position-dependent, and then used the Local Density Approximation (LDA) to fix these parameters. They extracted K(x) and v(x) from the Bethe ansatz solution of the homogeneous Lieb-Liniger model (see also Ref. [39] where LDA was used to calculate local correlation functions). This is exactly what we will do in Sec. 3 below. From there, they derived an expansion of the boson field which, in principle, allows to compute correlation functions. The same logic was followed by Ghosh in 2006 [40] and by Citro et al. in 2008 [41]. Some of these results have been reviewed in Ref. [22] (section V.E).
The same kind of approach was also developed in the context of multicomponent 1d Fermi gases. In 2003, following the spirit of [42], Recati et al. [43] investigated the spectrum and discussed experimental realizations of spinful ultra-cold Fermi gases; independently of [38], the authors assumed that the space-dependent parameters K(x) and v(x) could be fixed by LDA. This idea was later used by Liu et al. [44,45] to study the phase diagram of the 1d Hubbard model. As far as we are aware, this has not been explicitly used to compute correlation functions in this context, see Ref. [46] for a review.
The innovation of the present paper, compared to Refs. [35][36][37][38][39][40][41], is twofold. First, in Sec. 2 we discuss the IGFF and its correlation functions in full generality. To our knowledge, such a general and complete discussion has not appeared elsewhere, and it should be useful to some readers. Second, we believe that an important ingredient has been missed in Refs. [38,40,41], and that the results for correlation functions reported in those references are, in fact, not entirely correct. The reason is the following.
In general, local observables in a microscopic modelÔ(x) (say, the Lieb-Liniger model) are related to field theory operators φ(x) only through nonuniversal coefficients C. To elaborate, observablesÔ(x) are expected to have expansions of the formÔ where the sum in the r.h.s. runs over all possible local operators φ j in the field theory, and the non-universal coefficients C (Ô) j are dimensionful numbers. As usual, such an expansion is to be understood as a statement about correlation functions: correlations functions in the microscopic model are related to the ones of the field theory, providing asymptotic expansions of the former in the limit where all the points are well separated, In homogeneous systems, the non-universal coefficients merely contribute as global prefactors in the correlation functions (a useful and detailed discussion of those coefficients can be found in Refs. [47][48][49][50]). But, in inhomogeneous situations, those dimensionful coefficients C (O) j are themselves position-dependent, j (x), so they have a crucial impact on the correlators. This point seems to have been overlooked in previous works, see Fig. 8 in App. B for a plot comparing our result to the case where these coefficients are omitted.
In this paper, we use LDA to fix those dimensionful coefficients. We illustrate this in Sec. 3 in the Lieb-Liniger model. The prefactors are extracted from form factors formulae derived in the 1990s by algebraic Bethe ansatz [51][52][53], see App. B for more information. The method is then checked against numerical results in the Lieb-Liniger model obtained from DMRG, using the C++ library ITensor [54], see App. D for details about the simulation. The agreement is quite impressive, as can be seen in Figs. 5, 6 and 7 below. Figure 1: Cartoon illustrating the separation of scales in a trapped 1d gas. The typical length L on which the local chemical potential µ(x) = µ − V (x) varies is of the order of the total size of the system. This macroscopic scale is much larger than the microscopic scale corresponding to the inverse density ρ −1 . There exists a mesoscopic scale at which the system consists of fluid cells that are locally homogeneous, but still contain a very large number of particles.

The underlying assumption: separation of scales
The approach we adopt in this paper is valid in the limit where the system exhibits separation of scales, see Fig. 1. This is the limit where the confining potential V (x), and more generally all local thermodynamic quantities of the quantum gas -such as its particle density, energy density, momentum density, etc. -vary very slowly on the microscopic scale. That microscopic scale is naturally given by the inverse density ρ(x) −1 , so the condition that the density varies slowly reads The r.h.s. defines a macroscopic scale L, which is typically of the order of the length of the system. When the macroscopic scale is much larger than the microscopic one, there exists an intermediate -or mesoscopic -scale such that ρ Then a "mesoscopic fluid cell" of size is both homogeneous (because it is small compared to the scale L at which inhomogeneity becomes important) and contains a thermodynamically large number of particles (because it is large compared to ρ −1 ). This is the key assumption that underlies the Local Density Approximation used in Refs. [38][39][40][41], and more generally all hydrodynamic approaches [55] (LDA itself being nothing but a "hydrostatic" approach [42]).
The assumption is of course also a requirement for any effective field theory approach, because the fields themselves are supposed to describe coarse-grained versions of the microscopic degrees of freedom, and this makes sense only if there exist locally homogeneous cells over which coarse-graining can be performed.
In Sec. 3, we will explain in detail what limit we take in the trapped Lieb-Liniger model, and we will see that separation of scales holds exactly in our setup. The method we explore in this paper (which extends the previous results of Refs. [35][36][37][38][39][40][41]) should then be interpreted as a way of writing asymptotic expansions of the correlation functions in the N → +∞ limit, including not only the leading order, but also the first few finite-N corrections.

The effective Gaussian action
To conclude this introduction, we explain why the problem of a quantum gas of particles in a trap leads to the IGFF, defined by the action (1.8). The content of this subsection is very similar to arguments given in Refs. [17,18]; we repeat those here only for completeness.
There are several ways of showing the connection between Luttinger liquids and the GFF (homogeneous or not). Here we give an argument that is particularly short and is a good introduction to Secs. 2 and 3. More standard introductions can be found for instance in Refs. [4,5].
The argument consists of two steps.
Mapping on configurations of a height fieldin 1d, configurations of indistinguishable particles can be represented by an height function h(x) viâ namely h(x) is a real-valued function that is piecewise constant and jumps by 2π at the position of a particle. It is defined only up to a constant shift, h → h + const. To calculate ground state correlation functions, it is useful to imagine that the system evolves in imaginary time, and focus on correlation functions at arbitrary points (x, τ ) in spacetime, and then later specify that all points are taken at imaginary time τ = 0. For instance, the two-point function of the height field would be where H is the Hamiltonian of the system, and β is the inverse temperature, that one sends to zero. The fluctuating field h(x, τ ) is then viewed as a function on 2d spacetime.
Action for the height fieldthe second step consists in writing an action for the height field h(x, τ ). Our choice for the action is guided by two basic observations. First, assuming that the underlying microscopic model is described by a local Hamiltonian H, the action should be local. Second, physical observables should be invariant under constant shifts h → h + const., so the action must also possess that symmetry. This leads us to the general form where the dots stand for higher order derivatives. The Lagrangian density L cannot depend on h(x, τ ) itself, only on its derivatives, because we are asking that it is invariant under constant shifts. Finally, assuming that the action S[h] is minimized by a unique classical configuration of the field, h cl. , we can expand to second order around that minimum, In 2d, higher order terms have scaling dimensions larger than 2 and are RG irrelevant; we can therefore discard them. The only free parameters of the effective theory are then the three independent real components of the Hessian ∇ 2 L at h = h cl. , which is a positive 2 × 2 symmetric matrix that typically depends on position. It is convenient to interpret the inverse of that matrix as an emergent metric on spacetime and to rewrite the Gaussian action as . This is the action of the IGFF. It is the most general action for the height field h that is both local and invariant under constant shifts. We now study this theory in greater detail.

The 2d Inhomogeneous Gaussian Free Field
This section is devoted to the 2d Inhomogeneous Gaussian Free Field, which is the mathematical object that underlies inhomogeneous Luttinger liquids. The IGFF is a rather straightforward generalization of the Gaussian Free Field, parametrized by a function K : Ω → R >0 that represents the position-dependent coupling strength in the action S[h] = 1 8π Ω √ gd 2 x K(x) g ij ∂ i h∂ j h. The usual GFF is recovered when the function K is constant, simply by rescaling the height field h → √ Kh. In other words, while the usual GFF depends only on the domain Ω and on the conformal class of the metric g [11], the IGFF also depends on the function K.
For simplicity, we work on a simply connected, open subset Ω ⊂ R 2 . [Later, when we will apply the IGFF to inhomogeneous Luttinger liquids, Ω will be identified with spacetime.] Since every metric in 2d is conformally flat, and because the action (1.8) is invariant under Weyl transformations g ij → e 2σ g ij , w.l.o.g. we can work in the Euclidean metric such that the action of the IGFF becomes Here h is a real-valued function on the closure Ω, with Dirichlet boundary conditions,

Propagator of the IGFF, generalized Poisson equation
Correlation functions can be defined as path integrals, where we integrated by parts in the last line. Thus, the 2-point function is identified with the Green's function of a generalised Poisson operator ∇ · 1 is symmetric under exchange of x and x , and solves the linear differential problem (2.6) The superscript 'D' refers to the boundary conditions (Dirichlet), while the subscript [K] emphasizes the fact that the IGFF is parameterized by the function K : Ω → R >0 . Contrary to the GFF, where the Green's function is easily obtained by conformal mapping of the domain Ω onto the upper half-plane (leading to explicit formulas in a number of physically relevant problems), no such explicit expression is available in general for the IGFF. The Green's function of the generalized Poisson operator can, however, be efficiently calculated numerically. We note that the generalised Poisson operator is well-known from classical electrostatics [56]: it appears in the equation satisfied by the electrostatic potential V (x) in the presence of a spatially-varying dielectric constant ε(x): ∇ · ε(x)∇V (x) = 0. The analogy with electrostatics will be pushed further below.
In summary, n-point correlation functions h(x 1 ) . . . h(x n ) in the IGFF are all expressible in terms of the Green's function of a generalised Poisson operator; notice that the result is divergent when x i → x j , because the Green's function has a logarithmic singularity.
In applications to inhomogeneous Luttinger liguids, we need a few additional results about the IGFF, which provide natural generalisations of the ones that are well known for the GFF. First, we need to deal with vertex operators, which requires that we make sense of correlation functions of several insertions of h(x) at the same point. This is what we do in the next subsection. Second, we need to compactify the field h (meaning that we must view h as taking values in the circle R/2πZ instead of the real line R), which we do in subsequent subsections.

Correlations at equal points, regularized Green's function
As usual in field theory, one needs a regularization procedure to make sense of multiple insertions of the field h(x) at the same point, h n (x), n ≥ 2. This is provided by the normal order, noted : h n (x) :, which is conveniently defined as follows. For n = 0, : 1 : = 1, and for n = 1, and then, by induction on n, one defines : h n (x) : as (2.7b) The second term is introduced to cancel the divergence of the Green's func- With that definition, the expectation value : h n (x) : is finite, and is equal to This regularized Green's function will appear many times in the following. Notice that we use almost the same notation as for the Green's function itself, , but with a single argument instead of two.

Vertex operators, analogy with electric charges
Exponentials of the field h(x) define vertex operators, as in the usual GFF, Correlation functions of vertex operators can be computed directly from their definition, using Wick's theorem; this is a standard exercise of combinatorial nature which we leave to the reader. The result is In the literature, such vertex operators are sometimes referred to as electric charges, in analogy with 2d electrostatics [56]. A simple way of seeing the analogy is to interpret the expectation value of i∇h(x), as an electric field E(x) in the plane. [Here, the factor i is cosmetic; it is inserted in order to cancel the one in the exponential that defines the vertex operator, such that the expectation value of i∇h(x) is real.] In the presence of vertex operators, i∇h(x) acquires a non-zero expectation value, Thus, E(x) satisfies Maxwell's equations in a medium with dielectric constant ε(x) = 1/K(x), with pointlike electric charges at positions x j , (2.12) The first equation is the Gauss' law for the displacement field ε(x)E(x), and the second one is the Faraday's law (in the absence of magnetic flux through the plane) which is automatically satisfied here because E(x) is a gradient.
In fact, the logarithm of the correlation function (2.10) is nothing but the electrostatic energy of those pointlike electric charges, in the domain Ω with a local dielectric constant ε(x) = 1/K(x), surrounded by a perfect conductor (corresponding to the Dirichlet boundary condition), 1 8π The second term is of course the Coulomb interaction for all the pairs of particles, while the first term is the electrostatic energy of each independent particle that arises from its interaction with the medium and with the perfect conductor at the boundary. Notice that the integral in the l.h.s. needs to be properly regularized to recover the regularized Green's function in that first term.

Compactification of the height field, magnetic operators
So far, we have assumed that the height field h(x) was real-valued. From now on, we compactify the target space, and h(x) is viewed as a point in R/2πZ instead of R. This compactification has two important consequences on the theory. The first consequence is the quantization of electric charges: in order to be well-defined, the vertex operator V α (x) = : e iαh(x) : must be invariant under h → h + 2π. This implies that α is an integer.
The second consequence is the appearance of a new type of local operator O β (y), representing a puncture at point y ∈ Ω, around which the field h(x) has non-zero winding: h(x) jumps by 2πβ, for some integer β, when x is dragged around the puncture counterclockwise. In other words, where C y is a small oriented contour enclosing the point y. This identity holds when inserted inside correlation functions, e.g.
Due to Dirichlet boundary conditions that impose that the contour integral along the boundary ∂Ω vanishes, ∂Ω dx · ∇h(x) = 0, the set of operators O β1 (y 1 ), . . . , O βm (y m ) inserted inside a non-vanishing correlator must satisfy the neutrality condition The operators O β are often called "magnetic operators" in the literature. Again, this is an explicit reference to the electrostatic analogy. Indeed, the equations satisfied by the "electric field" E(x), namely the expectation value of ∇h(x) (here we drop the cosmetic i from the previous subsection, because the expectation value of ∇h(x) is real) with insertions of those operators, are: (2.18) The first constraint is the equation of motion for h(x) derived from the action (2.2). Again, we view it as the Gauss law in a medium with dielectric constant ε(x), this time without electric charges. The second is just a rewriting of Eq.
(2.15) using Stokes' formula, and we regard it as the Faraday law, imagining that the plane is transpierced by infinitely thin, constantly increasing, magnetic fluxes at positions y j .

Correlation functions of magnetic operators from electricmagnetic duality
We now turn to the calculation of correlators of magnetic operators. Again, such correlators are defined as a path integral where the path integral in the numerator is over functions h d from the punctured domain Ω\{y 1 , . . . , y m } to R/2πZ that have the correct winding β j around each puncture y j , We refer to those as "height configurations with defects". In this subsection (and only here), we use a subscript 'd' for configurations with defects. The denominator in Eq. (2.19) is the path integral on configurations without defects, namely the partition function of the IGFF on Ω.
The numerator can be evaluated by separating the configurations with defects into a classical part that satisfies the equation of motion, and a quantum, or fluctuating, part: is assumed to satisfy the equation of motion, the action splits,  23) and the remaining task is to calculate the integral where h 0 d (x) satisfies the constraint (2.20), the equation of motion ∇· 1 K(x) ∇h 0 d (x) = 0, and possesses Dirichlet boundary conditions.
The electrostatic analogy provides an elegant way of calculating that integral. Indeed, the integral is nothing but the electrostatic energy 1 8π , for an electric field created by constantly increasing fluxes that pierce the plane. If we could trade those magnetic fluxes for pointlike electric charges, then the answer would be given by Eq. (2.13).
This can be done by electric-magnetic duality. If we define a new fieldẼ with components In the last line, notice that we have replaced the superscript 'D' by 'N'. This is because Dirichlet boundary conditions are dual to Neumann boundary conditions. To see this, one can think of E as ∇h d , and ofẼ as the gradient ∇h of some other functionh. Because h d = 0 at the boundary ∂Ω, the component E that is tangential to the boundary vanishes. SinceẼ is obtained from a π/2-rotation of E, this implies that the normal componentẼ ⊥ vanishes. Hence, the dual fieldh has Neumann boundary conditions, instead of Dirichlet.
In summary, the result for the correlation function of magnetic operators is (2.27) where the Green's function (as well as its regularised version, defined exactly as in Eq. (2.9) above) is the one of the generalized Poisson operator ∇ · 4K∇, with Neumann boundary conditions. This Green's function G N [1/4K] (y, y ) is symmetric under exchange of y and y , and it solves the linear differential problem , In some applications of the IGFF, one expects that we will need correlation functions involving both electric and magnetic operators. Once again, the electrostatic analogy provides a convenient way of evaluating such "mixed" correlators. Indeed, the result must take the form such that its logarithm is the total electrostatic energy of a configuration of n pointlike electric charges and m punctures with insertions of fluxes. This total energy is a sum of n + n(n−1) terms for each pair of those. We have already encountered all those terms in previous subsections. The new n × m terms here are the ones that correspond to the energy of an electric charge at position x k in the electrostatic potential created by a magnetic flux inserted at y l .
This potential, which we call F D, , is a function of x and y with value in R/2πZ that satisfies a number of constraints, which we detail now. First, we need to choose a continuous function f : (2.30) Notice that, as a consequence, for any contour C y that encircles y.
It is important to stress that, while F D,N [K,1/4K] (x, y) depends on the choice of the function f , the correlation function (2.29) does not. Indeed, imagine that we have two functions f 1 and f 2 with winding number one, and that we look at the corresponding F 1 (x, y) and F 2 (x, y) defined by Eqs. (2.30). Then F 1 (x, y l )− F 2 (x, y l ) is a continuous function for x ∈ Ω, with no winding anywhere, that Then, summing over l from 1 to m and using the neutrality condition (2.16), one sees that Thus, it has to vanish everywhere. So the correlation function (2.29) is independent of the choice of f as claimed.
It is interesting to note that, when viewed as a function of y, F D,N [K,1/4K] (x, y) satisfies a set of constraints that are dual to Eqs. (2.30): (2.32) wheren y andt y are two unit vectors respectively normal and tangent to the boundary ∂Ω at position y. This is more easily seen by considering a discrete version of the compactified IGFF, in analogy with lattice electrostatics, see App. C.

Mixed electric-magnetic operators
Finally, it might also be convenient to deal directly with vertex operators that possess both an electric and a magnetic charge. The latter are obtained when one fuses an electric operator with a magnetic one, meaning that one takes the limit x, y → z in correlation functions involving O β (y) and V α (x). It is therefore convenient to introduce a new notation for vertex operators that carry both an electric and a magnetic charge: with two indices for the two charges, such that the previous "pure electric" or "pure magnetic" operators correspond to V α (x) = V α,0 (x) and O β (y) = V 0,β (y) respectively. The correlation function of such operators can be obtained by taking m = n and x i , y i → z i in Eq. (2.29). The result is Here the regularized function F D,N [K,1/4K] (z) is defined as

Application to the Lieb-Liniger model in a trap
We now turn to the problem of calculating correlation functions of trapped 1d Bose gases. This will illustrate how the machinery of the IGFF developed in Sec. 2 is useful in practice.
We focus on the Lieb-Liniger model of spinless bosons with repulsive delta interaction, defined by the Hamiltonian whereΨ † (x) (Ψ(x)) is the boson creation (annihilation) operator that satisfies the canonical commutation relation [Ψ(x),Ψ † (x )] = δ(x − x ), m is the mass of a boson, g = ḡ > 0 is the interaction strength, µ is the chemical potential and V (x) is a trapping potential. There are two main reasons for focusing on this Hamiltonian: it is the model that is experimentally relevant to describe Bose gases through the whole range of repulsion strength in one dimension [57], and, in the homogeneous case V (x) = 0, it is exactly solvable by Bethe Ansatz (for an introduction to the Bethe Ansatz solution of the Lieb-Liniger model, see e.g. Ref. [58]). Throughout this section, we consider that m, µ, V (x) andḡ are fixed parameters, and we focus on the limit → 0. Our goal is to study correlation functions in the ground state of H, and to understand how to get the first few terms of their asymptotic expansion in in that limit.

The limit → 0
Taking the limit → 0 while keeping all other parameters fixed is a particularly convenient way of taking the thermodynamic limit N → +∞. The reason is the following.
In the homogeneous case V (x) = 0, dimensional analysis shows that the particle density in the ground state must take the form for any positive value of the chemical potential µ > 0, where F (.) is some function that can be calculated from Bethe Ansatz. Thus, at least in the homogeneous case, the density of particles diverges as 1/ when → 0. Then, in the inhomogeneous case, one can rely on the following self-consistent argument. Assuming that the density of particles is sufficiently large at each point where the local chemical potential µ(x) = µ − V (x) is positive, one can rely on separation of scales, see Fig. 1. Under this assumption, the density is This is the Local Density Approximation. It shows that the density locally diverges as 1/ at every point where µ − V (x) > 0, thus separation of scales (see Fig. 1) becomes exact in the limit → 0.
Since the total number of particles is the integral of the density ρ(x) over the region where µ − V (x) > 0, it is clear that N ∝ 1/ , so that limit is a thermodynamic limit, as claimed. Importantly, in our setup, the local dimensionless parameter    Now, to apply the formalism of Sec. 2 to the trapped Lieb-Liniger gas, we need to fix the Luttinger parameter K and the (conformal class of the) metric g on the domain Ω = [−R, R] × R. To do this, we rely once again on separation of scales, and we use the exact solution from Bethe Ansatz that is available in the homogeneous case.
Thanks to separation of scales, we can imagine that we focus first on correlation functions within a single mesoscopic fluid cell, see Fig. 1. The mesoscopic cell is homogeneous and contains a thermodynamically large number of particles, so the correlation functions must be exactly the same as the ones of the homogeneous system in the thermodynamic limit. But, in the homogeneous problem, both K and g are known, and the metric is simply (with x = (x, τ )) (3.7) By dimensional analysis, the effective velocity v of gapless excitations above the ground state is of the form v = v F f 1 (γ) for some function f 1 , where v F = 2µ/m is the Fermi velocity. Similarly, the (dimensionless) Luttinger parameter K is of the form f 2 (γ). These functions f 1 and f 2 are known from Bethe Ansatz; they are plotted in Fig. 3. This fixes the metric g and the parameter K within each mesoscopic cell. Then, of course, the action of the IGFF in the entire domain Ω = [−R, R] × R is determined.
So, to sum up, we know what field theory needs to be solved: it is the IGFF in the metric (3.7), with a velocity v and a Luttinger parameter K that both depend on the position x through the local density ρ LDA (x) and the local dimensionless interaction parameter γ(x). Correlation functions can thus be expressed in terms of the Green's functions G D [K] and G N [1/4K] defined in Sec. 2, which are efficiently calculated numerically.
We conclude this subsection with a short remark about the coordinate system. In Sec. 2, we relied on a system of isothermal coordinates to simplify the expressions associated with the differential operators -the generalized Poisson operators -whose Green's functions appear in the IGFF correlators. Here, a system of isothermal coordinates is readily available [14]. Indeed, one can stretch the spatial coordinate x according tõ As a consequence, correlation functions can be written directly with the formalism of Sec. 2, by working in stretched coordinates x = (x, τ ). To get expressions of correlators in the physical coordinates (x, τ ), one simply has to keep track of Weyl factors: under the Weyl transformation g → e 2σ g, a local operator φ(x) with scaling dimension ∆ transforms as φ(x) → e −σ∆ φ(x). For instance, the two-point function of φ could first be calculated in the coordinate system (x, τ ) using the formalism of Sec. 2, and then be rewritten as

Expansion of the density operator
To relate correlation functions of a microscopic observableÔ(x) to the ones in the IGFF, we need to find an expansion of the form (1.1) forÔ(x) in terms of local operators in the field theory. This is what we do now, for the local densitŷ ρ(x). As in Sec. 2, we view the operators as evolving in imaginary time τ . So they are functions of the coordinate x = (x, τ ), and we will take τ = 0 at the end of the calculation, to get equal-time ground state correlations. The local operators in the IGFF are the derivative of the height field ∂ x h and the mixed electric-magnetic operators V α,β (x). Operators with non-zero magnetic charge β = 0 cannot appear in the expansion of the local densityρ(x, τ ), because they correspond to creation/annihilation processes at point (x, τ ); those will be discussed in more details in Sec. 3.5 below. So, the local density must have an expansion of the form where the C (ρ) p,0 are dimensionful coefficients that we need to determine, and the "descendents" terms correspond to derivatives of the local operators, which are less local and generate subleading corrections to correlation functions. For simplicity, in this paper we will discard them and keep only the terms p = ±1 in the sum: (3.12) Our task is now to identify the dimensionful coefficients C , where ρ is the particle density, and A(γ) is a real positive function of the dimensionless interaction parameter γ. It is also known (see e.g. Ref. [5]) that the phase of the coefficient C ±1,0 is e ±i2kFx where k F = π ρ is the Fermi momentum. The function A(γ) can be calculated from Bethe Ansatz, and is plotted in Fig. 4 (see App. B for details). Since the coefficients C (ρ) ±1,0 should depend only on the local properties of the gas, the expression found in the homogeneous case must remain valid also in the inhomogeneous case, replacing ρ and γ by ρ LDA (x) and γ(x). We then arrive at the expansion of the density operator Here, to lighten the notations, we write A(x) and K(x) instead of A(γ(x)) and K(γ(x)). The phase ϑ(x) is a WKB phase, given by (3.14)  [51][52][53]. See e.g. Refs. [47][48][49][50]59] for details on how to extract such coefficients from form factors. The dashed lines correspond to known asymptotics: It is obtained by requiring that ∂ x ϑ(x) equals the local Fermi momentum k F (x) = πρ LDA (x); the additive constant π 2 is fixed by an exact calculation in the free fermion case (i.e. the Tonks-Girardeau limit γ → ∞), see App. A.

Density profile, and density-density correlation
We now have all the ingredients that are necessary to calculate correlation functions of the local densityρ(x). Taking the expectation value of the r.h.s. in Eq. (3.13), and using the results of Sec. 2, one finds This follows from the fact that ∂ x h = 0 and V 1, , see Eqs. (3.10) and (2.34). In Fig 5, we compare this result to a direct DMRG simulation of the Lieb-Liniger gas. The agreement is excellent. [Another highly non-trivial check for formula (3.15) is the fact that, in the Tonks-Girardeau limit γ → +∞ and in a harmonic trap, the result is an exact match to the one obtained by evaluating the large-N asymptotics of the Hermite kernel, see App. A for details.] The oscillations of the density are well reproduced by the first subleading corrections from Eq. (3.12) and are usually interpreted as Friedel oscillations [61,62].
Next, we use the expansion (3.13) and the formulae of Sec. 2 to evaluate density-density correlations. [For a study of density-density correlations in the  Fig. 2). We see that the density profile obtained by including the first IGFF correction is in excellent agreement with the exact numerical profile, and that the Friedel oscillations are correctly reproduced at the edge of the trap. Bottom row: to show that the excellent agreement is not restricted to the case of harmonic potentials, we also display the density profile for the Tonks-Girardeau gas (i.e. γ → +∞) in a double-well potential.
homogeneous case, see e.g. Ref. [63].] We find, for the connected part, where x = (x, τ ) and we set τ = τ = 0. In Fig. 6, we display a comparison with the density-density correlation obtained from DMRG, as a function of x, for two positions x = 0 and x = −0.5R. Again, the agreement is excellent.   Fig. 2).

The one-particle density matrix
Finally, we will apply the IGFF to the computation of the one-particle density matrix Again, the first step consists in identifying the most relevant field theory operators that appear in the expansion of the creation and annihilation operatorŝ Ψ † (x) andΨ(x). Here, for simplicity, we restrict ourselves to the leading order, which is given by a single magnetic vertex operator, [Subleading terms will be investigated elsewhere.] The coefficient C 0,1 is identified in the same manner as for the density operator: we start by considering the case of (homogeneous) mesoscopic fluid cells, then go to the inhomogeneous case relying on the LDA.
Given that the creation/annihilation operator has dimension 1/2, and that the magnetic vertex operator has scaling dimension 1/4K, the amplitude of the coefficient must take the form C B(γ) for some function of the dimensionless interaction parameter B(γ). This function B(γ) is again calculated using form factors formulae, see Refs. [51][52][53]64] and Fig. 4. When going to the inhomogeneous case, we know from the homogeneous solution that the coefficient C (Ψ) 0,1 (x) does not have a space-dependent phase but can only have a global constant phase, which we can fix to zero, such that B(γ) is real and positive. We then have where we write B(x) instead of B(γ(x)).  The rest of the calculation is straightforward. We use formula (2.27), and, taking into account the Weyl factors (3.10), we obtain , (3.20) with x = (x, τ ), and τ = τ = 0. In Fig. 7, we check this formula against DMRG. Even though we only considered the leading order here, we find very good agreement.

Conclusion
The purpose of this paper was to develop the formalism of the IGFF, and provide an exhaustive study of correlations of local observables in that theory. We did so in Sec. 2. Then, in Sec. 3, we explained how, in practice, this formalism gives access to correlation functions of inhomogeneous systems, by focusing on ground state correlations of the Lieb-Liniger gas in a trapping potential.
To conclude this paper, let us mention four directions which, in our opinion, would deserve further investigation.
• in cold atoms experiments, where some correlation functions are measurable [22,[27][28][29]33], the gas is at finite temperature. Therefore, it would be very interesting to generalize the results of this paper to finite temperature. Usually, in field theory, working at finite temperature is relatively easy: one simply needs to compactify the imaginary time direction. However, there could be issues related to the boundary of the system: how to properly describe the fluctuations of the particles near the edge of the gas? Those won't be obtained simply by compactifying the time direction in the field theory. It would also be interesting to make the connection with other recent works on trapped 1d quantum gases at finite temperature, for instance Refs. [65][66][67].
• as mentioned in Sec. 1.1, the inhomogeneous Luttinger liquid also appears in the context of multi-component 1d Fermi gases. Motivated by recent experimental advances [68], it would be interesting to extend the results of Sec. 3 to the case of SU(N ) systems. In this case, the integrable model of interest (the one that replaces the Lieb-Liniger model) would be the Gaudin-Yang model [46].
• another natural extension of this work would be to tackle time-dependent problems. In the static case studied here, a key role is played by LDA, or hydrostatics, to fix the parameter K and the background metric g in the effective field theory. In a dynamical situation, for instance a breathing Lieb-Liniger gas in a trap [69][70][71][72], those parameters would have to be extracted from an hydrodynamic approach. It would be interesting to study how this works in practice, starting with the zero-temperature case. We note that, in a very inspiring paper, Abanov [73] has studied a related problem in imaginary time (see also Ref. [15], on a similar imaginary time problem).
• finally, perhaps the most challenging problem is to understand whether it is possible to have a more general theory of fluctuations and correlations in the recently developed theory of Generalized HydroDynamics (GHD) [74,75]. So far, the IGFF approach discussed here models only fluctuations of the particle density at zero temperature (and therefore corresponds only to a particular case in the more general GHD framework, dubbed "zeroentropy GHD" in Ref. [76]). In GHD, not only the particle density is expected to fluctuate, but all densities of conserved charges. Perhaps such a theory could take the form of a "fluctuating hydrodynamics" in the spirit of Ref. [77], or perhaps a multi-component version of the IGFF (possibly with arbitrarily large number of components). A step towards correlation functions in GHD has been taken very recently by Doyon in Ref. [78]; it would be a good starting point to understand if/how his results connect to inhomogeneous Luttinger liquids.

Acknowledgements
We The DMRG simulations were performed using the open-source ITensor library [54] .

A The Tonks-Girardeau limit
The Tonks-Girardeau (TG) regime is the limit of hard-core repulsion, i.e. γ → +∞. In the special case of an harmonic potential V (x) = 1 2 mω 2 x 2 , we will show that we recover some known results. But first, let us recall how the exact density can be computed in this case. To keep notations light, we set = m = ω = 1; Exact density in the harmonic trapfinding the eigenstates of a single boson confined in the harmonic trap V (x) = 1 2 x 2 is the same problem as solving the quantum harmonic oscillator. The eigenstates take the form where the n th. eigenstate has energy E n = n + 1 2 and H n is the Hermite polynomial of order n. Now, since bosons with infinite-repulsion map to free fermions [79], the groundstate for N bosons can be built by filling up the first N eigenstates. The density is then given by which is easily evaluated with the Christoffel-Darboux formula When N 1, this can be put in a more explicit form using the asymptotics of the Hermite polynomials, i.e.
Carrying out the asymptotic expansion, we arrive at where the phase θ(x) is the integral of the Fermi momentum k F (x), After some manipulation, we can also write down an asymptotic expression for the density-density correlation, Now, when γ → +∞, the Luttinger parameter is constant, K = 1 (that is the value for free fermion systems). As a consequence, the action (1.8) is conformally invariant, and the Green's functions G D [1] , G N [1/4] as well as the mixed function F D,N [1,1/4] can be obtained explicitly. Indeed, when the strip is conformally mapped to the upper-half plane, it boils down to an exercise in the method of images [80], see e.g. [17]. Below, we will show that we recover exactly the results for the average density and for the density-density correlation. The coordinatex is the stretched coordinate from Eq. (3.8); here, it has an explicit expression, namelyx = π 2 +arcsin x 2N . Finally, with the definition (3.14), the phase ϑ(x) reads Plugging everything into Eq. and its regularized part takes the form for x = (x, τ ) and τ = τ = 0. Plugging these in Eq. (3.20), we recover the celebrated result for the one-particle density matrix in an harmonic trap, see Refs. [17,81,82], where we know that lim γ→∞ B(γ) 2 = G 4 (3/2) In this appendix, we explain how the dimensionful coefficients C (Ô) j can be calculated in practice. For similar discussions that have appeared previously in the literature, see e.g. Refs. [47][48][49][50]59].
We work in the homogeneous, translation-invariant, problem, and the field theory is the usual GFF, which is conformally invariant. Thus, we will rely on conformal transformations and on the operator-state correspondence, namely that the operators φ j in the CFT correspond to eigenstates of the CFT hamiltonian |φ j .
In fact, Eq. (B.1) is strictly valid for an infinite system (x, τ ) ∈ R 2 . Since we will rely on numerical evaluation (i.e. we have finite system sizes L), our first task is to find a way of taking the limit L → ∞. To do so, we start by making the following assumptions: • for sufficiently large system sizes L, the low-energy excited states of the microscopic hamiltonian H can be unambiguously identified with the ones of the CFT Hamiltonian. In particular, the ground state of H for a system of size L, |0 L , is viewed as a microscopic version of the CFT vacuum |0 . Similarly, there is a unique eigenstate of H, noted |φ j L , that is viewed as a microscopic version of the CFT state |φ j .
• the form factor in the microscopic model, L φ j |Ô(x) |0 L , is known for arbitrary L.
With this at hand, the dimensionful coefficient C where ∆ φj is the scaling dimension of the CFT operator φ j . This formula is easily obtained as follows. First, we need to rewrite Eq. (B.1) for a periodic system (x, τ ) ∈ [0, L] × R with periodic boundary conditions in the x-direction. This is done by conformal mapping: the cylinder x + iτ ∈ [0, L] + iR of circumference L is mapped on the infinite plane with the conformal transformation z = e i2π x+iτ L . Then, the r.h.s. in Eq. (B.1) becomes Next, inserting the l.h.s. of Eq. (B.1) in L φ j | . |0 L and the r.h.s. in φ j | . |0 , one gets for τ = 0: where the denominator in the l.h.s. is the normalization of the two microscopic states and the CFT states and operators are normalized such that φ j | φ j (0) |0 = δ j,j in the plane. This is an approximation in finite size L, but it is expected to become exact in the thermodynamic limit L → ∞, hence the formula (B.2). [In the case where the operator φ j has non-zero spin (or equivalently, if the eigenstate |φ j L has non-zero momentum), the dimensionful coefficient possesses an x-dependent phase, which we dropped from Eq. (B.2) for simplicity.] In practice, for the Lieb-Liniger model, we evaluate the coefficients by solving the Bethe equations for a range of particle number N , simultaneously varying the length L = N/ρ such that the density ρ is fixed. We solve the Bethe equations numerically, The eigenstates of the Lieb-Liniger model are indexed by the configurations of Bethe (half-)integers {I 1 , I 2 , . . . , I N }. For instance, it is known that the ground state corresponds to the configura- More generally, any state in the CFT can be identified with a configuration of Bethe roots close the ground state one, with only a few I j 's that are shifted. See e.g. formula (9.18) in the first chapter of the book by Korepin et al. [58] for more information on the relation between the eigenstates of the LL model and those of the free boson CFT.
Given the ground state and an excited state for a given number of particles N , we evaluate the corresponding form factors using the formulae given in Refs. [51][52][53]. We do this for several system sizes N (or lengths L = N/ρ), then perform a polynomial fit in 1/N to get a numerical estimate of the limit N → ∞ in formula (B.2). This is how we obtain the functions A(γ) and B(γ) displayed in Fig. 4.
In Fig. 8 we display the result from our approach for the one-particle density matrix with the non-universal dimensionful coefficient B(x), against the the result where this coefficient is omitted, B = 1. The agreement with the DMRG calculation is much worse in the latter case.

C Electrostatics on the 2d lattice
In this appendix, we study classical electrostatics in an inhomogeneous dielectric medium in 2d, in close connection with the discussion of Sec. 2. This is the construction we use to compute the Green's functions G D [K] , G N [1/4K] and the mixed function F [K,1/4K] numerically; it should therefore help understanding the results of Sec. 2.
Let us start by considering a rectangular lattice whose nodes x are occupied by a discrete height field h, and the Luttinger parameter K lives on the edges, see Fig. 9. Equivalently, one can view this as a resistor network, with the field h viewed as an electrostatic potential V , and with K viewed as a resistor R on each edge. Figure 9: Classical electrostatics in a discretized 2d inhomogeneous medium.
Electric field in an inhomogeneous mediumas we did in the main text, we can look at the electric field E, on an edge xx between two neighboring sites x and x , To keep notations light, we restrict to a square lattice with spacing 1. Using Ohm's law, the current on the (oriented) edge xx is I xx = 1 K E x,x , and the Gauss' law at the vertex x then gives in the absence of an electric charge at site x i . If there is an electric charge α on site x, the r.h.s. is proportional to α; this is the discrete version of the Gauss' law in Eq. (2.12).
In the absence in the absence of a magnetic flux through the plaquettes, the curl of the field E vanishes. For instance, for the sites x i , x j , x o , x k drawn in Fig. 9, which is the discrete version of Faraday's law in Eq. (2.12), Finally, Dirichlet boundary conditions read h x = 0 for x ∈ ∂Ω; in terms of the electric field, this implies that the component tangential to the boundary vanishes on ∂Ω, In electrostatics, this corresponds to the domain Ω being surrounded by a perfect conductor [56].
Magnetic fluxes, electric-magnetic dualitynow, imagine that two plaquettes are pierced by two infinitely thin, constantly increasing, magnetic fluxes in their center, at positions y 1 and y 2 . The fluxes are topological defects around which the field h winds by a constant ±2πβ. In Fig. 9, this is represented by a defect line linking two defects. When h crosses the defect line, it jumps by 2πβ. Notice that we have inserted two defects (the two ends of the defect line) with opposite 'magnetic charge' ±β, in order to be compatible with the Dirichlet boundary conditions. In the absence of electric charges on the lattice sites, the electric field now satisfies where β y is the 'magnetic charge' through each plaquette, here equal to +β if y = y 1 , −β if y = y 2 , and 0 otherwise. On the 2d lattice, the electric-magnetic duality is easily constructed as follows. The dual fieldẼ is defined by a π/2-rotation of E, and a rescaling by 1/(2K), where E 1 and E 2 are the two components of E. This dual field lives on the edges of the lattice, as the original electric field E. But one can viewẼ as the discrete gradient of a dual height fieldh, which lives on the vertices of the dual lattice (i.e. the plaquettes of the original lattice), see Fig. 9. Then, the Gauss' law reads, for the dual fieldẼ, ∇ · 4K∇h = ∇ · 4KẼ = 4πβ y , on a plaquette y with magnetic flux β y . SinceẼ ⊥ ∝ E , it is also clear that the dual fieldẼ satisfies Neumann boundary conditions if E satisfies Dirichlet boundary conditions (E = 0).
x y Figure 10: Configuration used for the definition of the mixed function F x,y on the lattice. F x,y is defined as the electric potential felt on vertex x, knowing that the plaquette y is pierced by a flux: F x,y jumps by ±2π when x crosses the defect line that starts at y (dashed black).
Mixed electric-magnetic potential F x,yfinally, we discuss the mixed function F x,y on the lattice. It is defined as the potential felt by an electric charge at site x in the presence of a single magnetic monopole at site y on the dual lattice. More precisely, we start by fixing y, and a defect line that goes from y to an edge at the boundary, see Fig. 10. The height function h x that lives on the vertices has a 2π-discontinuity along the defect line. This means that the discrete gradient of h along an edge xx , which is usually defined as (∇h) xx = h x − h x , is replaced by (∇h) (d) xx = ±2π + h x − h x on all edges that cross the defect line. The ± sign is fixed by the orientation of the edge with respect to the defect line. One also fixes a function f x that lives on the vertices along the boundary ∂Ω, that has a 2π-discontinuity at the edge where the defect line crosses the boundary, see Fig. 10. The mixed function F x,y is then defined as the height function h x that has the right discontinuity along the defect line, and satisfies the Dirichlet boundary conditions h x = f x along the boundary.
In other words, the mixed function F x,y is defined as the solution of the linear problem where the definition of ∇F is replaced by (∇F ) (d) on edges crossed by the defect line. This is the lattice version of Eq. (2.30) in the main text. So far, we have regarded F x,y as a function of x, defined for some fixed y. But it is interesting to see that it also satisfies a set of dual constraints, as a function of the variable y, ∇ y · K∇ y F x,y = 0, (∇ y F x,y ) ⊥ = (∇f y ) if y ∈ ∂Ω, (C.7) which are the discrete version of Eq. (2.32). We now show that the first equation in (C.7) follows from (C.6); we leave the second one (the boundary condition) as an exercise to the reader. First, we note that, for two neighboring plaquettes y and y , the discrete gradient F x,y − F x,y is the electrostatic potential created by a short defect line on the dual edge yy , viewed at point x. This is illustrated in the following picture, Thus, the combination Q x := ∇ y · K∇ y F x,y corresponds to a sum of four terms, which we can view as the potential created by four short defect lines around y: Now comes the crucial observation: this combination Q x satisfies ∇ x · 1 K ∇ x Q = 0 for all sites x ∈ Ω. For sites that are sufficiently far from the plaquette y, this is obvious, and it simply follows from the fact that F x,y satisfies this equation. However, when x is one of the four corners of the plaquette y, one must be careful with the ±2π discontinuities. Writing the four terms appearing in the explicit expression of the discrete operator ∇ x · 1 K ∇ x , one sees that exactly two of them correspond to terms on edges that cross the defects: The ±2π jumps coming from those two crossings cancel, and the relation ∇ x · 1 K ∇ x Q = 0 holds, as claimed.
In addition, it is clear that Q x = 0 along the boundary x ∈ ∂Ω. Those two facts imply that Q x is identically zero, so ∇ y · K∇ y F x,y = 0 as claimed in (C.7).

D DMRG setup
In this work, Density Matrix Renormalization Group (DMRG) simulation was performed using the open-source C++ library ITensor [54]. The Lieb-Liniger model can be discretized in terms of the XXZ Heisenberg spin chain in its lowdensity regime [83], and in DMRG, this is the most usual way to simulate the LL model (along with the Bose-Hubbard model) [84][85][86][87]. Under this mapping,   Figure 11: Criterion we use to check that the low-density regime is reached. We see that the discretization must be increased as we want to simulation the LL model with smaller interaction parameter γ.

the XXZ Hamiltonian reads
(D.1) where j labels the sites, a 0 is the lattice spacing, n is the total number of sites, J = 2 /ma 2 0 and U = g/a 0 . In the low-density regime a 0 ρ −1 max , the (continuous) position corresponds to ja 0 → x.
We denote by |φ the ground state of Hamiltonian (D.1). The correlation functions of the LL model are then easily computed in terms of the Pauli matrices σ j . The connected part of the density-density correlation is given by ρ(x)ρ(x ) c = φ| σ z j σ z j |φ − φ| σ z j |φ φ| σ z j |φ .

(D.2)
Similarly, the one-particle density matrix can be computed in terms of the raising and lowering operators, In order to check that the low-density regime is correctly fulfilled, we can cook up some criterion. Indeed, performing DMRG in the homogeneous gas, we can extract numerically the prefactors A(γ) and C(γ) from the simulation, and check that they match with the (exact) ones calculated via algebraic Bethe ansatz. In Fig. 11(a), we see that, as γ gets smaller, the two results match as the density gets lower. This seems consistent with the fact that the Bethe ansatz form factors are calculated for ρ → 0. However, since we want to simulate systems with large numbers of particles, we can just as well increase the discretization. Concretely, we set the lattice spacing to a 0 = 1 for n = 512 sites; the, results seem to converge for a 0 decreased by at least one order of magnitude, see Fig. 11(b). In Sec. 3, simulations were performed on a lattice of n = 4096 sites.