Since the early days of Bose-Einstein condensation in ultracold gas experiments [1,2], the momentum distribution of the atoms has been a pivotal experimental observable [3]. Measured via time-of-flight imaging, the momentum distribution has allowed to observe and characterize a wide range of phenomena in a wide range of systems [1][2][3][4][5][6][7]. In recent years, momentum distribution measurements in ultracold gases in onedimensional (1D) and close-to-1D geometries have allowed to observe dynamical fermionization during the expansion in 1D [8], test the accuracy of generalized hydrodynamics [9], study the 2D-1D crossover [10], probe the effect of dipolar interactions in 1D gases [11], observe hydrodynamization after Bragg scattering pulses [12], unveil cooling by dimensional reduction [11,13], and characterize the dynamics of dipolarinteraction stabilized many-body quantum scars [14].
The momentum distribution f (p) is the Fourier transform, f (p) = dx dy e ip(x-y) g 1 (x, y),
of the equal-time correlation function
which is known as the one-body density matrix (OBDM). Here Ψ † (x) and Ψ(x) are the one-particle creation and annihilation operators, respectively, at position x. We set ℏ = 1 throughout * Attila Takács and Yicheng Zhang contributed equally to this work.
our analytical derivations, and reintroduce ℏ when comparing the analytical and numerical results in Sec. V B.
Nonlocal correlation functions like g 1 (x, y) are generally challenging to compute both analytically and numerically. Consequently, predicting the momentum distribution theoretically is difficult, specially in correlated gases out of equilibrium [6,8,9,12,14]. In the context of bosonic gases, this challenge has attracted significant attention from theorists over the years [15]. Analytical results have predominantly been restricted to 1D. Even in 1D, direct calculations for microscopic Hamiltonians are typically limited to hard-core bosons [16][17][18][19], which can be mapped onto noninteracting fermions through the Bose-Fermi mapping [20]. Numerically, the momentum distribution in equilibrium can be obtained using quantum Monte Carlo simulations [21][22][23] or, for integrable gases, by using sophisticated form-factor resummation methods [24,25]. Furthermore, numerical results in and out of equilibrium can be obtained for lattice hard-core bosons at zero [26,27] and finite [28,29] temperatures. Such lattice calculations have been used, in the low-density limit, to understand some of the recent experimental results in the continuum mentioned earlier [8,9,11,12,14].
There also exist well-known asymptotic results for the OBDM of 1D bosonic gases. In particular, in homogeneous (i.e., translational invariant) ground states, the long-distance asymptotic behavior of g 1 (x, y) is predicted by Luttinger liquid theory [30][31][32][33] to be of the general form:
where n is the 1D atom density, k F = πn is the associated Fermi wavevector, and K is the dimensionless Luttinger param-Sec. III, we study trapped 1D gases in equilibrium at finite but low temperature using the "inhomogeneous Luttinger liquid" approach [33,[42][43][44][45][46]. We derive analytical expressions for the OBDM in traps at finite temperature in the hard-core (Tonks-Girardeau) limit and in the weakly-interacting (Gross-Pitaevskii) limit. In Sec. IV, we provide a detailed discussion of our numerical method for evaluating the OBDM in the inhomogeneous Luttinger liquid for arbitrary repulsion strengths, generalizing the method of Ref. [45] to finite temperature. In Sec. IV, we benchmark our approach against exact numerical calculations and compare its predictions to recent experimental results. We summarize our results and discuss potential extensions in Sec. VI.
Throughout this paper, we focus on 1D gases of bosons with repulsive contact interactions. In the absence of an external potential, the corresponding model was introduced and solved by Lieb and Liniger [36]:
where Ψ † (x) and Ψ(x) are bosonic creation and annihilation operators, respectively, at position x in a ring of length L. We set the mass of the bosons m = 1, and c > 0 is the strength of the repulsive contact interaction.
The Lieb-Liniger model can be solved using the Bethe ansatz [36,47,48]. Focusing on the sector with N bosons, the Hamiltonian (5) can be written in the first-quantized form
with associated many-body eigenstates
⃗ x = (x 1 , . . . , x N ), whose explicit expression can be found, e.g., in Refs. [36,47,48]. Importantly, these eigenstates are labeled by a set of real spectral parameters λ = (λ 1 , . . . , λ N ) (or rapidities) whose allowed values are the solutions of the Bethe equations:
or, equivalently, in logarithmic form
Equation ( 9) uniquely specifies a rapidity set λ [hence, an eigenstate of the Hamiltonian (5)] for a set of distinct integers (half-integers) I 1 , . . . , I N for N even (odd). In fact, one may interpret the r.h.s. of Eq. ( 9) as the set of momenta of a noninteracting Fermi gas, and thus one can think of the corresponding rapidities as imposing a nontrivial quantization condition due to the contact interactions. For instance, the ground state set λ GS is obtained from the equally spaced configuration I j = -N +1 2 + j, i.e., by filling a Fermi sea for the associated noninteracting system.
In the thermodynamic limit, N → ∞ and L → ∞ at fixed density n = N/L, it is convenient to replace the rapidity set λ with a smooth density distribution ρ(λ j ) = lim N,L→∞ 1/[L(λ j+1 -λ j )]. In the ground state, the latter satisfies the following integral equation
with λ F fixed by n = λ F -λ F dλ ρ(λ). For later convenience, we also introduce the dressing operation of a generic function of rapidities h, defined as the solution to the integral equation [49] (11) in terms of which ρ(λ) ≡ [1/(2π)] dr . It is also possible to express Eq. ( 10) in terms of dimensionless variables α = c/λ F and g(u) ≡ ρ(λ F u):
normalized such that γ 1 -1 du g(u, α) = α. As a result, one finds that the equilibrium properties of the Lieb-Liniger gas depend uniquely on the dimensionless reduced coupling
with γ → 0 and γ → ∞ corresponding to the limits of weak and strong interactions, respectively. Although the Bethe ansatz approach provides an exact understanding of the spectrum of the Lieb-Liniger model ( 5), the calculation of correlation functions within this framework is a formidable challenge. Since determinant formulas for the field form-factors ⟨λ| Ψ † (0)|µ⟩ between two generic Bethe-ansatz eigenstates have been determined [50], one may express the OBDM (2) as
where λ is the rapidity set of the reference state (e.g., the ground state). The evaluation of Eq. ( 14) requires the summation over the intermediate Bethe-ansatz states |µ⟩, which beyond few-particle systems is a challenging task that needs to be tackled using sophisticated numerical algorithms, see, e.g., Ref. [25]. The numerical evaluation of Eq. ( 14) for the ground state, reported in Ref. [24], yielded the results plotted in Fig. 1.
A. Long-distance asymptotics of the OBDM Alternatively, a universal description of the system's correlations can be obtained at low energies using Luttinger liquid theory [31,32]. The main idea of this effective low-energy theory is to encode low-energy quantum fluctuations in terms of fluctuating bosons on top of the Fermi sea λ ∈ [-λ F , λ F ] obtained through the thermodynamic Bethe ansatz. Following this approach, it is possible to write the so-called harmonic fluid expansion of the field operator [15,[31][32][33]
with the density ∂ x φ and phase θ fluctuating fields satisfying
In Eq. ( 15), B m ≡ B m (γ) are dimensionless nonuniversal amplitudes associated to Umklapp scattering, whose values are obtained from field form factors in the thermodynamic limit, as detailed in Appendix A (see also Refs. [34,35,44,46]). The properties of such fluctuating bosons are determined by the Luttinger Hamiltonian
where v is the sound velocity and K is the Luttinger parameter, respectively. For the Lieb-Liniger model ( 5), these parameters are not independent, v = πn/K, with K = [1 dr (λ F )] 2 ≥ 1 for repulsive interactions. By establishing the two-point correlation of the fluctuating fields via Eq. ( 16), the Luttinger liquid theory enables the calculation of higher-order correlators by means of Wick's theorem. For the specific case of g 1 (x, y), this leads to Eq. ( 3) [15,32,33]. When evaluating the sum in Eq. ( 3), one can use the fact that each harmonic m contributes to the expansion in Eq. ( 15) as a short-distance correction that scales like n -∆m , with ∆ m = 2m 2 K + 1 2K [32]. Hence, by truncating the sum at its leading order (m = 0), one obtains the long-distance asymptotics of the OBDM [valid for |x| ≫ d(n)]
with |x| replaced by L sin(π|x|/L)/π in finite-size systems. Furthermore, the Luttinger liquid theory allows one to account for small thermal fluctuations. By incorporating the effect of a finite temperature T in the low-energy description of fluctuating fields, one obtains:
valid for a translationally invariant gas when ξ T ≪ d(n) ≪ |x|, with thermal length ξ T = πT /v (we set the Boltzmann constant k B = 1). See, e.g., Secs. 3 and 4 of Ref. [32] or the appendices of Ref. [33] for a derivation of Eq. ( 18).
The long-distance asymptotics in Eq. ( 3) exhibits an ultraviolet divergence in the limit |x -y| → 0. This singularity, absent in the microscopic model, is inherent to the Luttinger liquid description and must be regularized in order to obtain the momentum distribution via Fourier transform [cf. Eq. ( 1)]. To this end, the short-distance expansion of g 1 (x, y) for |x -y| ≪ d(n), reported in Eq. ( 4) and derived in Refs. [37][38][39], can be used. Retaining only the lowest orders, one has
with the coefficients expressed in terms of thermodynamic Bethe ansatz quantities
and
.
(22) Here we defined the functions
Strictly speaking, Eq.( 19) is valid for the ground state of the Lieb-Liniger model. However, while finite but low temperatures qualitatively change the long-distance decay of the OBDM from algebraic to exponential [cf. Eqs. ( 17) and ( 18)], they do not significantly affect its short-distance behavior. The validity of Eq. ( 19) for the finite-temperature gas was tested against quantum Monte Carlo simulations in Refs. [40,41].
As discussed in the introduction, Fig. 1 shows the groundstate OBDM for a translationally-invariant gas of size L, obtained by combining the asymptotic results of Eqs. ( 17) and (19). For the long-distance behavior, we retained only the leading-order term (m = 0), although we verified that subleading corrections (m = ±1) do not significantly affect the momentum distribution within the range of momenta that are shown. The matching at intermediate scales |x -y| ∼ d(n) is done by taking the minimum of the two asymptotic curves. While improved results could be obtained through, e.g., a polynomial interpolation between the two asymptotic regimes, we find that our approach provides accurate results without the need of further manipulations.
In experiments with ultracold bosonic gases in 1D geometries, which can be realized using 2D optical lattices [3,15] or atom chips [51], a confining potential V (x) is present, so the corresponding 1D gases are modeled using the Hamiltonian:
The confining potential breaks the Bethe ansatz solvability of the model. Yet, one can use the local density approximation (LDA) to describe local quantities in the inhomogeneous system using the corresponding Bethe ansatz results for the homogeneous gas with a local chemical potential µ -V (x), see, e.g., Refs. [3,15]. This simple approach provides an accurate description of inhomogeneous gases whenever the length scale associated to the changes in the density due to the external potential is much longer than the interparticle distance, namely, whenever
In what follows, we assume that V (x) = V (-x) so that in the ground state the gas is confined in a region x ∈ [-R, R] ⊆ [-L/2, L/2], with R being the "radius" of the atomic cloud. The generalization to non-symmetric traps is straightforward.
In the Bethe ansatz description, LDA is implemented through a position-dependent rapidity cutoff λ F (x) in Eqs. ( 10) and ( 11), fixed such that n LDA (x) = λ F (x) -λ F (x) dλ ρ(λ). Adopting a grand-canonical description for the gas locally allows one to determine λ F (x) from V (x) via the self-consistent equation e dr [±λ F (x)] = 0, with the dressed energy e dr (λ) satisfying
In Eq. ( 25), the chemical potential µ is chosen so that the inhomogeneous gas contains exactly N = dx n LDA (x) particles.
Implementing the LDA in the Bethe ansatz framework results in a position-dependent Fermi surface, on top of which low-energy fluctuations can be incorporated like in the standard Luttinger liquid theory reviewed in Sec. II. This leads to the so-called inhomogeneous Luttinger liquid Hamiltonian [42-44, 46, 52-60]
Notice that the Hamiltonian ( 26) is still quadratic in the fluctuating fields. Thus, exploiting Wick's theorem, it is possible to derive a generic expression for the long-distance asymptotics of the OBDM in the presence of confining potentials Green's functions of the phase fluctuating field θ(x). We also introduced the timescale
associated to the inhomogeneous Luttinger liquid model (26), which is the time needed by an excitation with velocity v LDA (x) to propagate from one edge to the other of the atomic cloud. We stress that ⟨ θ(x) θ(y)⟩ is the phase-phase expectation value computed on the inhomogenous equilibrium state (at either zero or finite temperature) of the trapped Lieb-Liniger gas (24), while G θθ (x, y) = ⟨ θ[s(x)] θ[s(y)]⟩ is the corresponding correlation after the change of coordinate x → s(x) that maps the modulated Fermi surface onto one with unit sound velocity and local Luttinger parameter K[s(x)] (see, e.g., Refs. [42,54,61] and the discussion below). One then needs to define a regularized Green's function for phase-phase correlations occurring at same position x [44,46]
where ultraviolet divergences are removed by exploiting the known result for the homogeneous gas, namely
Equation ( 27) readily provides the long-distance asymptotics of the OBDM in terms of G θθ (x, y). In general, analytical results for G θθ (x, y) are not available, so we treat this function as an input to our theory. In Sec. IV, we discuss an efficient numerical implementation to obtain G θθ (x, y) for arbitrary potentials, and at finite interaction strengths and temperature. Conversely, the limits of strong and weak interactions are exactly solvable and are discussed in the following paragraphs.
Lastly, we note that Eq. ( 27) has singularities in the limit x → y, which parallel those discussed in Sec. II for the homogeneous case. Therefore, the long-distance asymptotics given in Eq. ( 27) must be complemented with a short-distance expansion for |x -y| ≪ min[d LDA (x), d LDA (y)]
with ζ = (x + y)/2 and the coefficients C q (x) ≡ C q (c/n LDA (x)) obtained as simple LDA extensions of the results discussed in Sec. II. Combining together the two asymptotic results of Eqs. ( 27) and ( 30), we approximate the OBDM at all distances to determine the momentum distribution of the trapped gas (24).
A. Finite-temperature weakly interacting bosons in a harmonic trap
In the quasicondensate regime γ → 0 + , namely, at weak repulsive interactions c → 0 and high density n(x) = n(x)/c such that n(x) is finite, analytical results can be derived for the OBDM in a harmonic trap V (x) = 1 2 ω 2 x 2 , see also Ref. [52].
In this regime, our formalism below is equivalent to Bogoliubov theory, see, e.g., Ref. [62].
Starting from the equation of state, µ = cn (see, e.g., Ref. [63]), using the LDA one finds the local density
where R = √ 2µ/ω. Equivalently, R can be related to the number of particles through
In the quasicondensate regime, Bogoliubov theory predicts a divergent Luttinger parameter given by [33] K(x)
while the sound velocity is
In the special case of a harmonic potential, the Luttinger liquid Hamiltonian (26) can be diagonalized using the following mode expansion [52,57,64]
with [â p , â † q ] = δ p,q , [â p , âq ] = 0, and the phonon dispersion in the trap
Importantly, the mode amplitudes a p (u) and b p (u) entering in Eqs. (35) and ( 36) have a known analytical expression in terms of the Legendre polynomials L p (u) [52,57,64]
These functions satisfy the differential equation
and are normalized such that
Using Eqs. ( 35) and (36), the Hamiltonian (26) becomes diagonal in the mode operators and reads (up to an additive constant)
The two-point correlation functions of Luttinger fields computed in the ground state in harmonic traps thus have the following analytical expressions
Notice that the sum entering in Eqs. ( 43)-( 45) is fast converging in p. These formulas can be generalized straightforwardly to finite temperatures by replacing the ground state expectation values âp â † p = 1 and â † p âp = 0 by the thermal equilibrium ones, tr ρth âp â † p = 1 + n BE (ε) and tr ρth â † p âp = n BE (ε), for a the density matrix ρ th ∝ exp -β Ĥinh . The thermal occupation of a bosonic mode is n BE (ε) = 1/(e ε/T -1). One finds that, at finite temperature, Eqs. ( 43)-( 44) become
while ⟨ φ(x) θ(y)⟩ in Eq. ( 45) remains unchanged. Equation ( 46), together with the known result for the nonuniversal amplitude B 0 ≃ 1 for γ → 0 + [44], gives direct access to the asymptotic long-distance behavior of the OBDM (24) in the quasicondensate regime for a harmonic potential. Explicitly,
In Fig. 2, we plot the results obtained evaluating Eq. ( 48) at the center of the trap for different temperatures. The same results can be obtained from Eq. ( 27) by inserting G θθ (x, y) given in Eq. ( 46) (upon using the change of coordinates (50) specified below).
For completeness, we also report the expression of the Green's function for the density fluctuating fields, G ϕϕ (x, y) ≡ ⟨ φ[s(x)] φ[s(y)]⟩, entering, e.g., in the calculation of the density ripples and in the density-density correlations, see for instance Refs. [55,68], or Appendix B for a derivation:
In this section we discuss the numerical method used to determine the equilibrium two-point correlation functions of the Luttinger fields θ(x) and φ(x) in the Hamiltonian (26), valid for arbitrary strengths of the contact interaction at finite temperature. We follow Ref. [45], where a numerical method is described for the zero-temperature Green's functions (see also Ref. [46]), and extend that method to finite temperature.
Our starting point is the Hamiltonian (26), which for convenience we express in terms of the stretched coordinate s(x) in Eq. ( 50), with v LDA (x) = πn LDA (x)/K(x),
Here, we introduced the canonically conjugated momentum Π(s) = ∂ s θ(s)/π such that [ Π(s), φ(s)] = -iδ(s -s ′ ). This Hamiltonian can be readily discretized as follows
where M ≫ 1 is the number of sampling points in the unit interval, and K j ≡ K(s j ) is the discretized Luttinger parameter. The Luttinger fields are replaced with their discretized version Πj , φj satisfying [ Πj , φj ′ ] = -iδ j,j ′ . Open boundary conditions are imposed on the chain, implying that φ0 = φM+1 = 0 and K 0 = K M +1 = 1. In matrix form,
where
and h is the 2M × 2M Hamiltonian matrix having nonvanishing elements
for i, j = 1, . . . , M . It is convenient to change the operator basis from Φ to the bosonic modes
satisfying [ b+ j , bj ′ ] = δ j,j ′ and commuting otherwise, that we collect in the 2M -vector
with 2M × 2M matrix W having nonvanishing elements for j = 1, . . . , M W j,j ≡ W j,j+M = 1/ √ 2,
This notation allows us to recast the Hamiltonian (57) in the quadratic form
which can be diagonalized by a further unitary transformation U . Denoting η = U b, one has
with η † = (η † 1 , . . . , η † M , η1 , . . . , ηM ) and eigenvalues ε j ≡ ε j+M for j = 1, . . . , M , following from the symplectic structure of U required to preserve the canonical commutation relations of b± j operators [45]. Given the structure of the Fock space, it is convenient to consider the associated matrix
where Id is the M × M identity matrix and projects on the negative eigenvalue spectrum of h such that the degeneracy of the spectrum is removed. Denoting as γ j the eigenvectors, and as ω j < 0 the eigenvalues, of h restricted to the negative energy subspace, the two-point correlation matrix can be written as
where P(T ) is a M × M diagonal matrix that projects onto the target state over which the expectation value is computed. In a thermal state, it gives Bose-Einstein weights to the bosonic modes [P T (ω)] i,j = δ i,j 1 -e -ωj /T .
(69) Equation ( 68) fixes the structure of the desired correlation function up to a normalization of the eigenvectors γ j . In order to fix it, we consider the auxiliary (Hermitian) matrix
torus, which reads (see, e.g., Eq. (12.142) in the textbook [65]) ⟨ϕ(z, z)ϕ(0, 0)⟩ torus = -1 2 log ϑ 1 (z|iτ ) ∂ z ϑ 1 (0|iτ ) e -π (Im z) 2 τ .
(B9) Then, applying the method of images (see, e.g., Chapter 9 in Ref. [65] or Refs. [33,43]) to construct the two-point function on an annulus with Dirichlet boundary conditions on both sides, one arrives at Eq. (B8). The same exercise with Neumann boundary conditions leads to Eq. (B7).