Equations of state from individual one-dimensional Bose gases

We trap individual 1D Bose gases and obtain the associated equation of state by combining calibrated confining potentials with in-situ density profiles. Our observations agree well with the exact Yang-Yang 1D thermodynamic solutions under the local density approximation. We find that our final 1D system undergoes inefficient evaporative cooling that decreases the absolute temperature, but monotonically reduces a degeneracy parameter.


Introduction
Strongly interacting systems are ubiquitous in modern physics, from astrophysical objects such as neutron stars to the myriad of correlated electron systems in condensed matter. Experimental developments in ultracold atomic physics enable multiple avenues to explore interacting quantum matter, for example with optical lattices [1], Feshbach resonances [2] or mediated long-range interactions [3]. Furthermore, tailored potentials can reduce the effective dimensionality of cold atomic gases; for example, a 2D optical lattice can partition a 3D system into an array of 1D gases [4,5,6]. Remarkably, the theory of a 1D Bose gas (1DBG) with contact repulsive interactions is exactly solvable at all temperatures [7,8], making it an ideal system to benchmark experiment against theory.
In cold atomic gases, the ≈ 5 nm range [2] of the interatomic potential is vastly smaller than the 100 nm interatomic spacing, hence interactions are well described by a local contact potential with strength g. This gives the 1D Hamiltonian for N interacting bosons of mass m which in the absence of a potential, i.e. V (z) = 0, becomes the Lieb-Liniger [7] Hamiltonian. Lieb and Liniger showed that eigenstates of this Hamiltonian are parametrized by the dimensionless interaction parameter γ = mg/ 2 n, where n is the density. Here the relevance of interactions increases with decreasing density. For γ 1 mean-field theory accurately describes the system, while for γ 1 the atoms strongly avoid one another and behave much like weakly interacting fermions. C. N. Yang and C. P. Yang extended this solution to nonzero temperature [8] and cold atom experiments have validated the accuracy of the "Yang-Yang" thermodynamics [9,10,11].
Here we study the physics of individual 1DBGs using 87 Rb atoms in an optical "tube trap" (Figure 1(a)) and benchmark the thermodynamic equation of state (EoS) against Yang-Yang thermodynamics. Our individual-system realization bridges an existing gap in experiments, on the one hand avoiding the issue of ensemble averaging present in realizations using optical lattices [4,5,6,10] and on the other hand enabling the future study of 1D multi-component systems not viable using magnetic confinement potentials [12,13].
The physics of 1DBGs can be divided into three qualitative regimes [14] shown in Figure 2. For sufficiently high temperatures (green region) the EoS is dominated by thermal effects and approaches that of a non-interacting classical gas. Below the degeneracy temperature T d = 2 n 2 /2mk B , where k B is Boltzmann's constant, and for weak interactions (γ 1), the reduced thermal fluctuations allow Bose statistics to weigh in [15], creating a phase fluctuating degenerate gas. For lower temperatures where T /T d < 2γ (blue region), the thermal energy falls below the chemical potential and the system is well described by the Gross-Pitaevskii equation (GPE). In contrast, for systems with large interactions (γ 1) the EoS resembles that of an ideal Fermi gas [10], with the formation of a Fermi surface for T < T d (red region). Yang-Yang thermodynamics provides EoS's encompassing all of these regimes, relating quantities like the particle, entropy and pressure densities to the chemical potential µ and temperature T , e.g., n(µ, T ).
In trapped systems, the confining potential V (z) ≥ 0 can often be treated as an inhomogeneous chemical potential µ(z) = µ 0 − V (z) which allows for multiple regimes to be present in a single 1DBG. We define µ 0 as the local chemical potential at V (z) = 0.  Figure 2. Three regimes can be identified in the γ, T /T d parameter space which correspond to the strongly interacting degenerate (red), weakly interacting degenerate (blue) and non-degenerate (green) regimes. Specific experimental realizations of inhomogeneous systems cover a given range of these parameters (white circles).
This can be quantitatively understood [16] within the local density approximation (LDA) allowing the density profile n(z, T ) to be interpreted as an EoS n(µ(z), T ). As a consequence, the EoS can be experimentally determined from a well-calibrated trapping potential V (z) and the observed density profiles.
We extract this n(µ(z), T ) EoS from in-situ absorption images [17,18,19] of individual systems (Figure 1(b)), eliminating ensemble averaging. Because we obtain the 1D density directly, we do not apply the inverse Abel transform [18] thereby avoiding added noise. We benchmark our measurement against the Yang-Yang EoS (Figure 1(c)), from which other thermodynamic quantities become readily available (e.g. free energy, entropy and pressure).
This manuscript is organized as follows: in Sec.
(2), we describe our experimental setup and data acquisition protocol; in Sec. (3), we address the different calibration aspects of our analysis; in Sec. (4), we discuss the results; and in Sec. (5), we conclude.

Experiment
We prepare cold atoms beginning with a magneto-optical trap followed by vertical magnetic transport [20] into a magnetic quadrupole trap, ultimately loading a 1064 nm crossed optical dipole trap [21,22]. This gives N ≈ 2 × 10 5 atom Bose-Einstein condensates (BECs) of 87 Rb in the 5 S 1/2 |F = 1, m F = 0 hyperfine ground state with ≈ 70 Hz mean trapping frequencies. The atoms are then transferred into the composite high aspect-ratio trap shown in Figure 1(a). This trap includes a reddetuned (λ G = 1064 nm) Gaussian beam along e x with waist w G = 203(2) µm providing reduced longitudinal confinement owing to its larger waist as compared to the ≈ 70 µm crossed dipole beam waist. A transverse "tube trap" along e z is provided by a blue-detuned (λ LG = 532 nm) Laguerre-Gauss (LG 01 ) beam, tightly focused to a waist of w LG = 5.6(5) µm. In our standard configuration these beams have powers P G = 0.8(1) W and P LG = 1.0(1) W , giving a peak transverse trapping frequency ω ⊥ /2π = (ω x ω y ) 1/2 /2π = 17(2) kHz.
The transverse zero-point energy from ω ⊥ produces an anti-confining potential along e z due to the divergence of the LG beam. The anti-trapping potential shown in green in Figure 3 significantly alters the overall longitudinal potential where ω (0) ⊥ /2π denotes the peak transverse trapping frequency; z a is the center of the anti-trap; z R is the Rayleigh range of the LG beam; and V 0 is an energetic offset chosen such that the minimum of the potential is zero. The black shaded curve in Figure 3 shows the combined anharmonic potential of the longitudinal trap (red curve) and the anticonfining potential. Small amplitude longitudinal dipole oscillations in the combined potential have frequency ω z /2π = 12.1(2) Hz.   LG beam from zero in 250 ms until the tube trap can suspend atoms against gravity. Because the 3D system is always larger than 30µm, the ≈ 5 µm LG beam only captures a small fraction of the initial 3D ensemble. (ii) We then lower the intensity of the crossed dipole trap in 250 ms, allowing the atoms outside the tube trap to fall away. (iii) We then ramp up the final 1064 nm longitudinal trap in 250 ms. (iv) In the final 250 ms we simultaneously increase the intensity of the LG beam to its final value while removing the crossed dipole trap.
These 250 ms ramps were chosen to be adiabatic with respect to all the confining potentials. Monopole and dipole modes of the 1DBG can be induced by both beam misalignment and excessive ramp rates in this scheme. Our ramp times were chosen to mitigate these residual excitations. Intensity (a .u.) We control the temperature of the 1DBG by varying the temperature T 3D of the initial 3D Bose gas. We tune T 3D by adjusting the depth of the crossed dipole trap, covering the range from T 3D = 34 nK to 320 nK, with an observed BEC transition at T c 3D ≈ 160 nK. We determine T 3D with time-of-flight (TOF) measurements. The number of atoms in the 1DBG increases with decreasing T 3D due to the increasing 3D density as T 3D falls. Gravitational sag is a complicating factor: as the crossed dipole trap decreases, the 3D ensemble lowers due to gravity, but the vertical alignment of the tube trap does not. We mitigate this effect by increasing the crossed dipole power after the final evaporation such that the crossed dipole potential is in a fixed vertical position prior to loading the tube trap.

Imaging
We derive the density n(z) from in-situ absorption images. Our imaging system has a resolution of 1.85(5) µm and magnification that maps one 5.6 µm sensor pixel to 0.84(1) µm in the object plane. In preparation for imaging, we apply a 20 µs repump pulse to transfer the 5 S 1/2 |F = 1, m F = 0 atoms into the 5 S 1/2 |F = 2 hyperfine manifold. We then image [23] with a circularly polarized λ p = 780 nm probe beam resonant with the 5 S 1/2 |F = 2 to 5 P 3/2 |F = 3 transition for 20 µs at an average intensity of I = 2.5I sat , where I sat = 1.67 mW/cm 2 is the saturation intensity of the resonant atomic transition. An image of the probe beam following absorption I a , the probe without atoms present I p , and a dark frame with no probe present I d , are each recorded on a charge-coupled device (CCD) camera. From these images we obtain the absorbed fraction f = (I p − I a )/(I p − I d ). For each set of experimental parameters we repeat the experiment for ≈ 100 realizations.
Our image analysis is a multiple step process. We first preprocess the raw images to correct for background artifacts and improve the signal-to-noise ratio (SNR) by a factor of ≈ 10. We then extract the linear densities using an absorption model that includes a modest Lamb-Dicke suppression. As compared to a naïve model, n(z) increases by as much as 30%. This process leaves our qualitative results unchanged.
We reconstruct an optimal I opt p for each I a as a linear sum of I p from all realizations by minimizing the sum squared difference between I opt p and I a away from the atoms [24,25]. This reconstruction reduces fringing due to vibrational motion that occurred between acquiring I a and I p , along with shot noise present on each I p . We use similar techniques to remove a systematic difference in dark counts between I a and I p , as well as to account for structured read-out noise. We then compute mean absorbed fractionsf over each set of experimental parameters, and use uncentered principal component analysis to further suppress shot noise. Fromf and a detector calibrated [26,27,17] in units of I sat , we compute 'naïve' optical depths using the standard solution to the Beer-Lambert (BL) law [27], which we sum along columns to produce 'naïve' linear densities.
We take into account the fact that the transverse extent of our atom cloud (for a tube trap the radial confinement gives an extent of /mω ⊥ ≈ 110 nm) is below the resolution of both our imaging system and the optical scattering length 3λ 2 p /2π 2 ≈ 300 nm. We further incorporate the transverse diffusive motion that atoms undergo during the imaging pulse. Each of these effects violates the assumptions underlying the BL law. In comparison with the naïve BL law, our model for the density agrees at low density but deviates up to 30% at higher densities. This process is described in greater detail in the supplementary material.

Results
The results of our image processing are 1D density profiles n (j) (z) confined in the same trapping potential but with one of 24 different initial conditions labeled by j. In the local density approximation we expect that these density profiles can result from Yang-Yang thermodynamics. For each j, both the temperature T (j) and the overall chemical potential µ (j) 0 are in principle unknown because of the lack of suitable reservoirs. As a result we obtain these quantities from fits to the Yang-Yang EoS and assess their validity in terms of the reduced chi-squared.
For each j, the Yang-Yang EoS predicts the complete density profile as shown in Figure 1(c) with just two free parameters T (j) and µ (j) 0 . We constrain the fit to the trapping region between the local maxima of V (z). The potential is parametrized by the common set of parameters shown in Table 1. We include some of these as additional parameters in our fit shared between all j. In Table 1 we show the calibrations by other measurements along with their uncertainties; these are provided as initial guesses and bounds to the Yang-Yang fit. An additional uncalibrated parameter δz accounts for a tiny displacement of the 1DBG relative to the center of V (z) for times following the loading protocol. The inclusion of δz leaves the main results unchanged and its value lies within the relative alignment uncertainty of the trap centers. Different combinations of fixed parameters have no qualitative effect on the results. The third column in Table 1 shows the potential parameters derived from the Yang-Yang fit. We evaluate the goodness-of-fit with the reduced chi-squared χ 2 ν ≈ 1.5. Lastly, the residuals of the fit show systematic variations in I a , which are reflected by the gray pixels in Figure 1(c). Figure 5(a) shows the reduced density versus reduced chemical potential for two initial conditions, each plotting different cuts in the EoS n(µ, T ). The continuous curves in Figure 5(a) represent the Yang-Yang model with the T and µ 0 from our fits. For small chemical potential these density profiles are well described by the EoS of a noninteracting Bose gas while for µ > 0 they approach the predictions of GPE mean-field theory [14]. The Yang-Yang EoS accurately predicts both regimes. A sharp eye observes an apparent hysteresis loop visible in the trace labeled by A, this results from the spatial dependence of g that follows ω ⊥ . As shown in Figure 3, ω ⊥ is slightly off-centered, ultimately resulting in the observed behavior. The scattered white dots on Figure 2 represent these two traces in the γ, T /T d plane. These data are shown to be either in the interacting regime or below degeneracy, but not in both. Figure 5(b) summarizes the outcome of all our Yang-Yang fits in which we varied T 3D , the initial 3D temperature (left panels); or the hold time t in the 1D trap for our lowest T 3D (right panels) cloud. In the top left panel, we observe that as a function of decreasing T 3D , the 1D temperature T first remains constant and then counterintuitively increases. In contrast, as shown in the bottom left panel, the degeneracy parameter defined as T /T (0) d , is a monotonically increasing function of T 3D showing how the more degenerate 3D clouds result in more degenerate 1DBGs.
For the lowest achievable T 3D and as a function of hold time, we see that both the total atom number N and the 1D temperature T drop (top right panel in Figure 5(b) and Figure 5(c) respectively), yet the 1DBG doesn't become more degenerate (bottom right panel in Figure 5(b)). The simultaneous drop in T and N is consistent with evaporative cooling along the longitudinal axis of the tube-trap, which has depth of ≈ 700nK. The inability of such evaporative cooling to increase or even maintain degeneracy results from the slower relative decrease in T with respect to T d as driven by the atom number loss.
We explore the character of this loss by modeling the atom number decay with a model including one-body loss and three-body loss from photon scattering, background gas collisions and inelastic three-body collisions [5]. Figure 5(c) shows the measured atom number N (blue diamonds). We fit the decay model to the observed number (dark blue curve), giving a one-body loss coefficient K 1D 1 = 0.108(2)s −1 as well as a three-body loss coefficient K 1D 3 = 4.36(7) × 10 −29 cm 6 s −1 . The value of K 1D 1 is consistent with the combined vacuum-limited lifetime and estimated off-resonant scattering rate from the static dipole potentials. In contrast, the three-body loss coefficient from our fit is in excess of the intrinsic 3D three-body loss coefficient K 3D 3 = 5.8(3) × 10 −30 cm 6 s −1 [28] by a factor of ≈ 7.5. We attribute this enhancement to the difference in the three-body correlation function g (3) [5,14] from a purely coherent sample. The observed cooling is consistent with initial rapid evaporation as atoms with sufficient kinetic energy [29] overcome the longitudinal barrier of V (z).

Conclusions
We realized individually trapped 1DBGs in a crossed dipole trap formed by a bluedetuned LG 01 beam and a red-detuned Gaussian beam. We benchmarked the EoS computed from Yang-Yang thermodynamics against the measured density profiles. We found that evaporative cooling along the edges of the tube trap took place although this did not maintain or increase the system's degeneracy. Our approach enables future exploration of spinor 1DBGs associated with multi-component physics [30,31], including spin-orbit coupling [32]. This therefore presents a promising venue to study the limits of strongly interacting 1D systems in and out of equilibrium.

Acknowledgments
We thank Max Schemmer for insightful discussions. We acknowledge the support for this work provided by the AROs Atomtronics MURI, NIST and the NSF through the PFC at the JQI.
We then compute the mean dark frame I d over all shots, and perform principal component analysis (PCA) on the set of all dark frames, with two PCA eigenvectors revealing a source of correlated dark noise in the form of spatially sinusoidally varying dark counts with a different phase offset for each shot, which we also observed to be present in the PCA eigenvectors of the probe images (as eigenvectors four and five). We project each absorption image onto these eigenvectors I pca4 p and I pca5 p of the probe images in order to determine coefficients c pca4 and c pca5 . A reconstructed dark field image I recon d is then computed for each shot as

Absorbed fraction and saturation parameter
Absorbed fraction f and saturation parameter S images are then computed for each shot as where I sat is the saturation intensity in camera count units. The 'naïve' optical depth for a single shot can then be computed as: where α ≡ σ 0 /σ eff = I sat eff /I sat is the empirically determined ratio between the ideal twolevel and effective scattering cross sections due to imperfect polarization and magnetic field orientation. The average naïve optical depth over all repeated shots for each set of experimental parameters is computed as: where the mean product of absorbed fraction and f and saturation parameter S are taken over only the repeated shots for one set of experiment parameters, and where we compute the mean absorbed fraction within the log rather than the mean of the entire log term in order to avoid the systematic error that results from taking the mean of a nonlinear function of noisy data. This naïve optical depth is not accurate across our entire dataset due to our 1D system being narrower transversely than both the optical scattering length 3λ 2 p /2π 2 and our imaging resolution, both of which are violations of the assumptions of the Beer-Lambert law. We continue to compute further reconstructions of this naïve optical depth only for comparison with modeled linear densities which include a correction to the Beer-Lambert law to account for this, presented further below.

Appendix A.3. Dimensionality reduction of absorbed fraction
We dimensionally reduce the mean absorbed fractions f of each set of experiment parameters in order to reduce the effect of shot noise on column sums of the data. Since the point spread function resulting from diffraction in our imaging system is fixed from shot to shot, this has the effect of projecting the measured absorbed fractions onto the empirically observed point spread function and its most common variations, suppressing spurious apparent absorption due to shot noise in regions where the point spread function results in little absorption.
The dimensionality reduction proceeds as follows. First we crop each mean absorption image f (x, y) to the 75-pixel high region of interest that entirely contains our imaging system's point spread function to form f ROI (x, y). Then, for each x position x i in the image, we extract the vertical lines of f ROI (x, y) from all sets of experiment parameters, at that x pixel and the nearest four x pixels. Treating each vertical line as a vector § f ROI (x i ) = f ROI (x i , y), we obtain the set of vectors {f ROI (x i+j ), −2 ≤ j ≤ 2} and perform unentered PCA [33], keeping only the first four normalized eigenvectors {v n (x i ), n ∈ [1, 4]}. We then project the (also cropped to the ROI) vertical lines f ROI (x i ) of the absorbed fractions for each individual shot at the original x position onto the subspace spanned by these vectors: Thus we dimensionally reduce each vertical slice (within the ROI) of each shot's absorbed fraction onto a basis of four basis functions chosen by unentered PCA of the vertical slices of all mean absorbed fractions at that x pixel and the nearest four other x pixels. We then compute mean, dimensionally reduced absorbed fractions f red (x, y) for each set of experiment parameters within the ROI (Hereafter any images mentioned should be assumed to be cropped to the ROI). We observe that to within numerical rounding error, it makes no difference whether the mean absorbed fractions are dimensionally reduced into this space, or the individual absorbed fractions are, before being averaged together again. We do the latter in order to provide a statistical uncertainty estimate in the mean absorbed fraction of a given set of experiment parameters as where σ is the standard deviation over all repeated shots for one set of experiment parameters, {A red (x, y)} is the set of all dimensionally reduced absorbed fractions for those shots, and N is the number of repeated shots for that set of experiment parameters. § The mixed vector vs. function-of-space notation here is because we are treating the vertical slices of the images as vectors, performing dimensionality reduction on a slice-by-slice basis, whereas the x coordinate is merely a label selecting which vector we are referring to.

Saturation parameter at the position of the atoms
Due to the point spread function of our imaging system being larger than the vertical extent of our atomic cloud, the saturation parameter at the location of apparent absorption (after diffraction) does not correspond to the saturation parameter at the actual location of the atoms, which is where saturation effects are relevant. We estimate a saturation parameter S 0 (x) for each set of experiment parameters at the estimated y position of the atoms by interpolating the mean saturation parameter for that set of experiment parameters to the y position where there is maximum apparent absorption over all shots. The y position of maximum apparent absorption at each x position is taken to be a quadratic fit y 0 (x) = ax 2 + bx + c with parameters determined by maximizing the total absorption over all shots: where interp y is a one-dimensional spline interpolation function interpolating in the y direction only. The estimated saturation parameter at each x position for each set of experiment parameters is then

. Naïve linear density
We can now compute an improved naïve optical depth for each set of experiment parameters using the dimensionally reduced absorbed fractions and interpolated saturation parameter as: and then compute a naïve linear density n naïve (x) at each x position by dividing by the cross section and integrating along y: where ∆y is the pixel size.

Appendix A.5. Modeled linear density
We face three related problems in computing the column density n col (x, y) given a measured absorbed fraction f (x, y) and saturation parameter S(x, y) via the solution to the Beer-Lambert law [27]: 14) The first problem is that we do not measure f (x, y) directly-we measure it only after it has diffracted in the y direction, a difference which is not negligible given the size of our atom cloud in that direction. The second problem is that atoms do not only absorb light at their exact location in space, rather they absorb it from a surrounding region of space with cross sectional area given by the absorption cross section σ 0 [34]. The final problem is that our cloud is so small in the y direction that diffusion of atoms during imaging may not be negligible. These latter two problems mean that we cannot infer n col (x, y) from the usual solution to the Beer-Lambert law, we can only determine the convolution of n col (x, y) with some absorption profile g(y) that takes into account both the finite absorption region and the diffusion of atoms from their initial positions in the y direction, the direction in which g(y) is not small compared to our atom cloud. With this in mind, the solution to the Beer-Lambert law can be modified to read where the convolution is only along the y direction. Atomic diffusion and diffraction imply that the only quantity we have experimental access to is that is, we only observe a time average of absorption over the imaging pulse time τ , and we only observe the integral of the undiffracted absorbed fraction, since diffraction preserves this integral. If the second term of the solution to the Beer-Lambert law dominates, then the naïve linear density is accurate, since all three of diffraction, diffusion and convolution preserve integrals of the absorbed fraction. It is only the log term that causes a problem, since its integral is not conserved under diffraction.
Given a model for (n col * g)(x, y, t) with a single parameter n(x) for the linear density at each x position and a measurement f meas (x) for each x position, we can invert (A. 16) and (A.15) to obtain the linear density, under the assumptions of the model.
Our model is the following: The absorption profile g(y) is approximated by a Gaussian with unit integral and standard deviation σ y (t) equal at t = 0 to the optical scattering length σ 0 /π = 3λ 2 p /2π 2 and increasing due to atomic diffusion as time elapses. Since the atom cloud is narrower than this absorption profile, we approximate the convolution (n * g)(x, y, t) as: where σ 2 y (x, t) is increasing due to momentum diffusion: where the mean squared y velocity σ 2 vy is given by the scattering rate R scat and recoil velocity v rec : which is approximating isotropic scattering such that the per-scattering-event expected squared change in velocity is v 2 rec /3. The scattering rate, ignoring Doppler shifts away from resonance, is, in terms of the saturation parameter S: .
Putting this all together, the modeled y variance of the absorption profile is: Over the duration of our imaging pulse, the diffusion described by the second term results in an increase in the absorption profile's standard deviation by ≈ 30% compared to the effect of the non-zero optical scattering length alone. Using our absorption model (A.17) with an absorption profile with y variance given by (A.21), and saturation parameter S(x, y) given by our estimate S 0 (x), our modified Beer-Lambert law solution (A.15) becomes: which, if numerically inverted, defines a function that takes only n(x) as input and returns f (x, y, t) at any given time. Numerically integrating the result in t and y as per (A.16) extends this function into one which takes only n(x) and returns the expected f meas (x) for that linear density. Numerically inverting this function then yields our final aim, of a function that takes as input f meas (x) from our data and outputs a value of n(x) for the linear density implied by the measured data and the model. We perform the above computationally nontrivial calculation to extract modeled linear densities from our dimensionally-reduced mean absorbed fractions and interpolated saturation parameters for each set of experiment parameters.
The naïve and modeled linear densities agree at low densities but disagree by up to 20 percent at higher densities, with the naïve method underestimating linear densities compared to those obtained using the absorption model.

Appendix B. Yang-Yang thermodynamics
We use the Yang-Yang model [8] to describe our data. The exact diagonalization of the underlying Hamiltonian is carried out with the use of the thermodynamic Bethe ansatz (TBA) (T > 0 Bethe ansatz). From the TBA the following set of first-order integral equations can be derived where m is the mass, k B is the Boltzmann constant, T is the temperature, µ the chemical potential and c = mg/ 2 is the interaction wavenumber. Both k and q label momenta. The above equations can be solved recursively to compute n, the linear density, given the values for µ, T and g, the chemical potential, temperature and coupling constant. We implement a numerical solver for the YY equations within the LDA that takes the parameters µ, T as its primary input and computes c and g by using the appropriate values of the transverse trapping frequency f ⊥ , the 3D scattering length a 3D and mass m. We recursively solve for (k) and f (k) from which we ultimately compute the density n. We transform all the momentum and energy quantities We denote the Lieb-Liniger kernel (a normalized Lorentzian) as L(c, q). Our numerical solver performs a k-space convolution using the scipy.signal.fftconvolve method to evaluate the integrals. For this we use a N k = 1024 points grid that covers the range k = [−10 k th , 10 k th ], where k th = 2mk B T / 2 is the thermal wavenumber. We initialize 0 (k) = k 2 − µ and iterate over the following recursive relation j+1 (k) = 0 (k) − L(c 0 , q) ln(1 + e − j (q) ), (B.8) where denotes the Fourier convolution operator. Once the convergence condition i ( i,j+1 − i,j ) 2 /N k < tol is satisfied, we solve for f (k) with an initial guess f 0 (k) = [2π (1 + e (k) )] −1 and the recursive relation f j+1 (k) = f 0 (k) + L(c 0 , q) f j (q), (B.9) from which we get to evaluate (B.3). After watching all the unit conversions we get the linear density in particles per meter.