The granularity of weakly occupied bosonic fields beyond the local density approximation

We examine ground state correlations for repulsive, quasi one-dimensional bosons in a harmonic trap. In particular, we focus on the few particle limit N=2,3,4,..., where exact numerical solutions of the many particle Schroedinger equation are available employing the Multi-Configuration Time-dependent Hartree method. Our numerical results for the inhomogeneous system are modeled with the analytical solution of the homogeneous problem using the Bethe ansatz and the local density approximation. Tuning the interaction strength from the weakly correlated Gross-Pitaevskii- to the strongly correlated Tonks-Girardeau regime reveals finite particle number effects in the second order correlation function beyond the local density approximation.


Introduction
Observing strongly correlated atomic quantum gases in situ and in real time is quite an achievement [1]- [8]. By controlling the trapping geometry, one can effectively adjust the 'degree' of dimensionality, by feeding in more particles one can approach the thermodynamic limit, a central concept of our macroscopic world, and by controlling the coupling constant one can switch between universal classes of physical systems. Maybe all of this was once envisioned by the great minds who have conceived the very few exactly solvable models of many-body physics [9]- [14], but witnessing the merger of expectation and experiment proves to be an exciting period, today.
In the present paper, we have explored a particular aspect of this rich topic, by focusing on the quantum properties of the ground state containing only a few bosonic particles, i.e. N = 2, 3 and 4, inside a harmonic, one-dimensional (1D) trap. This situation is akin to the atomic or nuclear physics limit of a condensed matter system. There the implied granularity of fermionic matter appears as a shell structure in the energy configuration or on energy surfaces through the appearance of magic quantum numbers. This was unheard of in the field of uncharged gaseous matter until recently with the experimental use of neutral, repulsive, bosonic atoms ( [15] and references therein).
A striking example here is given by the 1D Bose gas. Its closed-form solution has been given for the homogeneous system using Bethe's ansatz in the pioneering work of Lieb and Liniger [11], which focused on the thermodynamic limit (N → ∞ at fixed density). These results have been extended to include finite particle numbers [16] and effects due to a slowly varying trapping potential [17,18], in the sense that the thermodynamic-limit results for the homogeneous system still hold locally even in the presence of an inhomogeneity. However, for small atom numbers N , the trapped system can be solved in a numerically exact fashion, without resort to such a local density approximation (LDA) [19]- [23]. This is the starting point of our paper, which aims at a detailed comparison of the exact correlation functions of N trapped bosons with those obtained in a finite-size homogeneous system under a LDA. That way we map out intrinsic confinement effects and discuss the validity of the LDA for small atom numbers. This is in clear contrast with previous publications, e.g. [22], where two of the authors have investigated, among other quantities, the numerically exact two-body density for N atoms 3 in a 1D harmonic trap. Here the two-body correlation function is normalized, so as to map out trap effects, and then compared with the LDA based on the finite-size Bethe ansatz. This allows us to investigate the transition from the Gross-Pitaevskii (GP) regime to the Tonks-Girardeau (TG) regime for a few particles not only with numerically exact results, but also using analytical expressions, thereby enhancing the insights.
Following this motivation, we will present the basic model of N particles trapped in 1D in section 2. Next, we introduce our physical observables and computational methods in section 3. The homogeneous limit of this system is the Lieb-Liniger model, which represents our benchmark. Its basic notions will be reviewed briefly in section 4 and used in section 5 to discuss our numerical results and analytical modeling of a few harmonically trapped particles.

Model
Let us consider a gas of N 1D bosons with repulsive, short range interactions in a harmonic trap [20], [23]- [25]. Then, the dynamics is given by the dimensionless Hamiltonian where we have measured energy in units ofhω, length in multiples of the harmonic oscillator length a 0 = √h /mω and used the short-hand notation ∂ j = ∂/∂ x j . For example, such an effective 1D description can be obtained by starting with real 3D bosonic atoms of mass m in a strongly anisotropic external trapping potential V (r) = 1 2 mω 2 x 2 + 1 2 mω 2 ⊥ (y 2 + z 2 ). If the transverse level spacing is much larger than in the axial direction β = ω ⊥ /ω 1, we can integrate out two dimensions by assuming that the 2D (yz)-subsystem only occupies the ground state. This procedure leads to the quasi-1D coupling constant g = 2βa s /a 0 , where a s denotes the s-wave scattering length of the bosons.
Nowadays, this situation can be realized experimentally [15] and it is possible to investigate quantum correlations in situ. In particular, we are interested in the properties of the ground state of a few interacting bosons, i.e. N = 2, 3, 4, . . . , and we will explore their quantum correlations. In a homogeneous system of length L = a 0 with a dimensionless scale , the linear number density n = N / is translation invariant and the qualitative behavior of the ground state strongly depends on the correlation parameter γ = g/n. This was first described by Lieb and Liniger [11]. In the thermodynamic limit lim N , →∞ N / = n, it turns out that the state of the homogeneous gas of bosons is completely characterized by this parameter. It is customary to call bosons weakly correlated for γ 1 (the GP regime) and strongly correlated for γ 1 (the TG regime).
The LDA extends this description to weakly inhomogeneous systems under the assumption that the variation of the ground state follows parametrically the spatial variation of the single particle density. In the following, this assertion will be probed by explicitly constructing the N -body wave function with the Bethe ansatz. We will analyze the behavior of the state over the whole range of interaction strengths for an increasing particle number (N = 2, 3 and 4) in order to investigate the transition towards the thermodynamic limit. The correlation functions of the ground state strongly depend on the particle number and this effect is most significant for a few bosons. This analysis will provide a profound understanding of the inhomogeneous system, which can only be solved numerically otherwise.

Observables and computational methods
Assuming we have complete knowledge of the symmetrized and normalized N -particle wave function (x) = (x 1 , . . . , x N ), then we need to extract relevant information about its behavior in terms of experimentally accessible observables [26,27]. Most relevant for this purpose are the number density n(x) = Nρ(x), which is proportional to the single particle density the first-order correlation function measuring phase coherence and the diagonal of the second-order correlation function measuring density fluctuations We will focus on the behavior of the second-order correlation function as it is a more sensitive probe for quantum statistical correlations in the system. Theoretically, much attention has already been directed towards second-order correlation functions [17,25], [28]- [33], but the discussions were mostly concerned with their properties in the thermodynamic limit. In contrast, our interest is directed towards the detailed behavior of the second-order correlation function for few-boson systems, which differs from the thermodynamic limit. The next subsection is devoted to a brief introduction to the multi-configuration timedependent Hartree (MCTDH) method, which is used to calculate the N -body ground state of the Hamiltonian in the presence of a trap in (1).

Multi-configuration time dependent Hartree method
The numerically exact MCTDH method [34,35] is a quantum-dynamics tool that has been applied successfully to systems of few identical bosons [21,22], [36]- [39]. It solves the timedependent N -body Schrödinger equation (ih∂ t − H ) (x, t) = 0 as an initial-value problem by expanding the solution in terms of direct (or Hartree) product states J ( and summing over all admissible configurations C = {J = ( j 1 , . . . , j N )|1 j i s}. In turn, the still unknown, best single-particle functions {φ j (x, t)|1 j s} are represented in a fixed primitive basis implemented on a grid and s denotes the maximum number of required basis functions. The permutation symmetry of (x, t) is ensured by the correct symmetrization of the expansion coefficients a J (t).
Using the Dirac-Frenkel variational principle, one can derive equations of motion for both a J (t) and φ j (x, t) [35]. Integrating this system of differential equations allows us to obtain the time evolution of the system via (5). This has the advantage that the basis { J (t)|J ∈ C} is variationally optimal at each time t. Thus it can be kept relatively small, rendering the procedure very efficient.
Although designed for time-dependent simulations, it is also possible to apply this approach to stationary states. This is done via the so-called relaxation method [40]. The key idea is to propagate an initial wave function (x, t = 0) by the non-unitary, imaginary time propagator U (τ ) = e −H τ . As τ → ∞, any excited state contribution is exponentially suppressed with e −(E m −E 0 )τ and we are left with the ground state. In practice, one relies on a more sophisticated scheme termed improved relaxation [41], which is much more robust, especially for excitations. Here, the energy E = |H | is minimized with respect to both the coefficients a J and the orbitals φ j (x). The effective eigenvalue problems thus obtained are then solved iteratively by first solving for a J with fixed orbitals and then optimizing φ j (x) by propagating them in imaginary time over a short period. That cycle will then be repeated.
Applying this procedure, we obtain the exact wave function for the ground state and can subsequently calculate the correlation functions according to (3) and (4). A profound understanding of the qualitative as well as quantitative details of this correlation function is the main topic of the present study. For this purpose, we compare the results of the MCTDH method to the Lieb-Liniger theory for the homogeneous Bose gas, which will be combined with an LDA. However, before doing so, we briefly review the main concepts of the Bethe ansatz and Lieb-Liniger theory, which can be solved exactly.

Lieb-Liniger theory and Bethe ansatz for the homogeneous system
As we have discussed, it is only possible to obtain the eigenstates and eigenvalues of the trapped Hamiltonian in (1) with a significant computational effort. However, if one can disregard the external trapping potential and supply the system with periodic boundary conditions instead, we obtain the Hamiltonian The corresponding eigenvalue problem has been solved analytically by Lieb and Liniger [11,16] using the Bethe ansatz [9] and they derived the ground state propertieseven for the thermodynamic limit. Before we discuss correlation functions and their functional behavior, we will briefly recall the quintessential steps.
To solve for the ground state and its energy, we have to consider the eigenvalue equation where the spatial region is {x = (x 1 , . . . , x N )| 0 x j } and the totally symmetric wave function satisfies periodic boundary conditions. The dimensionless system length determines the linear number density n = N / for a given particle number N . It still occurs in our dimensionless formulation of the Lieb-Liniger Hamiltonian because we want to allow for a straightforward comparison of the ground states of the Hamiltonians in (1) and (7) for equal particle numbers N and equal interaction strengths g. Furthermore, having the length as a free parameter we can tune the number density n such that the correlation parameter in the homogeneous system is equivalent to the correlation parameter at a certain position in the trapped system.
A subspace of the whole configuration space is the spatial simplex that contains only ascending coordinate N -tuples, i.e.
In this region, (7) together with its periodic boundary conditions are equivalent to  Due to the required symmetry of the wave function under particle exchange, knowledge of (x) in the region R is equivalent to knowing (x) in all other regions of the configuration space. In order to solve (8)-(11), one uses the ansatz where the summation extends over all N ! elements P of the permutation group S N . If we denote the wave vector of the N particles by k = (k 1 , . . . , k N ), then the permuted vector is k P = (k P(1) , . . . , k P(N ) ). For convenience we also introduce the scalar product k P x = N j=1 k P( j) x j . Immediately, one obtains, for the ground state energy of the N -particle system, The remaining task consists of determining the wave vector k such that all boundary conditions are fulfilled. After minor algebra, which is outlined in [11], one obtains the following equations for the components of the wave vectors of the ground state Furthermore, we note that the ground state solution possesses reflection symmetry, which means that for every positive component k j there exists a negative counterpart −k j . Before the ground state wave function can be calculated, we need to clarify how the factors a P are defined, such that the boundary conditions are fulfilled. This is done by first setting a 1 = 1. If P takes k into k P , then this is achieved by subsequent transpositions. For each transposition, the amplitude acquires a factor −e iθ l j , if k l and k j are transposed and k l is to the left of k j . The product of all these factors is a P . Thus, we end up with the following wave functions for N = 2, 3 and 4 particles in the region 0 x 1 · · · x N . In the thermodynamic limit it is possible to switch from a discrete distribution of k values to a continuous distribution and the ground state energy can then be written as [9,11] where e(γ ) is given in terms of the solutions of the Lieb-Liniger equations h(ξ, γ ) = 1 2π Using the Hellmann-Feynman theorem [42], we can directly obtain the second-order correlation function along the diagonal from the ground state energy which can be measured experimentally [5,6], [43]- [46]. In the thermodynamic limit, this simplifies to g (2) (x, x) ≡ g (2) (0, 0) = e (γ ).

Correlation functions of the homogeneous system
To get a feeling for the second-order correlation function and its dependence on the particle number N as well as the correlation strength γ , we first consider the homogeneous system, where this is straightforward. The results are depicted in figure 1 where we compare g (2) (0, 0) in the thermodynamic limit with results that were obtained by using finite values of the particle number N . As can easily be seen, the thermodynamic limit is approached very quickly and for N = 100 bosons the exact results obtained with the Bethe ansatz are almost indistinguishable from the thermodynamic limit.
The decrease of the second-order correlation function for large values of γ is due to an interaction-induced antibunching also known as fermionization of bosons. In the limit γ → ∞, this behavior was first predicted by Girardeau [10] and is due to a one-to-one correspondence between impenetrable bosons and spinless fermions in 1D systems. This can be understood  In figure 2, we plot the results for the wave vectors for an increasing particle number N = 2, 3, 4 and 10. For odd particle numbers, the components of the wave vectors approach even multiples of 2π, which is exactly the result that one would obtain for non-interacting, spinless fermions subject to periodic boundary conditions. However, for even particle numbers, the components of the wave vectors are also separated by a constant spacing of 2π if we go to strong interactions, but they approach odd multiples of π .
This does not correspond to the behavior of non-interacting, spinless fermions subject to periodic boundary conditions. This breakdown of the mapping between impenetrable bosons and spinless fermions has already been pointed out by Girardeau and originates from the periodic boundary conditions. Considering hard wall boundary conditions, then there are no restrictions on the particle number N . Furthermore, it can be noticed that the dependence of the wave vectors on the correlation parameter γ is very similar for the different particle numbers. Hence, the two components of the wave vectors with smallest magnitude for N = 4 and 10 are indistinguishable from the solid line, which depicts the components of the wave vector for N = 2. Apart from the behavior of the wave vectors themselves, a closer look at the wave function reveals that the factors a P appearing in (12) approach either +1 or −1 for γ → ∞. Hence the wave function approaches the limit of a Slater determinant as introduced by Girardeau. A simple way of understanding these behaviors is to consider the defining equations for the wave vectors (14) and (15) in the limit γ → ∞. Using a linear approximation for the arctan, one obtains which exactly yield the previously explained behavior for γ → ∞.
In the remainder of this section, we want to analyze the spatial dependence of the second-order correlation function for the homogeneous problem. Due to the translational symmetry of the system, the correlation functions only depend on the relative distance |x − y|. Hence, it suffices to just study the anti-diagonal of the correlation functions for y = −x. In particular, we will examine the particle numbers N = 2, 3 and 4, because they exhibit finite number effects, which vanish in the thermodynamic limit [23]. In figure 3, we depict the antidiagonal of the second-order correlation function g (2) (x, −x) for N = 2, 3 and 4 particles, respectively. For an increasing value of the correlation parameter γ , we again notice the strong antibunching at the origin, which reduces lim γ →0 g (2) (0, 0, γ ) = 1 − 1/N for ideal bosons to lim γ →∞ g (2) (0, 0, γ ) = 0.
Furthermore, the off-diagonal develops an oscillatory behavior, which becomes most pronounced for large values of the correlation parameter. The oscillatory behavior strongly depends on the particle number. For N strongly interacting particles in the system, we expect N − 1 maxima of g (2) (x, −x) in the interval 0 x /2. This can be understood by recalling the behavior of the wave vectors. In the limit γ → ∞ their components are all equally separated by 2π and therefore exactly N − 1 different combinations of the components exist in the calculation of g (2) (x, −x) when we insert (12) into (4). These combinations range from 2π to (N − 1)2π and lead to trigonometric functions with N − 1 different frequencies, thereby producing the behavior that was described above. Finally, it is also interesting to note that the anti-diagonal of the second-order correlation function has a kink at x = 0, which is due to the delta interaction of the particles. However, in the limit γ → ∞ the kink vanishes and g (2) (x, −x) approaches a smooth behavior at the origin. In general, it is possible to write down the analytic form of any correlation function because we have knowledge of the complete many-particle wave function. Using the particular form of the wave function for N = 2, 3 and 4, as given in (16)-(18), we have calculated the secondorder correlation function according to (4) by properly using the total symmetry of the respective wave functions. However, the lengthy analytical results are not enlightening and we refrain from writing them down in detail. Nevertheless, studying them in the limit γ → ∞, using (23), we find the following simple result where we have assumed that the results for N = 2, 3 and 4 extrapolate to general particle numbers N . This result is consistent with the known limit for non-interacting fermions, as can be seen by evaluating the geometric series This result has already been obtained by Girardeau [10]. However, he first considered the limit γ → ∞ and used the fermionic wave function for the evaluation of the second-order correlation function.

Correlation functions of the inhomogeneous system
With this understanding of the second-order correlation function and its dependence on the particle number N and correlation strength γ , we are now in a position to interpret the trapped results. For N = 2 particles, exact results for the ground state wave function are known [47,48], and for N = 3 and 4 particles we used the MCTDH method. The corresponding results for the spatial correlations and a fixed interaction strength of g = 10 are shown in figure 4.
To establish a common ground for the comparison of the homogeneous and the trapped system, we first of all restrict ourselves to the values of the second-order correlation function in the center of the trap g (2) (x = 0, 0). Tuning the interaction from g = 0.001 → 50, we cover the whole range of the correlation parameter γ from the weakly interacting GP regime to the TG regime of strong interactions. In figure 5, we compare the results in the trapped system for N = 2, 3 and 4 particles (diamonds, squares and circles) with the homogeneous results (solid, dashed and dashed-dotted lines) for the same parameters of g and N . The length of the homogeneous system is chosen such that the constant number density n = N / equals the number density of the trapped system in the center of the trap. In general, we can see good agreement of the results with slight deviations at the crossover from weakly to strongly interacting bosons around γ = 1.
Thus we conclude that the correlation parameter γ remains a valid parameter for the description of the inhomogeneous system as well. It is therefore possible to use a position-dependent correlation parameter γ (x) = g/n(x), which is the central hypothesis of the LDA.
In further considerations, we will apply this position dependent correlation parameter. Thus it is necessary to be acquainted with its spatial behavior. Hence, we plot the number density n(x) for N = 2, 3 and 4 particles for various strengths of the interaction in figure 6. With an increasing interaction strength the density of the ground state changes from a Gaussian shape at weak interactions to a broadened distribution with Friedel oscillations [49] for strong interactions.
From the previous discussion we learn that the second-order correlation function in the center of the trap, g (2) (0, 0), can well be understood by looking at the corresponding correlation function in the homogeneous system for the same parameters of g, N and γ . In the next step we want to investigate the behavior of the diagonal of the second-order correlation function, g (2) (x, x). For this purpose, we take the exact values of the diagonal of the second-order correlation function for the trapped system (obtained with an MCTDH calculation) and plot them as a function of the position-dependent correlation parameter γ (x). Comparing these values with the diagonal behavior of the second-order correlation function in the homogeneous case, presented in figures 1 and 5, basically corresponds to an LDA. As the density decreases by moving out of the center of the trap, increasing values of γ correspond to increasing values of x.
In figure 7, we plot the diagonal of the second-order correlation function g (2) (x, x) for N = 2, 3 and 4 at interaction strengths of g = 0.001, 0.01, 0.1, 1 and 10. Generally speaking, the results for N = 2, 3 and 4 for the trapped system (diamonds, squares and circles) agree rather well with the results of the Bethe ansatz for the homogeneous system in the weakly interacting regime. However, the trapped system deviates significantly from the homogeneous results for large correlation parameters. Starting with values of the correlation parameter around γ ≈ 1, we can see oscillations of the trapped results around the homogeneous curve. This means that the LDA gets less suited to describe the physical content of the trapped system.
The breakdown of the LDA can well be understood if one considers the behavior of the density in the trap. For large interaction strengths the density develops the previously mentioned Friedel-type oscillations which cannot be described within an LDA. This translates to the second-order correlation function where similar oscillations occur.
Whereas for N = 3 the Friedel-type oscillations lead to a peak of the density in the center of the trap, the opposite is true for N = 2 and 4, where we have a dip in the center. In terms of the second-order correlation function, this leads to oscillations around the homogeneous result that either start below the homogeneous curve (N = 3) or above it (N = 2 and 4), as one moves out of the center of the trap. This can be seen in figure 7(b) which magnifies the behavior of the second-order correlation function around γ = 15.
Apart from an oscillatory behavior for large interaction strengths, we have seen that the diagonal of the trapped system can be understood if we combine the Bethe ansatz with an LDA, in the sense described above. Finally, we want to investigate to what extent the full behavior of g (2) (x, y) can be analyzed by the same means. In this last step of the comparison, we calculate the value of g (2) (x, y) in the LDA in the following way. In analogy to the previous discussion, we first of all recall that the comparison of the homogeneous and inhomogeneous systems is made for the same values of the particle number N and the coupling constant g. In the original homogeneous system we have the translational symmetry and the properties of the correlation functions depend only on the relative distance |x − y|. In an inhomogeneous system the relative distance is not the only relevant property that determines the behavior of the second-order correlation function g (2) (x, y). Instead we have to incorporate the spatial dependence of the density to arrive at an LDA. For the diagonal part of the correlation function we have already seen that this combination of the Bethe ansatz and the LDA leads to good agreement with the results for the inhomogeneous system. We can extend this procedure to non-diagonal coordinate pairs (x, y) by choosing the density at the center of the mass coordinate (x + y)/2 for the LDA. As previously, this density is given by the exact MCTDH results and in the next step we again adjust the dimensionless length for the homogeneous system, such that we obtain the same density for the fixed particle number N . Thereafter we solve the Bethe ansatz for these parameters and extract the anti-diagonal value of the second-order correlation function for a relative distance of |x − y|. This result is eventually used to represent g (2) (x, y) for the combination of the Bethe ansatz and the LDA. This procedure is repeated for every coordinate pair (x, y) that is used to plot the second-order correlation function.
In this fashion, we obtain the counterparts of the exact trapped results in figure 4 and depict them in figure 8. The most striking difference that can be noticed by comparing the homogeneous with the trapped results is the large dip in the off-diagonal for N = 2 and 3 in the homogeneous case with LDA. This is merely due to the fact that the periodic boundary conditions in the Bethe ansatz prevent one from going to large distances. For a larger particle number and consequently a larger number density, this difference begins to be negligible in the region of interest, as can be seen in the plot for N = 4. Apart from this artifact the general features are similar and we conclude that we can also understand the overall behavior of the second-order correlation function in terms of an LDA and the Bethe ansatz.
Having a closer look at the anti-diagonal in the center of the trap in figure 9, we get a clearer illustration of these conclusions. While the periodicity prevents one from matching the exact physical behavior for any x, at least the initial slope and the approximate shape for x < 1 are modeled well.

Conclusion
We have examined the ground state correlations for repulsive, quasi 1D bosons in a harmonic trap. In particular, we have focused on the few particle limit N = 2, 3, 4, . . ., where exact numerical solutions of the many particle Schrödinger equation are available with the multiconfiguration Hartree method. These numerical results for the inhomogeneous system are modeled with the analytical solution of the homogeneous problem using the Bethe ansatz and the LDA. Tuning the interaction strength from the weakly correlated GP to the strongly correlated TG regime reveals finite number effects in the second-order correlation function beyond the LDA.
To conclude our contribution, we want to comment on the experimental realization of the particular few-body system that has been discussed and the chances of seeing the deviations of the LDA in the second-order correlation function. The experiments on strongly correlated 1D Bose gases [4,50,51] have dealt with effective atom numbers down to N ∼ 15. However, other experiments have for example measured the dynamics of only two atoms isolated in a superlattice [52]. Moreover, experimental techniques nowadays allow for the control and processing of even single atoms, e.g [53]. Nevertheless, observing effects as in figure 7 requires a high-precision measurement of g (2) (x, x) in combination with the use of quantum-noise detection tools [54], which has not been reported in the papers cited above.