Correlation properties of a one-dimensional repulsive Bose gas at finite temperature

We present a comprehensive study shedding light on how thermal fluctuations affect correlations in a Bose gas with contact repulsive interactions in one spatial dimension. The pair correlation function, the static structure factor, and the one-body density matrix are calculated as a function of the interaction strength and temperature with the exact ab-initio Path Integral Monte Carlo method. We explore all possible gas regimes from weak to strong interactions and from low to high temperatures. We provide a detailed comparison with a number of theories, such as perturbative (Bogoliubov and decoherent classical), effective (Luttinger liquid) and exact (ground-state and thermal Bethe Ansatz) ones. Our Monte Carlo results exhibit an excellent agreement with the tractable limits and provide a fundamental benchmark for future observations which can be achieved in atomic gases, cavity quantum-electrodynamic and superconducting-circuit platforms.


I. INTRODUCTION
Understanding the role of strong correlations in complex quantum many-body systems such as gases of interacting atoms or electrons is one of the most important challenges in modern condensed matter physics, material research, and chemistry. In this context, one-dimensional (1D) Bose gases composed by particles with contact repulsive interactions [1] have proven to be a versatile testbed for the study of strong correlations in quantum many-body physics, where solidstate, low-temperature, and atomic physics converge. They offer a simple yet nontrivial model which can still be captured with reasonable theoretical effort [2], by allowing for the direct comparison between exact and analytical results.
The 1D Bose gas model may be precisely simulated by cavity quantum-electrodynamic devices, by shedding light on the behavior of its correlation properties [3]. Most importantly, it has been experimentally realized in superconducting circuits where correlations have been indeed measured from weakly to strongly interacting regime [4]. However, the best experimental platform so far was provided by ultracold atoms as they offered an exquisite and precise control over many system parameters [5], including the interaction strength between atoms. These ultracold atomic gases are so dilute that the indepth understanding of the physics at stake is possible in most cases. Finally, they can be realized in different spatial dimensions, including 1D geometry.
At finite temperature, many momentum modes in the longitudinal 1D direction are excited. This leads to markedly different behavior of a 1D Bose gas compared to the threedimensional Bose-Einstein condensate (BEC), where the lowest momentum mode remains macroscopically occupied at low temperature. The excitation of momentum modes with temperature is the origin of highly enhanced density and phase fluctuations which prevent the off-diagonal long-range order [6,7] and then the formation of a true BEC. Instead, an extraordinarily rich landscape of many different degenerate regimes emerges [8][9][10][11][12][13][14], which are separated by smooth crossovers and they might or not share the typical features of a BEC. They can be explored by changing the interaction strength and temperature.
Distinct quantum regimes can be conveniently classified by the characteristic behavior of thermodynamic properties. Asymptotic analytical limits have been described by a variety of physical models [13,14] and compared with the exact solution provided by the thermal Bethe-Ansatz (TBA) method within the Yang-Yang theory [15,16]. A new intriguing quantum regime in a 1D Bose gas has been recently predicted by the same authors of the present manuscript [14]. It is signalled by a thermal anomaly in the temperature dependence of the specific heat for any finite value of the interaction strength. Calculations of the dynamic structure factor with the ab-initio Path Integral Monte Carlo (PIMC) method showed that, at temperatures similar to the anomaly threshold, the excitation pattern experiences the breakdown of the quasiparticle description holding instead at low temperature [14]. Both behaviors of the specific heat and the dynamic structure factor suggest that the novel so-called hole anomaly, which is induced by the presence of unpopulated states in the excitation spectrum of the system, is a reminiscence of the superfluidnormal phase transition at the critical temperature, although any phase transition is not allowed in 1D geometry [17].
Different quantum regimes in a 1D Bose gas can be also identified with the corresponding analytical limits of the pair correlation function (also called two-body distribution function) [8][9][10]. However, the full function cannot be readily obtained with TBA and require instead ab-initio techniques. At zero temperature, several correlation properties have been calculated for arbitrary interaction strength with quantum Monte Carlo [18,19] and density matrix renormalization group (DMRG) [20] techniques. Previous finite-temperature studies only focus on very limited regimes of interaction and temperature and are based on different approaches. In the weaklyinteracting limit, some mean-field methods, including Bogoliubov theory [21,22], and the stochastic Gross-Pitaevskii equation [23] can be employed at low temperature, while the semiclassical expansion [24] is valid only at high temperature.
In the strongly-interacting regime, available methods are the Bose-Fermi mapping [25] and the numerical stochastic gauge technique whose results are only reported for a single value of temperature [26]. At temperature below the hole anomaly, correlation properties have been computed with Bethe-Ansatz method [27], PIMC technique using worm updates [28], by relying on the mapping with the sinh-Gordon model [29,30] and on the form factor approach [31]. At high temperature, imaginary time stochastic gauge-P simulations can be employed [10]. Low-momentum [32] and large-distance [33][34][35][36] correlations have been also explored.
The present work fills an important existing gap by providing a complete study of the stationary microscopic correlation properties of a 1D Bose gas for a wide range of interaction strength, all over from the weakly-to the strongly-repulsive regime, and for temperatures crossing from the quantum to the classical gas. We calculate the spatial dependence of the one-body density matrix and the pair correlation function, as well as the static structure factor, as a function of the momentum in a very broad range of values. This calculation is carried out using the PIMC technique which allows us to obtain unbiased results with a controllable accuracy for any temperature and interaction strength values. Our results show an excellent agreement with exact thermal Bethe-Ansatz solution [15,16] and analytical limits in regimes of their validity. Our PIMC findings are fundamental for benchmarking the theoretical predictions and to stimulate future experimental measurements.
The structure of the paper is as follows. In Sec. II we introduce the studied model of the 1D Bose gas at finite temperature. Section III is devoted to a brief description of the Path Integral Monte Carlo method used for the numerical results of the correlation functions. Our findings of the pair correlation function, the static structure factor, and the one-body density matrix are reported in Secs. IV, V, and VI, respectively, and they exhibit an excellent agreement with analytical and thermal Bethe-Ansatz limits. In Sec. VII, we discuss possible upcoming experimental observations of our predictions. Finally, in Sec. VIII, we draw the conclusions and future perspectives of our work.

II. MODEL
The many-body Hamiltonian of a 1D homogeneous gas of N Bose particles interacting via repulsive contact pseudopotential (Lieb-Liniger model) is given by where m is the atomic mass, g = −2 2 /(ma) > 0 is the positive 1D coupling constant [37], and a < 0 is the negative 1D s-wave scattering length. The dimensionless interaction strength γ = −2/(na) is related to the gas parameter na, where n = N/L is the linear density with L being the length of the system. We are interested here in the thermodynamic properties which are obtained by taking the N, L → ∞ limit while keeping the density n constant. There is a continuous crossover which encompasses different quantum degeneracy regimes. In the Gross-Pitaevskii (GP) limit of weak repulsion γ 1 and high density n|a| 1 the gas admits a mean-field description [38]. In the Tonks-Girardeau (TG) regime [39,40] of very strong repulsion γ 1 and low density n|a| 1, bosons become impenetrable and they cannot cross each other. This constraint, together with the spatial peculiarity of 1D systems, acts as an effective Pauli exclusion principle, resulting in a dramatic suppression of three-body losses and making the TG gas stable. In this highly correlated limit, the wavefunction of strongly repulsive bosons can be mapped onto that of an ideal Fermi gas, resulting in identical thermodynamics and excitation spectrum [40]. Seminal experiments have explored the GP-TG crossover in the past years [41][42][43][44][45][46][47].
At zero temperature, the Bethe-Ansatz method allows one to calculate exact energetic properties of the Lieb-Liniger model such as the ground-state energy E 0 , chemical potential µ 0 = (∂E 0 /∂N ) a,L and sound velocity v = n/m(∂µ 0 /∂n) a which are all functions of the interaction strength γ [38,[48][49][50]. The speed of sound v monotonically increases with γ from the mean-field v GP = gn/m to the Fermi v F = πn/m value in the TG limit.
At finite temperature, the complete thermodynamics within the canonical ensemble can be inferred from the Helmholtz free energy A = E − T S, where E is the internal energy and S = − (∂A/∂T ) a,N,L the entropy. Its knowledge can be used to calculate the pressure and the inverse isothermal compressibility where µ = (∂A/∂N ) T,a,L is the chemical potential at finite temperature and we have employed the Gibbs-Duhem relation dP = ndµ + sdT with s = S/L the entropy density. The quantity (∂n/∂µ) is directly related to the atom number fluctuations which can be measured in ultracold atom experiments. In this way, one has access to the isothermal compressibility, κ T , equation of state and temperature T [51][52][53][54][55]. The experimental estimate of temperature is also possible through the use of Eq. (3) and the combination of the measurements of the density n, the isothermal compressibility κ T and the pressure P [56,57].

III. PATH INTEGRAL MONTE CARLO TECHNIQUE
The PIMC is a powerful stochastic method which allows to obtain numerically structural and energetic properties of a quantum system at finite temperature [66] which is described here by the microscopic Hamiltonian (1). Some information on the structure of the excitation spectrum can be inferred as well from the dynamic structure factor by using the inverse Laplace transform [14]. PIMC method relies on the Feynman path-integral description of an ensemble of N quantum atoms in terms of a set of N classical polymers, each of them reproducing a quantum delocalized particle [66]. In this way, the thermal average or expectation value O of the quantum observable O is expressed as a multidimensional integral which can be efficiently computed via Monte-Carlo sampling of the coordinates R: where we have introduced the thermal density matrix n T = e −βH /Z, the partition function Z = Tr e −βH , the inverse temperature β = (k B T ) −1 and H is the Hamiltonian, Eq. (1). In the following, we work in coordinate representa- . , x N,i } denotes a set of the coordinates of the N atoms of the system and G(R 1 , R 2 ; β) is the Green function propagator describing the evolution in the imaginary time β from the initial R 1 to the final R 2 configuration. The configuration PR in Eq. (5) is obtained by applying a permutation P of the particle labels to the initial configuration R and the sum P over the N ! permutations allows to take into account the quantum statistics of the identical bosonic atoms.
The key aspect of the Path Integral formalism is the convolution property of the propagator which can be generalized to a series of intermediate steps R 2 . . . R M (technique known as trotterization) defining a path with M configurations and total time β = εM , where ε is the time step. For a finite value of M , the path is discrete in time. In the limit of large M , the path becomes continuous and ε approaches zero, so each step corresponds to high temperatures T . In this classical gas limit, the propagator admits an analytical approximation where the quantum effects of the non-commutativity between the kinetic and interaction potential operators in the Hamiltonian are neglected. The thermal expectation value, Eq. (5), can be then approximated as where the Bose-Einstein statistics imposes the boundary condition R M +1 = PR 1 . The probability distribution is positive definite and its integral over the space of configurations is equal to the unity. The PIMC approach is based on a stochastic evaluation of the integral in Eq. (7) by sampling N × M degrees of freedom according to the probability distribution p(R 1 , . . . , R M ). In order to improve the permutation sampling, we use the worm algorithm [67]. Decomposition (7) becomes exact in the limit M → ∞ where the imaginary time ε is small and the analytical high-temperature approximation for the propagator G is accurate, allowing for an exact calculation of the thermal averages O with the PIMC method. We approximate the propagator within a pair-product scheme which is based on the exact solution of the two-body problem for the Hamiltonian in Eq. (1) [68,69]. The optimization of the number of the convolution terms M in Eq. (7) has been achieved by benchmarking the PIMC expectation values of the internal energy per particle E/N and the isothermal compressibility κ T as a function of temperature and interaction strength against the exact thermal Bethe-Ansatz results [14].
The PIMC method requires an average time of several days for performing a single calculation as the computational cost scales quadratically with the number of particles N . This scaling is much more favorable as compared to the diagonalization methods which typically have an exponential cost. PIMC thus allows us to efficiently calculate various observables, including correlation functions, by obtaining a mean value and an error bar. The errors are characterized by two different contributions. The first one is the statistical error which can be reduced by increasing the number of iterations. It typically scales with the inverse of the square root of the number of iterations. The second one is the systematic error which is controlled by decreasing the time step ε and increasing the number of configurations M . By fixing ε, it is easy to see that M is proportional to the inverse of temperature M ∼ β ∼ 1/T . The lower is the temperature, the larger is the required M to keep the systematic error sufficiently small. In the classical limit of high temperatures, the effects of quantum delocalization of particles are less relevant and a small value of M is sufficient for an accurate description of the many-body system. In the quantum regime of low temperatures, it is necessary instead a large M which makes the calculation more time demanding. By fixing the temperature, the computational cost also depends on the interaction strength γ which scales with the inverse of the linear density γ ∼ 1/n where n ∼ N . In the Gross-Pitaevskii regime of weak interactions γ 1 and high density, a large number of atoms N is needed, making calculations challenging. By approaching the opposite Tonks-Girardeau regime of strong interactions γ → ∞, N can be relatively small. However, the importance of quantum correlations is enhanced by increasing γ, thus requiring a large value of M . To summarize, within the PIMC calculations, it is possible to keep the error bars under control and to decrease them in a systematic way at the expenses of longer simulation times.
In the following Sections, we report the PIMC results for the pair correlation function, the static structure factor, and the one-body density matrix, whose expectation values have been computed with Eq. (7) for a broad range of values of interaction strength γ and temperature.

IV. PAIR CORRELATION FUNCTION
The pair correlation function or normalized density-density correlator quantifies the probability of finding two particles separated by a distance x [38], according to whereψ (x) is the bosonic field operator and · · · denotes an average over an ensemble at thermal equilibrium at a given temperature T . Here n = ψ †ψ is the diagonal density where the operators are evaluated at the same spatial coordinate. The pair correlation function provides the characteristic length scale over which the density-density fluctuations decay. In a 1D Bose gas, g 2 (x) plays a key role in the study of interference properties of atom lasers in a 1D waveguide [70], by fostering practical applications such as matter-wave interferometry [71]. In the experimental image of an expanding gas cloud, g 2 (x) can be employed to probe complex many-body states of trapped ultracold atoms [72]. The local pair correlation function g 2 (0), at x = 0, gives the probability of two particles to overlap. Its value is related to the derivative of the free energy with respect to the coupling constant (∂A/∂g) T,N,L [8] according to the Hellmann-Feynman theorem [73]. By using Eq. (4) and the definition of the interaction strength γ, one finds the relation with the Tan's contact C [38]: Its value can be calculated exactly in a 1D Bose gas through the thermal Bethe Ansatz [8,13,15,16]. Since g 2 (0) and C are proportional, Eq. (9), they both connect the short-range behavior of correlations with the thermodynamic properties. Analytical limits of g 2 (0) and C in various physical regimes are reported in Refs. [8,74] and [13], respectively. At large relative distances in a gas phase, atoms become uncorrelated and the pair correlation function approaches the g 2 → 1 value.
In the decoherent classical (DC) regime, an excellent approximation for the pair correlation function is [9,10] where σ = λ/ √ 2π is proportional to the thermal de Broglie wavelength λ = 2π 2 / (mk B T ) and erfc (x) is the complementary error function. The DC regime holds for k B T max 2 2 /(ma 2 ), 2 n 2 /(2m) . That is, the temperature must be large compared both to the characteristic energy associated with the s-wave scattering and to the quantum degeneracy temperature T d = T F /π 2 related to the Fermi energy E F = k B T F = 2 π 2 n 2 / (2m). For T T d , the chemical potential µ is large and negative, so the bosonic occupation number is small n(k) 1 and can be approximated by the Maxwell-Boltzmann distribution n(k) MB ≈ e −β[ 2 k 2 /(2m)−µ] , where k is the momentum. Such an approximation yields Eq. (10). In the DC regime, both phase and density fluctuations are large, thermal effects always dominate over the interactions and the system approaches the ideal Bose gas behavior at very high temperature even if quantum effects still play a role [8][9][10]14], as signalled by the corrective term in Eq. (10) depending on the interaction strength γ. The parameter σ is thus always smaller than the absolute value of the 1D s−wave scattering length: σ < |a|.
In the non-interacting limit γ → 0 of Eq. (10), the classical ideal gas result [75] is recovered and characterized by the Gaussian Maxwell-Boltzmann decay at large distance with a correlation length fixed by σ which determines the approaching to the uncorrelated limit g 2 → 1.
In Figs. 1-3, PIMC results for the pair correlation function g 2 (x) are reported with symbols for several values of the interaction strength γ and temperature in units of the quantum degeneracy value T d = 2 n 2 / (2mk B ). Horizontal lines denote the results of exact TBA values for the local pair correlation function g 2 (0), Eq. (9), which is a monotonic increasing function with temperature at fixed γ. PIMC results are in a perfect agreement with the TBA predictions for the local correlations at x = 0 and provide a full description for non-local Below the quantum degeneracy temperature T T d and for any interaction strength γ, the amplitude of the longrange decay of the pair correlation function g 2 (x) is proportional to 1/x 2 and corresponds to the Fourier transform of the linear phonon expression for the static structure factor, S(k) = |k|/(2mv). Another main feature observed in g 2 (x) is the antibunching effect that is voiding the region around x = 0 as the repulsive atoms tend to avoid each other. In the limit of strong repulsion (large γ), pair correlation function approaches that of an ideal Fermi gas (IFG), g 2 (x) IFG ≈ 1 − sin 2 (πnx)/(πnx) 2 , for which the voiding is complete, g 2 (0) IFG = 0 [9,10,19,25,76]. Here, g 2 (x) IFG features Friedel oscillations manifested as a series of maxima located at the multiples of the mean interparticle distance ∼ 1/n. Friedel oscillations can be interpreted as a interference effect between incident and reflected wave in the twobody collisions on an impenetrable potential. Friedel oscillations are also found in the density profile of a 1D interacting electron gas with an impurity [77]. Since local maxima in the pair correlation function imply the existence of more likely separations between atoms, Friedel oscillations can be also understood as a quasicrystalline order (with a period ∼ 1/n) in the two-particle sector of the many-body wavefunction even though the density of the gas is uniform. The ideal Fermi gas limit g 2 (x) IFG is equally valid for strong attractive interactions, i.e., for large and negative γ, by describing the pair correlations of a metastable state known as super Tonks-Girardeau gas [78][79][80]. For weaker interactions (small and intermediate values of γ), g 2 (x) is a monotonically increasing function at very low temperatures.
Thermal fluctuations lead to a qualitatively different behavior due to the bunching effect in which the probability of two particles to overlap is enhanced, so that the local value for g 2 (x) increases with temperature until the maximum value g 2 (0) = 2 of the classical gas is reached. As a result of the competition between the thermal bunching and antibunching due to the interparticle repulsion, the pair correlation function shows a global maximum larger than unity g 2 (x max ) > 1 located at a finite interparticle distance x max > 0. This global maximum is formed above a threshold temperature which is higher for larger values of γ. The value in the maximum, g 2 (x max ), increases while its position x max decreases with temperature. At high temperatures, all pair correlation functions, calculated at different values of γ, approach the Gaussian shape as predicted by the decoherent classical theory and given by Eq. (10). We find an excellent agreement between the PIMC results and the DC theory. In this regime, the local pair correlation function is always larger than unity, g 2 (0) > 1. The maximum value is achieved for local correlations x → 0, where strongest bunching effect is observed. The limiting value g 2 (0) → 2 is the same as in an ideal Bose gas and is fully driven by the thermal energy. The DC regime is achieved at higher temperatures for stronger interactions [14] and the value of the global maximum in the pair correlation function is approximated by [9,10] where τ = T /T d . This maximum is more clearly visible by increasing γ, Figs. 1-3.
In Fig. 1, we show the characteristic examples of the pair correlation function in the weakly-interacting Gross-Pitaevskii regime with γ = 10 −1 . Each panel corresponds to different thermal regimes in order of increasing temperature values: hole anomaly (upper), intermediate (middle) and DC regime (lower) [14]. In the limit of zero temperature, g 2 (0) ≈ 0.7 for this value of γ, that there is a weak suppression of the probability of two particles to overlap, as compared to the uncorrelated (ideal Bose gas) value, g 2 (0) = 1. At very low temperature, there is a quasicondensate in the system, which is characterized by suppressed density fluctuations due to repulsive interactions [81] and enhanced longwavelength phase fluctuations due to thermal excitations [82]. The resulting local pair correlation function indicates secondorder coherence g 2 (0) ≈ 1, for which there exists a finite probability that two particles come close to each other. So that, the existence of the quasicondensate is fully driven by the correlations induced by very weak interatomic repulsion, similarly to the case of the DC regime g 2 (0) ≈ 2 where interaction effects can be neglected.
In Fig. 2, we report the pair correlation function for the intermediate value of the interaction strength γ = 1 in the hole anomaly (upper panel) and in the DC (lowest panel) regimes [14]. This is the most difficult regime to describe as it cannot be treated perturbatively and one has to rely on numerics in order to obtain accurate results. At zero temperature, g 2 (0) ≈ 0.5, in agreement with thermal Bethe-Ansatz calculations. The potential energy of two-body interactions, E pot /N = g 2 (0)gn/2, is maximal in this regime, as compared to the Gross-Pitaevskii regime with vanishing coupling constant, g → 0, and the Tonks-Girardeau regime of impenetrable particles, g 2 (0) = 0. As the temperature is increased, antibunching effect becomes less prominent, until temperature T ≈ 4T d is reached with g 2 (0) ≈ 1. Further increase in the temperature leads to bunching effect.
In Fig. 3, we show the pair correlation function in the case of strong repulsion, γ = 10. Different thermal regimes are reported in each panel: hole anomaly (upper), virial hardcore (HC) and beyond (middle) and DC regime (lower) [14]. At zero temperature, the local pair correlation function approaches zero g 2 (0) ≈ 0.1, by signalling the antibunching effect. The gas is then fermionized as the strong repulsion mimics the Pauli exclusion principle for intrinsic fermions at zero relative distance and pairs of bosons are thus never at the same spatial position, g 2 (0) = 0. The antibunching is washed out by thermal effects. At the lowest temperature here reported T = 2T d , Friedel oscillations are not visible at larger distances, by confirming that they exist only below the quantum degeneracy point T T d . In the virial hard-core regime (we refer to the middle panel of Fig. 3) of strong interactions and high temperature π 2 < τ γ 2 [13,14] the local pair correlation function does not manifest complete antibunching g 2 (0) = 0, as witnessed by exact TBA and PIMC results. Instead, its value remains finite (0.4 g 2 (0) 0.6) which is in contrast to what is claimed in Refs. [9,10].
The crossover from the virial HC to the DC regime can be induced by increasing the temperature [14]. It is signalled by the global maximum in the pair correlation function emerging at temperature τ ≈ γ 2 corresponding to σ ∼ |a| [10] where a is the 1D scattering length and σ, up to some multiplicative constant, corresponds to the thermal de Broglie wavelength. Such a crossover is explained by the interplay between the thermal bunching and repulsion antibunching effects acting on comparable scales σ ∼ |a| [10]. In the crossover regime γ 2 /τ 0.1 − 0.4 (middle and lowest panels of Fig. 3), the  local pair correlation function exhibits a behavior similar to the quasicondensate g 2 (0) ≈ 1, indicating local second-order coherence. However, unlike the quasicondensate regime, the non-local correlations on distance scale |x| ∼ σ ∼ |a| are not coherent but bunched as witnessed by the presence of the maximum [10].
In the high temperature regime (lowest panel of Fig. 3), we find an excellent agreement between the PIMC results and the analytical limit provided by the decoherent classical theory, Eq. (10), except for the lowest value of temperature for which the deep DC regime is not yet achieved [8][9][10]14].
It is important to point out that other known analytical limits for the pair correlation function g 2 (x) (weakly-and stronglyinteracting at high temperatures and decoherent quantum regimes) [9,10] are not reported in the present work as they do not show agreement with our exact PIMC results. Comparison with the approximate theories of the weakly-and stronglyinteracting quantum regimes [9,10] is not possible due to the growing computational cost of PIMC method going down to extremely low temperatures. While our exact PIMC and TBA results for the local pair correlation function g 2 (0) agree for any considered value of the interaction strength and temperature, they do not match with the corresponding analytical limits [8] in some regimes. Given the proportionality with Tan's contact parameter, Eq. (9), g 2 (0) is an additional thermodynamic property. New analytical descriptions have been built and they exhibit an excellent agreement with TBA in their regimes of validity and for several thermodynamic quantities [13,14,50].

V. STATIC STRUCTURE FACTOR
The static structure factor S (k) is directly related to the Fourier transform of the pair correlation function, Eq. (8) [38]: and quantifies the spatial correlations in momentum space, where k = k 1 − k 2 is the relative momentum of atoms. At small k, the static structure factor is sensitive to thermal and dynamical correlations and at k = 0 is related to the isothermal compressibility κ T , Eq. (3) [38] S (0) = k B T κ T .
Its value can be computed exactly through the thermal Bethe-Ansatz calculation of κ T [14][15][16]27] which has been recently compared with corresponding PIMC results, by showing an excellent agreement [14]. Equation (12) is a general result following from the fluctuation-dissipation theorem and holds for any value of temperature and interaction strength.
For large momenta, k → ∞, the static structure factor approaches the uncorrelated value, S (|k| n) → 1 [38]. If the excitation spectrum (k) is exhausted by a coherent single-mode (SM) quasiparticle [14] then the static structure factor can be approximated as follows [38] The prefactor in Eq. (13) recovers the Feynman relation [83] and corresponds to the static structure factor at zero temper- In a previous publication [14], PIMC results for S (k) were used to estimate the spectrum within the Feynman approximation and to show that it is valid for temperatures below the hole anomaly value. The low-temperature correction in Eq. (13) is derived from the principle of the detailed balance. By using the T = 0 Bogoliubov (BG) spectrum (k) = ( kv) 2 + ( 2 k 2 /2m) 2 [48,49] in Eq. (13), one can analytically obtain the static structure factor at low temperature in the weakly-interacting regime τ √ γ 1 [8-10, 13, 14]. For large momenta, one can use the expansion of the BG spectrum ( |k| mv) ∼ 2 k 2 /(2m) + mv 2 in Eq. (13) and obtain at the level of the Feynman approximation [38]: which recovers the model-independent limit S (|k| n) → 1 after a further expansion |k| mv. For arbitrary values of the interaction strength γ, the first term in Eq. (14) is recovered by assuming the free-particle excitation spectrum, (k) = 2 k 2 /(2m) in Eq. (13) within the Feynman approximation. Equations (13) and (14) are very general as they are also valid for a Bose gas in three spatial dimensions [38].
At small momenta, |k| mv, the excitation spectrum is dominated by linear phonons (k) = v |k| propagating with sound velocity v. This allows to obtain from Eq. (13) a generic expression for the static structure factor at small momenta, The speed of sound itself is related to the zero-temperature compressibility, and can be obtained from the equation of state, µ 0 (n), which is expressed in terms of the chemical potential at zero temperature µ 0 . Equation (15) is then a reliable small-k approximation for the static structure factor which corresponds to the Luttinger Liquid (LL) regime holding at very low temperatures k B T µ 0 and for any interaction strength γ [50,84]. Eq. (15) recovers the classical result, Eq. (12), where the isothermal compressibility is calculated at zero temperature In Figs. 4-6, we report the static structure factor S (k) calculated for characteristic values of the interaction strength γ and temperature in units of the quantum degeneracy value T d = 2 n 2 / (2mk B ). In Fig. 4, we report S(k) in the weaklyinteracting GP regime with γ = 10 −1 . Each panel corresponds to a different physical thermal regime ordered by increasing temperature values: hole anomaly (upper) and intermediate (lower) [14]. In the upper panel, black dashed lines denote Eq. (13) calculated with the Bogoliubov spectrum and black solid line represents the corresponding large-k behavior, Eq. (14). We show S(k) calculated in the regime of the intermediate interactions, γ = 1, in Fig. 5. The case of strong correlations, γ = 10, is considered in Fig. 6.
Horizontal lines denote the exact zero-momentum static structure factor S (0) as calculated by TBA, Eq. (12). In the limit of zero temperature, zero value is recovered, S (0) = 0, for any γ. S (0) quantifies the interplay between the quantum and thermal fluctuations. The effect of the temperature is more prominent in weakly-interacting systems. For example, at the degeneracy temperature T = T d , the thermal effects are large for weak interactions, γ = 10 −1 , and moderate for strong interactions, γ = 10. In fact, in the weakly-repulsive  regime, the inverse isothermal compressibility κ −1 T , Eq. (3), is small and S (0) grows rapidly as the temperature is increased. In the opposite fermionized regime with γ 1, κ −1 T is large and S (0) remains small even at relatively high temperature. A characteristic thermal crossover temperature, defined by the condition S(0) = 1, is an increasing function of γ and reaches values as large as τ = T /T d ∼ 10 2 (not reported in Fig. 6). Above this temperature threshold, the static structure factor develops a global maximum S(k max ) = S(0) > 1 which defines the most likely momentum k max of atoms. The zero-momentum value, S (0), is a non-monotonic function of temperature and it exhibits a maximum occurring at a certain temperature τ max . The temperature τ max itself increases monotonically with increasing the strength of interactions γ. This temperature is higher than that of the hole anomaly and smaller than that of the DC regime [14] for any interaction strength. Such universal result holds even for the largest con-   sidered interaction strength, γ = 10, where the maximum of S (0) occurs at temperature around τ max ∼ 10 4 not shown in Fig. 6. It is instructive to identify the temperatures and momenta for which the Luttinger Liquid description, Eq. (15), is applicable. This theory is based on assuming a linear phonon dispersion relation and has the same underlying assumption as the single-mode result given by Eq. (13). By making a comparison of this effective Luttinger Liquid theory (coloured dashed lines in Figs. [4][5][6] with the exact numerical PIMC result, one finds that the single-mode phonon assumption holds well at low momenta and low temperatures. For a fixed value of γ, such low-k approximation holds for a broader range of k at lower temperatures. The range of applicability gets broader by increasing γ provided the Luttinger-Liquid condi-tion |k| mv being v a monotonic increasing function with γ [49,50]. Our results also confirm the temperature-range of validity of the Luttinger-Liquid description k B T µ 0 . Such a range gets broader by increasing the interaction strength as the zero-temperature chemical potential µ 0 is an increasing function with γ [48,50].
For each set of results at fixed γ, the decoherent classical regime at high temperature is not reported as the static structure factor is trivially a constant function for any value of k and approaches the model-independent value S (k) ≈ 1.

VI. ONE-BODY DENSITY MATRIX
The one-body density matrix (OBDM) quantifies the coherence and is defined as [38] It is proportional to the amplitude of the process in which one particle is annihilated at position x 2 and the same state is recovered by creating a particle at position x 1 . Definition (17) applies to any value of the interaction strength and temperature. For a zero displacement, x = 0, one recovers the diagonal density of the system which, in the uniform configuration studied here, exactly provides the linear density n. The momentum distribution n (k) which is a function of the momentum k, is the Fourier transform of the OBDM [38]. At very high temperatures, T T d , the system behaves as a classical gas and follows the Maxwell-Boltzmann (MB) statistics. The OBDM exhibits the Gaussian behavior and vanishes at distances much larger than the standard deviation σ = λ/ √ 2π which is proportional to the thermal wavelength λ, already defined in Eq. (10).
At small distance x, the OBDM can be expanded as a sum of analytic and non-analytic terms [85] and holding for any value of the interaction strength and temperature: The coefficients c i of the Taylor expansion of the analytic part of the OBDM are the corresponding moments of the momentum distribution n (k) [85], they diverge for i > 3 and the odd coefficients vanish c 1 = c 3 = · · · = 0 due to reflection symmetry. From the Hellmann-Feynman theorem [73] one finds: It is easy to show that the second coefficient, Eq. (20), can be rewritten in terms of the average kinetic energy H kin [38], see Appendix A. The non-analytic part of the OBDM expansion starts as |nx| 3 with the coefficient b 3 which is related to the high-momentum tail of n (k), as it depends on the Tan's contact C Both the internal energy and Tan's contact per particle, E/N and C/N respectively, in Eqs. (20)- (21) have been evaluated exactly with thermal Bethe Ansatz [13][14][15][16]. Equations (19)- (21) for the short-range expansion of the OBDM have been derived at T = 0 [58] and here they are generalized to a finite temperature whose dependence enters in E and C.
In higher spatial dimensions, where the Bose-Einstein condensate emerges below the critical temperature, the presence of Off-Diagonal Long-Range Order in homogeneous systems is manifested by a finite value of the OBDM at large distances, which corresponds to the condensate fraction [38]. In our system instead, the OBDM always vanishes at large distances, g 1 → 0, i.e., there is no Bose-Einstein condensation even at zero temperature.
Within the Luttinger Liquid theory, applicable at very low temperature k B T µ 0 and for any interaction strength γ, the OBDM exhibits a power-law decay [33,38,84,86]: holding for distances much larger than the healing length ξ = /( √ 2mv) and smaller than the thermal correlation radius x T = v/ (k B T ). The Luttinger Liquid parameter K = v F /v can be calculated exactly with Bethe Ansatz for any interaction strength γ, via the corresponding results of the speed of sound v. The power-law decay (22) has been shown to be accurate for any value of K by comparison with Diffusion Monte Carlo results at zero temperature [33], where quantum fluctuations depend explicitly on the value of v [38]. The vanishing asymptotic value of the one-body density matrix at large distances excludes the existence of Bose-Einstein condensation in infinite systems [87]. However, due to the smallness of the exponent 1/(2K) in the weakly-interacting GP regime, γ 1, the OBDM does not vanish at macroscopic distances |x| ξ, signalling the presence of a quasicondensate exhibiting features of superfluids [88].
In Figs. 7-9, we show with symbols the PIMC results of the one-body density matrix g 1 (x) for different values of the interaction strength γ and temperature. Solid lines denote the short-range expansion, Eq. (19), calculated from TBA quantities. In the upper panel of the sets of results for γ = 10 −1 and 1, the black dashed line denotes the Luttinger-Liquid quasi long-range decay, Eq. (22), valid in the limit of zero temperature. Instead, the power-law decay (22) is not reported for γ = 10 as in this strongly-interacting regime only temperatures beyond the validity of the LL theory have been explored. The DC regime of high temperatures reported in the lower panels of Figs. 7-9 allows a comparison with the Gaussian behavior expected for a classical gas, Eq. (18), and the corresponding predictions are shown with black dot-dashed lines. We notice that for any interaction strength, a better agreement with PIMC results is found at small distances by   pointing out that, in the DC regime, quantum statistics still play a role at large interatomic separation despite the high temperature. All analytical limits show a fair agreement with the corresponding PIMC results calculated at the same value of temperature. In Fig. 7, we show the OBDM in the weakly-interacting   GP regime with γ = 10 −1 . Each panel corresponds to different physical thermal regimes in order of increasing temperature values: hole anomaly (upper), intermediate (middle) and DC regime (lower) [14]. It is worth noticing that although the OBDM vanishes at large distances, in agreement with the absence of the Bose-Einstein condensation in one dimension, still the OBDM can be significantly different from zero at low temperatures, even at distances large compared to the mean interparticle distance x > n −1 . In the considered example, the coherence is preserved at distances of the order of hundreds mean interparticle distances. Such a behavior is often referenced as a quasicondensation and it is responsible for   making the Gross-Pitaevskii and Bogoliubov theories still applicable in 1D. Instead, as the temperature is increased, thermal fluctuations significantly suppress the coherence, making it disappear for x ∼ n −1 (middle panel) and even smaller distances (lower panel). In the latter case, the comparison with the Maxwell-Boltzmann prediction (black dot-dashed lines in the lower panel) makes it clear that at high temperatures the coherence is lost at distances of the order of σ.
The regime of the intermediate interaction strength, γ = 1, is the most difficult one for an analytic description as there is no separation of scales and the scattering length is of the same order as the mean interparticle distance a ∼ n −1 . Furthermore, at low temperature, the kinetic energy is of the same order as the potential interaction energy. The OBDM in this interaction regime is shown in Fig. 8. We observe that at low temperatures the OBDM is significantly different from zero up to distances of the order of tens of the interparticle distances, making the quasicondensate concept applicable even to this case. The coherence is instead lost when the temperature becomes comparable with the rest of the equal energy scales, i.e. for k B T ∼ k B T d = 2 n 2 / (2m) ∼ 2 /(ma 2 ). Thus, τ = T /T d = 1, separates the quantum low-temperature and decoherent classical regimes. In the upper panel of Fig. 8, we report results at very low temperature and large distances where the short-distance expansion, Eq. (19), is not shown.
The regime of large values of the interaction strength γ is very interesting as the Girardeau's mapping applies here, that is the wavefunction of bosons with infinite γ equals the absolute value of the wavefunction of the ideal fermions in 1D [40]. From this mapping, diagonal properties (pair correlation function, static structure factor, etc.) and the energy are the same in the two systems. Instead, the momentum distribution and the OBDM are manifestly different. Indeed, the momentum distribution of an ideal Fermi gas at zero temperature is given by the step function, n(k) IFG = ( n) −1 for |k| < k F (k F = πn is the Fermi momentum) and zero elsewhere. The corresponding OBDM is given by its Fourier transform and is an oscillating where the maxima of the Friedel oscillations occur at multiples of k −1 F . Also, the fermionic OBDM changes its sign (for example, for a displacement which interchanges two nearest fermions, n −1 < x < 2n −1 ). Obviously, the OBDM for bosons is different, Fig. 9, as g 1 (x) cannot be negative in the bosonic case, and a significant effort has been devoted to the calculation of the OBDM of the Tonks-Girardeau gas [58,87,[89][90][91][92]. Eventually, all coefficients of the zero-temperature short-range expansion, Eq. (19), have been calculated analytically in Ref. [85]. Even at low temperature (upper panel of Fig. 9), the coherence is rapidly lost, as The loss of coherence is exponentially fast at a finite temperature with the Fermi energy being the only relevant scale in the γ → ∞ limit. It is interesting to note, that the Maxwell-Boltzmann regime, Eq. (18), is reached only when the temperature is of the order of thousands of the quantum degeneracy temperature, see the lower panel of Fig. 9.

VII. EXPERIMENTAL CONSIDERATIONS
Cold atoms provide a powerful and clean platform for studying the physics of correlations in a widely tunable range of interaction and temperature. One-dimensional geometry can be realized in highly anisotropic trap configurations, where the confinement in the two radial y and z spatial directions is so strong that the motion of atoms is restricted to the x-axis. A reliable description of the system is then based on an effective 1D model, where both the temperature T and the chemical potential µ are much smaller than the first excited energy level of the trapping potential [93]. This condition can be expressed as k B T, µ ω ⊥ where ω ⊥ = ω y = ω z de-notes the harmonic trap frequency in the strongly confined directions.
Highly anisotropic geometries are created in optical dipole traps [82], two-dimensional optical lattices generating a set of 1D tubes [11,[41][42][43][44][94][95][96][97], on an atom chip [81,93,98,99] and with a single optical tube trap [12]. The latter two setups allow for the realization of a single 1D Bose gas which avoids the issue of ensemble averaging over many non-identical copies of the system taking place in lattice configurations. In addition, the exploration of even the classical gas regime of high temperature has been achieved in the optical tube trap [12], by keeping satisfied the condition k B T, µ ω ⊥ for the 1D geometry.
An important advantage of the one-dimensional geometry is that three-body losses are strongly suppressed [43], making it possible to reach the limit of a diverging coupling constant (unitarity). In three dimensions, unitary Bose gas is instead prone to instability due to the formation of Efimov trimers, which are absent in one dimension. A spatial uniform density can be achieved in a flatbox potential [100] which can be advantageous as compared to the measurements in the axially trapped configuration where the density is dependent on the spatial coordinate and, for example, the simple power law of the OBDM in the Luttinger Liquid regime gets totally blurred by the trap. The interaction strength γ ∼ 1/ (na) can be tuned at will in experiments by adjusting the linear density n via the strongly confining potential in the radial directions [37,42,44] or by applying Fano-Feshbach resonances to tune the 3D scattering length, and correspondingly the 1D length a in the Confinement Induced Resonance [96,101].
The measurement of temperature (thermometry) of onedimensional Bose gases can be achieved in various ways, which are often based on the comparison between experimental measurement statistics and theoretical predictions. For example, the temperature can be extracted from the pair correlation function of a single quasicondensate which is released and expands in the time-of-flight experiment along the 1D axis [102,103]. As a result, many repetitions of the same measurement are needed in order to obtain an accurate comparison with the theory. Instead, thermometry from a single absorption image measurement during the time-of-flight expansion has been recently made possible with neural network [104]. Compared to the standard method based on fitting the pair correlation function, the network can achieve the same precision needing only half the amount of experimental images.
In a 1D Bose gas, the pair correlation function g 2 (x) can be detected using spatially resolved in-situ single-atom counting as proposed in Ref. [9], standard absorption imaging [103] and observing intensity correlations in the interference pattern of two identical spatially displaced quasicondensates [105,106] during time-of-flight expansion. The local value g 2 (0) has been measured by photoassociation [44]. The OBDM g 1 (x) has been detected with the matter-wave interferometry of two copies of expanding quasicondensates [107].
Beyond the implementation in ultracold atoms, the 1D Bose gas has been experimentally realized in superconducting circuits where the OBDM and the pair correlation functions have been measured [4].

VIII. CONCLUSIONS
In conclusion, we have carried out a detailed study of correlation functions in a one-dimensional Bose gas with contact interactions at finite temperatures. The Path Integral Monte Carlo method is employed to perform exact calculations of the one-and two-body correlation functions, and the static structure factor. We find a good agreement with the predictions of a variety of perturbative and effective descriptions, as well as exact theories in regions of their applicability.
In particular, the short-range correlations have been extracted from several thermodynamic properties evaluated with Thermal Bethe Ansatz. The value of the pair correlation function g 2 (0) is related to the probability of two particles to overlap. It shows a competition of two opposite phenomena, that is the antibunching, i.e. tendency of two repulsive atoms to avoid each other, and a bunching effect, that is an enhanced probability of being close to each other due to thermal fluctuations. The short-range expansion of the One-Body Density Matrix contains both analytic and nonanalytic terms. The analytic contributions (x n ) are related to n-th moments of the momentum distribution and have non-vanishing coefficients for even powers, the leading one (x 2 ) can be expressed in terms of the average kinetic energy. The nonanalytic contributions come as an expansion of powers of |x|, being the leading one proportional to |x| 3 whose coefficient depends on the Tan's contact and thus captures the amplitude of the 1/k 4 -decay of the momentum distribution. By comparison with PIMC results, we show that the short-range expansion of the OBDM is valid for any value of the interaction strength and temperature.
At zero temperature, the coherence is lost in a slow longrange power-law decay for the OBDM allowing for the application of condensate-based theories (Gross-Pitaevskii, Bogoliubov, etc.) in the weakly-repulsive regime. At temperature higher than the hole-anomaly threshold, the Decoherent Classical perturbative theory well reproduces PIMC results for the pair correlation function. It captures both the bunching effect, described by a Gaussian term whose width is provided by the de Broglie thermal wavelength, and the antibunching of strength proportional to the interaction parameter γ, Eq. (10). By further increasing the temperature, correlations are no longer affected by Bose or Fermi statistics but rather are described by the Maxwell-Boltzmann Gaussian decay which is independent on γ.
The low-momentum behavior of the static structure factor is well captured by a single-mode quasiparticle approximation below the hole-anomaly temperature. A linear-phonon single-mode dispersion lies in the base of the Luttinger Liquid theory whose applicability is confirmed for any interaction strength. In the weakly-repulsive regime, the quasiparticle excitation is provided by the Bogoliubov spectrum which is linear at small momentum but deviates at larger momenta. The zero-momentum static structure factor S(0) is related to the isothermal compressibility and its value can be extracted from the Thermal Bethe Ansatz. It quantifies the interplay between the quantum and thermal fluctuations and shows that the thermal effects are more important for weak interactions.
Our exact PIMC and thermal Bethe-Ansatz results do not agree with some known analytical limits of the pair correlation function g 2 (x) [9,10] and of its local value g 2 (0) [8], respectively, the latter of which is a thermodynamic property. Improved analytical descriptions have been recently developed and they exhibit an excellent matching with the exact TBA method for several thermodynamic quantities [13,14,50]. We expect that our PIMC calculations of the correlation functions will serve as a universal benchmark for future theories and experimental measurements. Our results may stimulate further theoretical and experimental investigations aimed at the complete understanding at the microscopic level of 1D Bose gases through quantum many-body correlations which allow for the identification of different quantum regimes. The PIMC results for the momentum distribution will be presented in an upcoming publication [108].
Looking forward, interesting perspectives include the calculation of n-body correlations [109] and the study in optical lattices [110]. Particularly appealing is the extension to the harmonically axially-trapped configuration [28,60,[111][112][113][114][115] which is suitable for the investigation of breathing modes [116][117][118] whose collective frequencies change by crossing different quantum regimes [119]. Recently, a temperatureinduced hydrodynamic-to-collisionless transition has been predicted with the emergence of the excitation of two different frequencies in the time-evolution of breathing modes in both the GP and TG limits [120]. Such a transition may shed light on the absence of thermalization in 1D Bose gases even after thousands of interatomic collisions [95].
Our results can be extended to systems with positive 1D s−wave scattering length a > 0: (i) gases with short-range interactions like the strongly correlated metastable state (super Tonks-Girardeau gas) [78,80,121] and the hard-core model [39,[122][123][124][125], (ii) finite-range interacting systems such as dipolar [126,127] and Rydberg atoms [128], 4 He liquid [129], 3 He gas [130], and a single chain of 4 He atoms confined in a nanopore even with disorder [131]. Other interesting extensions of our work include spin models [132,133] and multicomponent systems [134][135][136][137]. Our predictions are important for the exploration of the properties of impurities immersed in helium [138] and in a Bose gas [139,140] as a function of the interaction strength of the bath and temperature. Possible extensions of our work also include quantum liquids in binary bosonic mixtures at finite temperature and 1D geometry [141] which can be realized in current experiments [142]. Finally, our findings can be also applied to 1D liquids of 4 He which have been recently observed [143] in a wide range of interaction and temperature regimes.
Eq. (A1) relates the Tan's contact, Eq. (4), the pressure, Eq. (2) and the total internal energy E = H kin + H int where the average interaction contribution can be calculated from the Hellmann-Feynman theorem [73,120]: and where we have used Eq. (4). By combining Eqs. (A1)-(A2), one calculates the average kinetic energy: By using Eq. (A2) in Eq. (20) of the main text, one finally finds [38]: where the average kinetic energy per particle can be expressed as H kin /N = 2 k 2 / (2m) with k 2 = N −1 ´+ ∞ −∞ dkn(k)k 2 calculated from the momentum distribution n (k) which is a function of the momentum k [38].