PAPER • OPEN ACCESS

Equations of state from individual one-dimensional Bose gases
To cite this article: F Salces-Carcoba et al 2018 New J. Phys. 20 113032

View the article online for updates and enhancements.

This content was downloaded by ispielman from IP address 129.2.180.226 on 26/11/2018 at 14:35

New J. Phys. 20 (2018) 113032

https://doi.org/10.1088/1367-2630/aaef9b

PAPER

Equations of state from individual one-dimensional Bose gases
OPEN ACCESSF Salces-Carcoba1, C J Billington1,2, A Putra1, Y Yue1, S Sugawa1,3 and I B Spielman1
RECEIVED

1

17 August 2018
REVISED2

12 October 20183

ACCEPTED FOR PUBLICATION

9 November 2018
PUBLISHED

23 November 2018

Joint Quantum Institute, National Institute of Standards and Technology, and University of Maryland College Park, Gaithersburg, MD
20899, United States of America
School of Physics and Astronomy, Monash University, Victoria 3800, Australia
Present address: PRESTO, Japan Science and Technology Agency (JST), Saitama 332-0012, Japan, and Graduate School of Science, Kyoto
University, Kyoto 606-8502, Japan.

E-mail: pacosalces@gmail.com
Keywords: quantum gases, atomic physics, one-dimensonal Bose gas, equation of state

Original content from thiswork may be used under
the terms of the Creative
Commons Attribution 3.0
licence.

Abstract
We trap individual 1D Bose gases and obtain the associated equation of state by combining calibrated
Any further distribution of conﬁning potentials with in situ density proﬁles. Our observations agree well with the exact Yang–
this work must maintain
Yang 1D thermodynamic solutions under the local density approximation. We ﬁnd that our ﬁnal 1D
attribution to theauthor(s) and the title ofsystem undergoes inefﬁcient evaporative cooling that decreases the absolute temperature, but
the work, journal citation
monotonically reduces a degeneracy parameter.
and DOI.

1. 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–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⎤⎡  2 ¶2()
H = å ⎢+Vz
⎥ + g å d (z i - z j ) ,i
2⎦i ⎣ 2m ¶z i
i¹j

(1)

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 dimensionlessinteraction parameter g = mg  2n , where n is the density. Here the relevance of interactions increases withdecreasing density. For γ=1 mean-ﬁeld theory accurately describes the system, while for γ?1 the atoms
strongly avoid one another and behave much like weakly interacting fermions. Yang and Yang extended thissolution to non-zero temperature [8] and cold atom experiments have validated the accuracy of the ‘Yang–Yang’
thermodynamics [9–11].
Here we study the physics of individual 1DBGs using 87Rb atoms in an optical ‘tube trap’ (ﬁgure 1(a)) andbenchmark the thermodynamic equation of state (EoS) against Yang–Yang thermodynamics. Our individualsystem realization bridges an existing gap in experiments, on the one hand avoiding the issue of ensemble
averaging present in realizations using optical lattices [4–6, 10] and on the other hand enabling the future studyof 1D multi-component systems not viable using magnetic conﬁnement potentials [12, 13].
© 2018 The Author(s). Published by IOP Publishing Ltd on behalf of Deutsche Physikalische Gesellschaft

New J. Phys. 20 (2018) 113032

F Salces-Carcoba et al

Figure 1. (a) Schematic of the trap formed by a tightly focused blue-detuned laser beam in the LG01 mode propagating along ez and a
red-detuned Gaussian laser beam propagating along ex. (b) In situ absorption image of the trapped 1DBG. (c) Linear densities from asingle 1DBG (white circles), the average of 100 realizations (green circles), and the ﬁt of the Yang–Yang model to the average (blacksolid curve) giving T=152 nK and μ0/ÿ=2π×1.92 kHz, along with the residuals (gray squares).

Figure 2. Three regimes can be identiﬁed in the γ, T/Td parameter space which correspond to the strongly interacting degenerate(red), weakly interacting degenerate (blue) and non-degenerate (green) regimes. The different qualitative regimes are approximately
delineated by the black lines, but these are simply crossovers. Speciﬁc experimental realizations of inhomogeneous systems cover agiven range of these parameters (white circles).

2New J. Phys. 20 (2018) 113032

F Salces-Carcoba et al

The physics of 1DBGs can be divided into three qualitative regimes [14] shown in ﬁgure 2. For sufﬁciently
high temperatures (green region) the EoS is dominated by thermal effects and approaches that of a noninteracting classical gas. Below the degeneracy temperature Td =  2n2 2mkB, where kB is Boltzmann’s constant,and for weak interactions (γ = 1), the reduced thermal ﬂuctuations allow Bose statistics to weigh in [15],creating a phase ﬂuctuating degenerate gas. For lower temperatures where T/Td<2γ (blue region), the thermalenergy 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<Td (red region). Yang–Yang thermodynamics provides EoS’s
encompassing all of these regimes, relating quantities like the particle, entropy and pressure densities to thechemical potential μ and temperature T, e.g. n(μ, T).
In trapped systems, the conﬁning 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 deﬁne μ0 as
the local chemical potential at V(z)=0. This can be quantitatively understood [16] within the local density
approximation (LDA) allowing the density proﬁle n(z, T) to be interpreted as an EoS n(μ(z), T). As aconsequence, the EoS can be experimentally determined from a well-calibrated trapping potential V(z) and theobserved density proﬁles.
We extract this n(μ(z), T) EoS from in situ absorption images [17–19] of individual systems (ﬁgure 1(b)),
eliminating ensemble averaging. Because we obtain the 1D density directly, we do not apply the inverse Abeltransform [18] thereby avoiding added noise. We benchmark our measurement against the Yang–Yang EoS(ﬁgure 1(c)), from which other thermodynamic quantities become readily available (e.g. free energy, entropy
and pressure).
This manuscript is organized as follows: in section 2, we describe our experimental setup and dataacquisition protocol; in section 3, we address the different calibration aspects of our analysis; in section 4, we
discuss the results; and in section 5, we conclude.

2. Experiment
We prepare cold atoms beginning with a magneto-optical trap followed by vertical magnetic transport [20] intoa magnetic quadrupole trap, ultimately loading a 1064 nm crossed optical dipole trap [21, 22]. This givesN≈2×105 atom Bose–Einstein condensates (BECs) of87Rb in the 5 S1 2 ∣F = 1, mF = 0ñ hyperﬁne groundstate with ≈70 Hz mean trapping frequencies. The atoms are then transferred into the composite high aspectratio trap shown in ﬁgure 1(a). This trap includes a red-detuned (λG=1064 nm) Gaussian beam along ex withwaist wG = 203 (2) mm providing reduced longitudinal conﬁnement owing to its larger waist as compared tothe »70 mm crossed dipole beam waist. A transverse ‘tube trap’ along ez is provided by a blue-detuned(λLG=532 nm) Laguerre–Gauss (LG01) beam, tightly focused to a waist of w LG = 5.6 (5) mm . In our standard
conﬁguration these beams have powers PG=0.8(1) W and PLG=1.0(1) W, giving a peak transverse trapping
frequency w^ 2p = (wx wy )1 2 2p = 17 (2) kHz .
The transverse zero-point energy from ω⊥ produces an anti-conﬁning potential along ez due to thedivergence of the LG beam. The anti-trapping potential shown in green in ﬁgure 3 signiﬁcantly alters the overalllongitudinal potential
V (z ) =⎛ 2z 2 ⎞w (^0)+Vexp⎜ - 2 ⎟ - V0,t
1 + (z - z a )2 zR2⎝ wG ⎠(2)(0)
where w^
2p denotes the peak transverse trapping frequency; za is the center of the anti-trap; zR is the Rayleighrange of the LG beam; and V0 is an energetic offset chosen such that the minimum of the potential is zero. The
black shaded curve in ﬁgure 3 shows the combined anharmonic potential of the longitudinal trap (red curve) andthe anti-conﬁning potential. Small amplitude longitudinal dipole oscillations in the combined potential havefrequency ωz/2π=12.1(2) Hz.
Figure 4 depicts our four step loading scheme. (i) We ﬁrst ramp up the intensity of the 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 mm ,the »5 mm 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ﬁnal 1064 nm longitudinal trap in 250 ms. (iv) In the ﬁnal 250 ms we simultaneously increase the intensity of the
LG beam to its ﬁnal value while removing the crossed dipole trap.
These 250 ms ramps were chosen to be adiabatic with respect to all the conﬁning potentials. Monopole anddipole 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.

3New J. Phys. 20 (2018) 113032

F Salces-Carcoba et al

Figure 3. Trapping potential along ez with both contributions from the anti-trap (solid green curve) and red-detuned beam (solid red
curve). The shaded region marking the total trapping potential illustrates the uncertainty from the parameters entering into
equation (2). This includes the covariance matrix for the parameters of V(z) from our global Yang–Yang ﬁt discussed in section 4.

Figure 4. Adiabatic loading procedure. Each curve shows the intensity of a laser beam. The dashed–dotted blue curve depicts the
crossed dipole trap; the solid green curve denotes the LG tube trap; the dashed red curve marks the longitudinal trapping laser.

We control the temperature of the 1DBG by varying the temperature T3D of the initial 3D Bose gas. We tune T3D
by adjusting the depth of the crossed dipole trap, covering the range from T3D = 34–320 nK, with an observed BEC
c» 160 nK . We determine T3D with time-of-ﬂight measurements. The number of atoms in thetransition at T3D
1DBG increases with decreasing T3D due to the increasing 3D density as T3D 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 ﬁnal evaporation such thatthe crossed dipole potential is in a ﬁxed vertical position prior to loading the tube trap.

3. Imaging
We derive the density n(z) from in situ absorption images. Our imaging system has a resolution of 1.85 (5) mm
and magniﬁcation that maps one 5.6 mm sensor pixel to 0.84 (1) mm in the object plane. In preparation for
imaging, we apply a 20 ms repump pulse to transfer the 5 S1/2 ∣F = 1, mF = 0ñ atoms into the 5 S1/2 ∣F = 2ñhyperﬁne manifold. We then image [23] with a circularly polarized λp=780 nm probe beam resonant withthe 5 S1/2 ∣F = 2ñ to 5 P3/ 2 ∣F = 3ñ transition for 20 ms at an average intensity of I=2.5Isat, where Isat =1.67 mW cm-2 is the saturation intensity of the resonant atomic transition. An image of the probe beamfollowing absorption Ia, the probe without atoms present Ip, and a dark frame with no probe present Id, are eachrecorded on a charge-coupled device camera. From these images we obtain the absorbed fraction
4

New J. Phys. 20 (2018) 113032

F Salces-Carcoba et al

Table 1. Table summarizing the different parameters of V(z). The calibrated values and theiruncertainties were used as the central values and bounds for the ﬁt.
Parameter

Calibrated value

Value from ﬁt

Calibration method

w^ 2p
za
zR
Vt kB
wG
wz 2p V(z)dz

17 (4) kHz
0 (10) m m
185 (29) m m-1.17 (25) m K
203 (2) m m
12.13 (20) Hz

17 (2) kHz-7.670 (8) m m
185 (5) m m-1.37 (6) m K

Transverse expansion in TOF
Alignment precision
Intensity proﬁle of LG beam
Intensity proﬁle of Gaussian beam
Intensity proﬁle of Gaussian beam
Small amplitude dipole oscillations

8.19 (30) μm

f = (Ip - Ia ) (Ip - Id ). For each set of experimental parameters we repeat the experiment for ≈100
realizations.
Our image analysis is a multiple step process. We ﬁrst preprocess the raw images to correct for backgroundartifacts and improve the signal-to-noise ratio by a factor of ≈10. We then extract the linear densities using anabsorption 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 popt for each Ia as a linear sum of Ip from all realizations by minimizing the sum
squared difference between I popt and Ia away from the atoms [24, 25]. This reconstruction reduces fringing due tovibrational motion that occurred between acquiring Ia and Ip, along with shot noise present on each Ip. We use
similar techniques to remove a systematic difference in dark counts between Ia and Ip, as well as to account forstructured read-out noise. We then compute mean absorbed fractions f¯ over each set of experimental
parameters, and use uncentered principal component analysis (PCA) to further suppress shot noise. From f¯ anda detector calibrated [17, 26, 27] in units of Isat, 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 radialconﬁnement gives an extent of  mw^ » 110 nm ) is below the resolution of both our imaging system andthe optical scattering length 3l2p 2p 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.

4. Results
The results of our image processing are 1D density proﬁles n( j) (z ) conﬁned in the same trapping potential butwith one of 24 different initial conditions labeled by j. In the LDA we expect that these density proﬁles can result
from Yang–Yang thermodynamics. For each j, both the temperature T ( j) and the overall chemical potential m(0j)are in principle unknown because of the lack of suitable reservoirs. As a result we obtain these quantities from ﬁts
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 proﬁle as shown in ﬁgure 1(c) with just two freeparameters T ( j) and m(0j) . We constrain the ﬁt 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 asadditional parameters in our ﬁt 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 ﬁt. An
additional uncalibrated parameter dz accounts for a tiny displacement of the 1DBG relative to the center of V(z)for times following the loading protocol. The inclusion of dz leaves the main results unchanged and its value lies
within the relative alignment uncertainty of the trap centers. Different combinations of ﬁxed parameters have noqualitative effect on the results. The third column in table 1 shows the potential parameters derived from theYang–Yang ﬁt. We evaluate the goodness-of-ﬁt with the reduced chi-squared cn2 » 1.5. Lastly, the residuals of
the ﬁt show systematic variations in Ia, which are reﬂected by the gray pixels in ﬁgure 1(c).
Figure 5(a) shows the reduced density versus reduced chemical potential for two initial conditions, eachplotting different cuts in the EoS n(μ, T). The continuous curves in ﬁgure 5(a) represent the Yang–Yang modelwith the T and μ0 from our ﬁts. For small chemical potential these density proﬁles are well described by the EoSof a non-interacting Bose gas while for μ>0 they approach the predictions of GPE mean-ﬁeld 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 ﬁgure 3, ω⊥ isslightly off-centered, ultimately resulting in the observed behavior. The scattered white dots on ﬁgure 2
5

New J. Phys. 20 (2018) 113032

F Salces-Carcoba et al

Figure 5. Results from Yang–Yang ﬁt. (a) EoS for two different realizations (ﬁlled colored circles) along with Yang–Yang EoS (solid
color curves) and mean-ﬁeld prediction (solid gray curve). (b) Output parameters describing the state of the different T3D realizations(left side panels) and the subsequent time evolution of the lowest T3D realization (right side panels). The upper panels display the 1D
absolute temperature T while the bottom panels display the same temperature in units of the peak degeneracy temperature Td(0) .
(c) Total number of atoms as a function of time (blue diamonds) along with the ﬁt to the one and three-body loss model (solid dark
blue curve).

represent these two traces in the γ, T/Td 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 ﬁts in which we varied T3D, the initial 3D
temperature (left panels); or the hold time t in the 1D trap for our lowest T3D (right panels) cloud. In the top left
panel, we observe that as a function of decreasing T3D, the 1D temperature T ﬁrst remains constant and thencounterintuitively increases. In contrast, as shown in the bottom left panel, the degeneracy parameter deﬁned as
T Td(0), is a monotonically increasing function of T3D showing how the more degenerate 3D clouds result in
more degenerate 1DBGs.
For the lowest achievable T3D 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 ﬁgures 5(b) and (c), respectively), yet the 1DBG does not become
more degenerate (bottom right panel in ﬁgure 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 ≈700 nK. The inability ofsuch evaporative cooling to increase or even maintain degeneracy results from the slower relative decrease in Twith respect to Td 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-bodyloss and three-body loss from photon scattering, background gas collisions and inelastic three-bodycollisions [5].
Figure 5(c) shows the measured atom number N (blue diamonds). We ﬁt the decay model to the observed
number (dark blue curve), giving a one-body loss coefﬁcient K11D = 0.108 (2) s-1 as well as a three-body loss
coefﬁcient K31D = 4.36 (7) ´ 10-29 cm6 s-1. The value of K11D 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
coefﬁcient from our ﬁt is in excess of the intrinsic 3D three-body loss coefﬁcient K33D = 5.8 (3) ´ 10-30 cm6 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 sufﬁcient kinetic energy [29] overcome the longitudinal barrier of V(z).
6New J. Phys. 20 (2018) 113032

F Salces-Carcoba et al

5. Conclusions
We realized individually trapped 1DBGs in a crossed dipole trap formed by a blue-detuned LG01 beam and a reddetuned Gaussian beam. We benchmarked the EoS computed from Yang–Yang thermodynamics against themeasured density proﬁles. We found that evaporative cooling along the edges of the tube trap took placealthough 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.

Appendix A. Image preprocessing and extraction of modeled linear densities
We perform preprocessing of our raw images in order to improve signal to noise and correct for known
systematic error, before extracting linear densities from each set of experimental parameters using an absorption
model. The analysis pipeline from raw absorption, ﬂat ﬁeld and dark ﬁeld images to linear densities is as follows.
In the below descriptions of our analysis pipeline we alternately use bold symbols such as u when we are treating
an image array as a vector for the purposes of linear algebra, and ordinary symbols with spatial dependenceu(x, y) when we are treating images as functions of space.
A.1. Probe image reconstruction
For each shot we reconstruct an optimal probe image I popt as a linear sum of probe images from all shots:
I popt = Fc ,(A.1)

where I popt is I popt unwrapped into a column vector, F is a matrix containing all probe images as columns and c isa vector of coefﬁcients. The optimal coefﬁcients are determined by weighted linear regression:(FTWF ) c = FTWIa,(A.2)

where W is a diagonal matrix of weights equal to zero in a region of interest (ROI) about the atoms and one
otherwise, and Ia is the image (as a column vector) from the shot in question of the probe with atoms present.
The vector of coefﬁcients c is determined by numerically solving the linear system, leading to an I popt thatminimizes the sum squared error with Ia in the region outside the ROI. This probe reconstruction both reduces
fringing due to vibrational motion that occurs between exposures within a shot, and reduces shot noise present
on each reconstructed probe image on account of the dimensionality reduction entailed by linear
regression [24, 25].
A.2. Dark ﬁeld reconstruction
We correct for a small spatially inhomogeneous systematic difference in counts (≈1.2 max, ≈ 0.2 typical)
between absorption and probe images, which we attribute to variation in ambient brightness over the 60 Hz
mains power cycle (this is systematic rather than random, as each shot is synchronized to begin at the same pointin the 60 Hz cycle). We ﬁt a candidate two-dimensional function to the measured average difference between
absorption and probe images, which we include as a reference image Idsys in the above linear regression in order
to extract a coefﬁcient csys for each absorption image for how much of this offset was present (we obtain csys≈1
in all cases indicating little shot-to-shot variation in the offset).
We then compute the mean dark frame Id over all shots, and perform 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 ﬁve). We project each absorption image onto these
eigenvectors Ippca4 and Ippca5 of the probe images in order to determine coefﬁcients cpca4 and cpca5. A
reconstructed dark ﬁeld image Idrecon is then computed for each shot as
Idrecon = Id + c sysIdsys + c pca4I ppca4 + c pca5I ppca5.

7

(A.3)

New J. Phys. 20 (2018) 113032

F Salces-Carcoba et al

A.2.1. Absorbed fraction and saturation parameter.
then computed for each shot asf=

S=

Absorbed fraction f and saturation parameter S images are
I popt - Ia

I popt - Idrecon

,

(A.4)

1 opt
(I p - Idrecon) ,
Isat

(A.5)

where Isat is the saturation intensity in camera count units. The ‘naïve’ optical depth for a single shot can then becomputed as:
ODna¨ve= - a log (1 - f ) + Sf ,
ı

(A.6)

where a º s0 seff = Isat eff Isat is the empirically determined ratio between the ideal two-level and effective
scattering cross sections due to imperfect polarization and magnetic ﬁeld orientation. The average naïve opticaldepth over all repeated shots for each set of experimental parameters is computed as:
ODnaı̈ve » - a log (1 - f ) + SA ,(A.7)

where the mean product of absorbed fraction and f and saturation parameter S are taken over only the repeatedshots 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 themean 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 3l2p 2p 2 and our imaging
resolution, both of which are violations of the assumptions of the BL 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 BL law to account for this, presented further below.
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 theeffect of shot noise on column sums of the data. Since the point spread function resulting from diffraction in our
imaging system is ﬁxed from shot to shot, this has the effect of projecting the measured absorbed fractions onto theempirically 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 the75-pixel high ROI that entirely contains our imaging system’s point spread function to form fROI (x , y ). Then, foreach x position xi in the image, we extract the vertical lines of fROI (x , y ) from all sets of experiment parameters, atthat x pixel and the nearest four x pixels. Treating each vertical line as a vector4 fROI (x i ) = fROI (x i, y ), we obtainthe set of vectors { fROI (x i + j ), -2  j  2} and perform uncentered PCA [33], keeping only the ﬁrst four
normalized eigenvectors {vˆn (x i ), n Î [1, 4]}. We then project the (also cropped to the ROI) vertical lines fROI (x i )of the absorbed fractions for each individual shot at the original x position onto the subspace spanned by these
vectors:
4

fred (x i , y ) = fred (x i ) = å [vˆn (x i ) · fROI (x i )] vˆn (x i ).

(A.8)

n= 1

Thus we dimensionally reduce each vertical slice (within the ROI) of each shot’s absorbed fraction onto a basis of fourbasis functions chosen by uncentered PCA of the vertical slices of all mean absorbed fractions at that x pixel and thenearest four other x pixels. We then compute mean, dimensionally reduced absorbed fractions fred (x , y ) for each setof 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 asDfred (x , y ) =4s [{ fred (x , y )}]
N

,

(A.9)

The mixed vector versus 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.

8

New J. Phys. 20 (2018) 113032

F Salces-Carcoba et al

where σ is the standard deviation over all repeated shots for one set of experiment parameters, {fred(x, y)} is theset of all dimensionally reduced absorbed fractions for those shots, and N is the number of repeated shots for thatset of experiment parameters.
A.3.1. Saturation parameter at the position of the atoms. Due to the point spread function of our imaging systembeing 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 theatoms, which is where saturation effects are relevant. We estimate a saturation parameter S0 (x ) for each set ofexperiment parameters at the estimated y position of the atoms by interpolating the mean saturation parameterfor 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 ﬁty0 (x ) = ax 2 + bx + c with parameters determined by maximizing the total absorption over all shots:⎛⎞[a , b , c ] = argmax ⎜⎜ å å interp [S (x i , y )](ax i2 + bx i + c ) ⎟⎟ ,ya, b, c ⎝shots x i⎠(A.10)

where interp is a one-dimensional spline interpolation function interpolating in the y direction only. Theyestimated saturation parameter at each x position for each set of experiment parameters is thenS 0 (x ) = interp [S (x, y )]( y0 (x )).

(A.11)yA.4. Naïve linear density
We can now compute an improved naïve optical depth for each set of experiment parameters using thedimensionally reduced absorbed fractions and interpolated saturation parameter as:
ODnaı̈ve red (x , y ) » - a log (1 - fred (x , y )) + S 0 (x ) fred (x , y ) ,(A.12)and then compute a naïve linear density nnaı̈ve (x ) at each x position by dividing by the cross section andintegrating along y:nnaı̈ve (x ) =Dyå ODnaive red (x, y ) ,s0 y(A.13)

where Δy is the pixel size.
A.5. Modeled linear density
We face three related problems in computing the column density ncol(x, y) given a measured absorbed fraction
f (x, y) and saturation parameter S(x, y) via the solution to the BL law [27]:s0 n col (x , y ) = - a log (1 - f (x , y )) + S (x , y ) f (x , y ).

(A.14)

The ﬁrst 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 asurrounding region of space with cross sectional area given by the absorption cross section σ0 [34]. The ﬁnal
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 ncol(x, y) from the usual solution to the BL law,
we can only determine the convolution of ncol(x, y) with some absorption proﬁle g(y) that takes into accountboth the ﬁnite absorption region and the diffusion of atoms from their initial positions in the y direction, thedirection in which g(y) is not small compared to our atom cloud. With this in mind, the solution to the BL law
can be modiﬁed to reads0 (n col * g )(x , y , t ) = - a log (1 - f (x , y , t )) + S (x , y ) f (x , y , t ) ,(A.15)

where the convolution is only along the y direction. Atomic diffusion and diffraction imply that the only quantity
we have experimental access to isfmeas (x ) =1tt

ò0 dt ò dy f (x, y, t ) = åy fred (x, y ) Dy,(A.16)

that is, we only observe a time average of absorption over the imaging pulse time τ, and we only observe theintegral of the undiffracted absorbed fraction, since diffraction preserves this integral.
If the second term of the solution to the BL 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 termthat causes a problem, since its integral is not conserved under diffraction.
9

New J. Phys. 20 (2018) 113032

F Salces-Carcoba et al

Given a model for (ncol * g )(x , y , t ) with a single parameter n(x) for the linear density at each x position anda measurement fmeas(x) for each x position, we can invert (A.16) and (A.15) to obtain the linear density, under theassumptions of the model.
Our model is the following: The absorption proﬁle g(y) is approximated by a Gaussian with unit integral andstandard deviation σy(t) equal at t=0 to the optical scattering length s0 p = 3l2p 2p 2 and increasing due
to atomic diffusion as time elapses. Since the atom cloud is narrower than this absorption proﬁle, we
approximate the convolution (n * g )(x , y , t ) as:
⎤
⎡
y2
⎥,
exp ⎢ - 2⎢⎣ 2s y (x , t ) ⎥⎦2ps 2y (x , t )n1D (x )(n * g )(x , y , t ) »(A.17)

where s 2y (x , t ) is increasing due to momentum diffusion:s0
1+ s 2vy (x , t ) t 2,
(A.18)
p3
where the mean squared y velocity s 2vy is given by the scattering rate Rscat and recoil velocity vrec:
12
s 2vy (x , t ) = (2p )-1R scat (x ) vrect,
(A.19)
3
which is approximating isotropic scattering such that the per-scattering-event expected squared change in2
velocity is vrec
3. The scattering rate, ignoring Doppler shifts away from resonance, is, in terms of the saturation
parameter S:
s 2y (x , t ) =R scat (x ) =G S (x ).
2 1 + S (x )(A.20)

Putting this all together, the modeled y variance of the absorption proﬁle is:s 2y (x , t ) =s0
S (x )
G
2 3
+
vrec
t .
p
36p 1 + S (x )(A.21)

Over the duration of our imaging pulse, the diffusion described by the second term results in an increase in theabsorption proﬁle’s standard deviation by ≈30% compared to the effect of the non-zero optical scattering lengthalone.
Using our absorption model (A.17) with an absorption proﬁle with y variance given by (A.21), and saturation
parameter S(x, y) given by our estimate S0(x), our modiﬁed BL law solution (A.15) becomes:
s0 n (x )⎡⎤y2exp ⎢ - 2s 2 (x , t ) ⎥⎣⎦y
2ps 2y (x , t )= - a log (1 - f (x , y , t )) + S 0 (x ) f (x , y , t ) ,(A.22)

which, if numerically inverted, deﬁnes 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 fmeas(x) for that linear density. Numerically inverting this function then yields our
ﬁnal aim, of a function that takes as input fmeas(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 non-trivial calculation to extract modeled linear densities from ourdimensionally-reduced mean absorbed fractions and interpolated saturation parameters for each set ofexperiment 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 theabsorption model.

Appendix B. Yang–Yang thermodynamics
We use the Yang–Yang model [8] to describe our data. The exact diagonalization of the underlying Hamiltonianis carried out with the use of the thermodynamic Bethe ansatz (TBA) (T>0 Bethe ansatz). From the TBA thefollowing set of ﬁrst-order integral equations can be derived (k ) = 2k 2k T ¥2c-m- B òln (1 + e- (q ) kB T ) dq ,2m2p -¥ c 2 + (k - q)2

2pf (k)(1 + e (k) kB T ) = 1 +

10

¥

2c

ò-¥ c 2 + (k - q)2 f (q) dq,

(B.1)
(B.2)

New J. Phys. 20 (2018) 113032

F Salces-Carcoba et aln=¥ò-¥ f (q) dq,(B.3)

where m is the mass, kB is the Boltzmann constant, T is the temperature, μ the chemical potential, andc=mg/ÿ2 is the interaction wavenumber. We use g = 2w^ a3D (1 - Ca3D l^), where a3D is the threedimensional s-wave scattering length, C is a constant of order unity [35], and l^ =  mw^ is the oscillator
length. Both k and q label momenta. The above equations can be solved recursively to compute n, the linear
density, given the values for μ and T, the chemical potential and temperature.
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 l⊥, a3D and m. We recursively solve forò(k) and f (k) from which we ultimately compute the density n. We transform all the momentum and energy
quantitiesk˜ = k

2mkB T  2 ,(B.4)

E˜ = E kB T,

(B.5)

1
c˜
ln (1 + e-˜ (q˜)) dq˜ ,2
p c˜ + (k˜ - q˜)2

(B.6)

so that the ﬁrst two YY equations read2˜ (k˜ ) = k˜ - m˜ - ò

¥

-¥

˜
2pf (k˜ )(1 + e˜ (k )) = 1 +

¥

1

c˜

ò-¥ p c˜ 2 + (k˜ - q˜)2 f (q˜) dq˜.

(B.7)

We denote the Lieb–Liniger kernel (a normalized Lorentzian) as L(c, q). Our numerical solver performs a kspace convolution using the scipy.signal.fftconvolve method to evaluate the integrals. For this we
use a Nk=1024 points grid that covers the range k = [-10 k th, 10 k th ], where k th = 2mkB T  2 is thethermal wavenumber. We have explicitly veriﬁed that for the parameters in our experiment using different grids
give no changes to the predicted EoS. We initialize ò0(k)=k2−μ 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 Nk <  tol issatisﬁed, we solve for f (k) with an initial guess f0 (k ) = [2p (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.

References[1] Bloch I, Dalibard J and Zwerger W 2008 Many-body physics with ultracold gases Rev. Mod. Phys. 80 885–964[2] Chin C, Grimm R, Julienne P and Tiesinga E 2010 Feshbach resonances in ultracold gases Rev. Mod. Phys. 82 1225–86[3] Lahaye T, Menotti C, Santos L, Lewenstein M and Pfau T 2009 The physics of dipolar bosonic quantum gases Rep. Prog. Phys. 72 126401[4] Paredes B, Widera A, Murg V, Mandel O, Folling S, Cirac I, Shlyapnikov G V, Hansch T W and Bloch I 2004 Tonks–Girardeau gas ofultracold atoms in an optical lattice Nature 429 277–81[5] Laburthe Tolra B, O’Hara K M, Huckans J H, Phillips W D, Rolston S L and Porto J V 2004 Observation of reduced three-bodyrecombination in a correlated 1d degenerate Bose gas Phys. Rev. Lett. 92 190401[6] Kinoshita T, Wenger T and Weiss D S 2004 Observation of a one-dimensional Tonks–Girardeau gas Science 305 1125–8[7] Lieb E H and Liniger W 1963 Exact analysis of an interacting Bose gas: I. The general solution and the ground state Phys. Rev. 130
1605–16[8] Yang C N and Yang C P 1969 Thermodynamics of a one-dimensional system of bosons with repulsive deltafunction interaction
J. Math. Phys. 10 1115–22[9] van Amerongen A H, van Es J J P, Wicke P, Kheruntsyan K V and van Druten N J 2008 Yang–Yang thermodynamics on an atom chipPhys. Rev. Lett. 100 090402[10] Vogler A, Labouvie R, Stubenrauch F, Barontini G, Guarrera V and Ott H 2013 Thermodynamics of strongly correlated onedimensional Bose gases Phys. Rev. A 88 031603[11] Yang B, Chen Y-Y, Zheng Y-G, Sun H, Dai H-N, Guan X-W, Yuan Z-S and Pan J-W 2017 Quantum criticality and the Tomonaga–
Luttinger liquid in one-dimensional Bose gases Phys. Rev. Lett. 119 165701[12] Rauer B, Grišins P, Mazets I E, Schweigler T, Rohringer W, Geiger R, Langen T and Schmiedmayer J 2016 Cooling of a one-dimensional
Bose gas Phys. Rev. Lett. 116 030402[13] Bongs K, Burger S, Dettmer S, Hellweg D, Arlt J, Ertmer W and Sengstock K 2001 Waveguide for Bose–Einstein condensates Phys. Rev.
A 63 031602[14] Kheruntsyan K V, Gangardt D M, Drummond P D and Shlyapnikov G V 2005 Finite-temperature correlations and density proﬁles ofan inhomogeneous interacting one-dimensional Bose gas Phys. Rev. A 71 053615[15] Armijo J, Jacqmin T, Kheruntsyan K and Bouchoule I 2011 Mapping out the quasicondensate transition through the dimensionalcrossover from one to three dimensions Phys. Rev. A 83 021605[16] Kohn W and Sham L J 1965 Self-consistent equations including exchange and correlation effects Phys. Rev. 140 A1133–8

11New J. Phys. 20 (2018) 113032

F Salces-Carcoba et al[17] Yefsah T, Desbuquois R, Chomaz L, Günter K J and Dalibard J 2011 Exploring the thermodynamics of a two-dimensional Bose gas
Phys. Rev. Lett. 107 130401[18] Ku M J H, Sommer A T, Cheuk L W and Zwierlein M W 2012 Revealing the superﬂuid lambda transition in the universal
thermodynamics of a unitary Fermi gas Science 335 563–7[19] Mukherjee B, Yan Z, Patel P B, Hadzibabic Z, Yefsah T, Struck J and Zwierlein M W 2017 Homogeneous atomic Fermi gases Phys. Rev.
Lett. 118 123401[20] Greiner M, Bloch I, Hänsch T W and Esslinger T 2001 Magnetic transport of trapped cold atoms over a large distance Phys. Rev. A 63
031401[21] Lin Y-J, Perry A R, Compton R L, Spielman I B and Porto J V 2009 Rapid production of 87Rb Bose–Einstein condensates in a combined
magnetic and optical potential Phys. Rev. A 79 063631[22] Starkey P T, Billington C J, Johnstone S P, Jasperse M, Helmerson K, Turner L D and Anderson R P 2013 A scripted control system forautonomous hardware-timed experiments Rev. Sci. Instrum. 84 085111[23] Genkina D, Aycock L M, Stuhl B K, Lu H-I, Williams R A and Spielman I B 2016 Feshbach enhanced s-wave scattering of fermions:direct observation with optimized absorption imaging New J. Phys. 18 013001[24] Ockeloen C F, Tauschinsky A F, Spreeuw R J C and Whitlock S 2010 Detection of small atom numbers through image processing Phys.
Rev. A 82 061606[25] Li X, Ke M, Yan B and Wang Y 2007 Reduction of interference fringes in absorption imaging of cold atom cloud using eigenface method
Chin. Opt. Lett. 5 128–30[26] Hueck K, Luick N, Sobirey L, Siegl J, Lompe T, Moritz H, Clark L W and Chin C 2017 Calibrating high intensity absorption imaging of
ultracold atoms Opt. Express 25 8670–9[27] Reinaudi G, Lahaye T, Wang Z and Guéry-Odelin D 2007 Strong saturation absorption imaging of dense clouds of ultracold atoms Opt.
Lett. 32 3143–5[28] Burt E A, Ghrist R W, Myatt C J, Holland M J, Cornell E A and Wieman C E 1997 Coherence, correlations, and collisions: what one
learns about Bose–Einstein condensates from their decay Phys. Rev. Lett. 79 337–40[29] Fedichev P O, Reynolds M W and Shlyapnikov G V 1996 Three-body recombination of ultracold atoms to a weakly bound s level Phys.
Rev. Lett. 77 2921–4[30] Yang L, Guan L and Pu H 2015 Strongly interacting quantum gases in one-dimensional traps Phys. Rev. A 91 043634[31] Yang L and Pu H 2016 Bose–Fermi mapping and a multibranch spin-chain model for strongly interacting quantum gases in one
dimension: dynamics and collective excitations Phys. Rev. A 94 033614[32] Cole W S, Lee J, Mahmud K W, Alavirad Y, Spielman I B and Sau J D 2017 Spin–orbit coupled bosons in one dimension: emergent
gauge ﬁeld and Lifshitz transition arXiv:1711.05794[33] Cadima J and Jolliffe I T 2009 On relationships between uncentred and column-centred principal component analysis Pak. J. Stat. 25
473–503[34] Striebel M, Wrachtrup J and Gerhardt I 2017 Absorption and extinction cross sections and photon streamlines in the optical near-ﬁeld
Sci. Rep. 7 15420[35] Dunjko V, Lorent V and Olshanii M 2001 Bosons in cigar-shaped traps: Thomas–Fermi regime, Tonks–Girardeau regime, and in
between Phys. Rev. Lett. 86 5413–6

12