Emergence of coherence in a uniform quasi-two-dimensional Bose gas

Phase transitions are ubiquitous in our three-dimensional world. By contrast most conventional transitions do not occur in infinite uniform two-dimensional systems because of the increased role of thermal fluctuations. Here we explore the dimensional crossover of Bose-Einstein condensation (BEC) for a weakly interacting atomic gas confined in a novel quasi-two-dimensional geometry, with a flat in-plane trap bottom. We detect the onset of an extended phase coherence, using velocity distribution measurements and matter-wave interferometry. We relate this coherence to the transverse condensation phenomenon, in which a significant fraction of atoms accumulate in the ground state of the motion perpendicular to the atom plane. We also investigate the dynamical aspects of the transition through the detection of topological defects that are nucleated in a quench cooling of the gas, and we compare our results to the predictions of the Kibble-Zurek theory for the conventional BEC second-order phase transition.

case ζ 1 (with ζ = k B T /hν z ), the transverse condensation phenomenon induces an extended in-plane coherence. Second, we explore the quench dynamics of the gas prepared in an initial state such that ζ 1 and observe density holes associated to vortices. We study the relation between the cooling rate and the number of vortices that subsist after a given relaxation time, and we compare our results with the predictions of the Kibble-Zurek theory.

Results
Production of uniform gases in quasi-2D geometries. We prepare a cold 3D gas of rubidium ( 87 Rb) atoms using standard laser and evaporative cooling techniques. Then we transfer the gas in a trap formed with two orthogonal laser beams at wavelength 532 nm, shorter than the atomic resonance wavelength (780 nm), so that the atoms are attracted towards the regions of low light intensity (Fig. 1 a). The strong confinement along the z direction (vertical) is provided by a laser beam propagating along the x direction. It is prepared in a Hermite-Gauss mode, with a node in the plane z = 0, and provides a harmonic confinement along the z direction with a frequency ν z in the range 350−1500 Hz. For the confinement in the horizontal xy plane, we realize a box-like potential by placing an intensity mask on the second laser beam path, propagating along the z direction ( Fig. 1 b). Depending on the study to be performed, we can vary the shape (disk, square, double rectangle) and the area A (from 200 to 900 µm 2 ) of the region accessible to the gas in the plane. The relevance of our system for the study of 2D physics is ensured by the fact that the size of the ground state along the z direction a z = h/(mν z )/2π ∼ 0.3 -0.6 µm is very small compared to the in-plane extension √ A ∼ 15 -30 µm. The number of atoms N that can be stored and reliably detected in this trap ranges between 1000 and 100 000. We adjust the temperature of the gas in the interval T ∼ 10 − 250 nK by varying the intensity of the beam creating the box potential, taking advantage of evaporative cooling on the edges of this box. The ranges spanned by ν z and T allow us to explore the dimensional crossover between the thermally frozen regime (ζ 1) and the unfrozen one (ζ 1). Examples of in situ images of 2D gases are shown in Fig. 1 c,d,e.
Phase coherence in 2D geometries. For an ideal gas an important consequence of Bose-Einstein statistics is to increase the range of phase coherence with respect to the prediction of Boltzmann statistics. Here coherence is characterized by the one-body correlation function G 1 (r) = ψ † (r)ψ(0) , whereψ(r) [resp.ψ † (r)] annihilates (resp. creates) a particle in r, and where the average is taken over the equilibrium state at temperature T . For a gas of particles of mass m described by Boltzmann statistics, G 1 (r) is a Gaussian function ∝ exp(−πr 2 /λ 2 T ), where λ T = h/(2πmk B T ) 1/2 is the thermal wavelength.
Consider the particular case of a 2D Bose gas (e.g., ζ 1). When its phasespace-density D ≡ ρλ 2 T becomes significantly larger than 1 (ρ stands for the 2D spatial density), the structure of G 1 (r) changes. In addition to the Gaussian function We slice a horizontal sheet from a 3D cold gas of 87 Rb atoms using a blue detuned laser beam propagating along x and shaped with an intensity node in the z = 0 plane. It creates an adjustable harmonic confinement along z of frequency ν z = 350-1500 Hz. We superimpose a hollow beam propagating along z and producing a uniform confinement in the xy plane (see b). The power P box of this beam is ramped up in 10 ms to its maximal value P max box corresponding to a potential barrier U box ∼ k B ×3 µK for the 87 Rb atoms. After holding P box at P max box for 0.5 s, we lower it linearly to its final value P f box in a typical time of t evap = 2 s and keep it constant for a typical t hold = 0.5 s. We vary P f box to adjust the final temperature of the gas via evaporative cooling. (b) The in-plane (xy) confinement is provided by a blue-detuned laser beam shaped by placing a dark intensity mask on its path and imaging it at the position of the atoms. (c, d, e): In-situ density distributions of uniform gases trapped in a disk of radius R = 12 µm, a square box of length L = 30 µm, and two coplanar and parallel rectangular boxes of size 24 × 12 µm 2 , spaced by d = 4.5 µm. These distributions are imaged using a high intensity absorption imaging technique (see methods). mentioned above, a broader feature ∝ exp(−r/ ) develops, with the characteristic length that increases exponentially with D (see [2] and supplementary material) Usually two main effects amend this simple picture: • In a finite system, when the predicted value of becomes comparable to the size L of the gas, one recovers a standard Bose-Einstein condensate, with a macroscopic occupation of the ground state of the box potential [38]. The G 1 function then takes non-zero values for any r ≤ L and the phase coherence extends over the whole area of the gas.
• In the presence of weak repulsive interactions, the increase of the range of G 1 for D 1 is accompanied with a reduction of density fluctuations, with the formation of a "quasi-condensate" or "pre-superfluid" state [15,16,39]. This state is a medium that can support vortices, which will eventually pair at the superfluid BKT transition for a larger phase-space-density, around D ∼ 8 − 10 for the present strength of interactions [39]. At the transition point, the coherence length diverges and above this point, G 1 (r) decays algebraically.
Role of the third dimension for in-plane phase coherence. When the thermal energy k B T is not negligibly small compared to the energy quantum hν z for the tightly confined dimension, the dynamics associated to this direction brings interesting novel features to the in-plane coherence. First, we note that the function G 1 (r) can be written in this case as a sum of contributions of the various states j z of the z motion (see supplementary material). The term with the longest range corresponds to the ground state j z = 0, with an expression similar to (1) where D is replaced by the phase-space-density D 0 associated to this state. Now, consider more specifically the unfrozen regime ζ 1. In this case one expects that for very dilute gases only a small fraction f 0 of the atoms occupies the j z = 0 state; Boltzmann statistics indeed leads to f 0 = 1 − e −1/ζ ≈ 1/ζ 1. However for large total phase-space densities D tot , Bose-Einstein statistics modifies this result through the transverse condensation phenomenon (BEC ⊥ ) [1]: The phase-space-density that can be stored in the excited states j z = 0 is bounded from above, and D 0 can thus become comparable to D tot . This large value of D 0 leads to a fast increase of the corresponding range of G 1 (r), thus linking the transverse condensation to an extended coherence in the xy plane. This effect plays a central role in our experimental investigation.
Phase coherence revealed by velocity distribution measurements. To characterize the coherence of the gas, we study the velocity distribution, i.e., the Fourier transform of the G 1 (r) function. We approach this velocity distribution in the xy plane by performing a 3D time-of-flight (3D ToF): We suddenly switch off the trapping potentials along the three directions of space, let the gas expand for a duration (a-f ) Surface density distribution ρ(x, y) (first row) and corresponding radial distributions (green symbols) obtained by azimuthal average (second row). The distribution is measured after a 12 ms time-of-flight for a gas initially confined in a square of size L = 24 µm, with a trapping frequency ν z = 365 Hz along the z direction. The temperatures T and atom numbers N for these three realizations are a,d: (155 nK, 28 000), b,e: (155 nK, 38 000), c,f: (31 nK, 19 000). The continuous red lines are fits to the data by a function consisting in the sum of two Gaussians corresponding to N 1 and N 2 atoms (N = N 1 +N 2 ). The Gaussian of largest width (N 2 atoms) is plotted as a blue dashed line. The bimodal parameter ∆ = N 1 /N equals a,d: 0.01, b,e: 0.12 and c,f: 0.60. (g) Variation of ∆ with N for a gas in the same initial trapping configuration as a-f and for T = 155 nK (red symbols). Error bars are the standard errors of the mean of the binned data set (with 4 images per point on average). The solid line is a fit to the data by the function f (N ) = (1 − (N c /N ) 0.6 ) for N > N c , and f (N ) = 0 for N ≤ N c , from which we deduce N c (T ). Here N c = 3.2 (1) × 10 4 , where the uncertainty range is obtained by a jackknife resampling method, i.e. fitting samples corresponding to a randomly chosen fraction of the global data set. τ , and finally image the gas along the z axis. In such a 3D ToF, the gas first expands very fast along the initially strongly confined direction z. Thanks to this fast density drop, the interparticle interactions play nearly no role during the ToF and the slower evolution in the xy plane is governed essentially by the initial velocity distribution of the atoms. The time-of-flight (ToF) duration τ is chosen so that the size expected for a Boltzmann distribution τ k B T /m is at least twice the initial extent of the cloud. Typical examples of ToF images are given in Fig. 2 a-f. Whereas for the hottest and less dense configurations, the spatial distribution after ToF has a quasi-pure Gaussian shape, a clear non-Gaussian structure appears for larger N or smaller T . A sharp peak emerges at the center of the cloud of the ToF picture, signaling an increased occupation of the low-momentum states with respect to Boltzmann statistics, or equivalently a coherence length significantly larger than λ T .
In order to analyze this velocity distribution, we chose as a fit function the sum of two Gaussians of independent sizes and amplitudes, containing N 1 and N 2 atoms, respectively (see Fig. 2 d-f ). We consider the bimodality parameter ∆ = N 1 /N defined as the ratio of the number of atoms N 1 in the sharpest Gaussian to the total atom number N = N 1 + N 2 . A typical example for the variations of ∆ with N at a given temperature is shown in Fig. 2 g for an initial gas with a square shape (side length L = 24 µm). It shows a a sharp crossover, with essentially no bimodality (∆ 1) below a critical atom number N c (T ) and a fast increase of ∆ for N > N c (T ). We extract the value N c (T ) by fitting the function ∆ ∝ (1 − (N c /N ) 0.6 ) to the data. We chose this function as it provides a good representation of the predictions for an ideal Bose gas in similar conditions (see methods).
Phase coherence revealed by matter-wave interference. Matter-wave interferences between independent atomic or molecular clouds is a powerful tool to monitor the emergence of extended coherence [4,14,[40][41][42]. To observe these interferences in our uniform setup, we first produced two independent gases of similar density and temperature confined in two coplanar parallel rectangles, separated by a distance of 4.5 µm along the x direction (see Fig. 1 e). Then we suddenly released the box potential providing confinement in the xy plane, while keeping the confinement along the z direction (2D ToF). The latter point ensures that the atoms stay in focus with our imaging system, which allows us to observe interference fringes with a good resolution in the region where the two clouds overlap. A typical interference pattern is shown in Fig. 3 a, where the fringes are (roughly) parallel to the y axis, and show some waviness that is linked to the initial phase fluctuations of the two interfering clouds.
We use these interference patterns to characterize quantitatively the level of coherence of the gases initially confined in the rectangles. For each line y of the pixelized image acquired on the CCD camera, we compute the x-Fourier transform ρ(k, y) of the spatial density ρ(x, y) (Fig. 3 b). For a given y this function is peaked at a momentum k p (y) > 0 that may depend (weakly) on the line index y. Then we consider the function that characterizes the correlation of the complex fringe contrastρ[k p (y), y] along two lines separated by a distance d Here * denotes the complex conjugation and the average is taken over the lines y that overlap with the initial rectangles. If the initial clouds were two infinite, parallel lines with the same G 1 (y), one would have γ(d) = |G 1 (d)| 2 [43]. Here the non-zero extension of the rectangles along x and their finite initial size along y make it more difficult to provide an analytic relation between γ and the initial G 1 (r) of the gases. However γ(d) remains a useful and quantitative tool to characterize the fringe pattern. For a gas described by Boltzmann statistics, the width at 1/e of G 1 (r) is λ T / √ π and remains below 1 µm for the temperature range investigated in this work. Since we are interested in the emergence of coherence over a scale that significantly exceeds this value, we use the following average as a diagnosis tool Γ = γ(d) , average taken over the range 2 µm < d < 5 µm.
For the parameter Γ to take a value significantly different from 0, one needs a relatively large contrast on each line, and relatively straight fringes over the relevant distances d, so that the phases of the different complex contrasts do not average out. For a given temperature T , the variation of Γ with N shows the same thresholdtype behaviour as the bimodality parameter ∆. One example is given in Fig. 3 c, from which we infer the threshold value for the atom number N c (T ) needed to observe interference fringes with a significant contrast.
Scaling laws for the emergence of coherence. We have plotted in Fig. 4 the ensemble of our results for the threshold value of the total 2D phase-space-density D tot,c ≡ N c λ 2 T /A as a function of ζ = k B T /hν z , determined both from the onset of bimodality as in Fig. 2 g (closed symbols) or from the onset of visible interference as in Fig. 3 c (open symbols). Two trapping configurations have been used along the z direction, ν z = 1460 Hz and ν z = 365 Hz. In the first case, the z direction is nearly frozen for the temperatures studied here (ζ 2). In the second one, the z direction is thermally unfrozen (ζ 8). All points approximately fall on a common curve, independent of the shape and the size of the gas: D tot, c varies approximately linearly with ζ with the fitted slope 1.4 (3) for ζ 8 and approaches a finite value ∼ 4 for ζ 2.
In the frozen case, a majority of atoms occupy the vibrational ground state j z = 0 of the motion along the z direction, so that D tot essentially represents the 2D phase-space-density associated to this single transverse quantum state. Then for D tot ≥ 1, we know from Eq. (1) and the associated discussion that a broad component arises in G 1 with a characteristic length that increases exponentially with the phase-space-density. The observed onset of extended coherence around D tot ∼ 4 can be understood as the place where starts to exceed significantly λ T . The regime around D tot ∼ 4 is reminiscent of the presuperfluid state identified in [15,16]. It is different from the truly superfluid phase, which is expected at a higher phase-space-density (D tot ∼ 8) for our parameters [39]. Therefore the threshold violet open diamonds. Error bars show the 95% confidence bounds on the N c parameter of the threshold fits to the data sets. The black solid line shows a linear fit to the data for ζ > 8, leading to D tot, c = 1.4 (3) ζ. The black dash-dotted lines show contours of identical ratios of the coherence range to the thermal wavelength λ T . The coherence range is evaluated by the value of r at which G 1 (r) = G 1 (0)/20 (see text) and we plot (in increasing D tot order) ratios equal to 1, 1.2, 1.5, 2, 3 and 8. Boltzmann prediction corresponds to a ratio of ∼ 0.98. D tot, c is not associated to a true phase transition, but to a crossover where the spatial coherence of the gas increases rapidly with the control parameter N .
For ν z = 365 Hz, the gas is in the "unfrozen regime" (ζ 1), which could be naively thought as irrelevant for 2D physics since according to Boltzmann statistics, many vibrational states along z should be significantly populated. However thanks to the BEC ⊥ phenomenon presented above, a macroscopic fraction of the atoms can accumulate in the j z = 0 state. This happens when the total phase-space-density exceeds the threshold for BEC ⊥ (cf. supplementary material): In the limit ζ → ∞, BEC ⊥ corresponds to a phase transition of the same nature as the ideal gas BEC in 3D. In the present context of our work, we emphasize that although BEC ⊥ originates from the saturation of the occupation of the excited states along z, it also affects the coherence properties of the gas in the xy plane.
In particular when D tot rises from 0 to D tot, c , the coherence length in xy increases from ∼ λ T (the non-degenerate result) to ∼ a z , the size of the ground state of the z motion. This increase can be interpreted by noting that when BEC ⊥ occurs (Eq. 4), the 3D spatial density in the central plane (z = 0) is equal to g 3/2 (1)/λ 3 T , where g s is the polylogarithm of order s and g 3/2 (1) ≈ 2.612. For an infinite uniform 3D Bose gas with this density, a true Bose-Einstein condensation occurs and the coherence length diverges. Because of the confinement along the z direction, such a divergence cannot occur in the present quasi-2D case. Instead, the coherence length along z is by essence limited to the size a z of the j z = 0 state. When D tot = D tot, c the same limitation applies in the transverse plane, giving rise to coherence volumes that are grossly speaking isotropic. When D tot is increased further, the coherence length in the xy plane increases, while remaining limited to a z along the z direction. The results shown in Fig. 4 are in line with this reasoning. For ζ 1, the emergence of coherence in the xy plane occurs for a total phase-space-density D tot, c ∝ ζ, with a proportionality coefficient α = 1.4 (3) in good agreement with the prediction π 2 /6 ≈ 1.6 of Eq. (4).
We have also plotted in Fig. 4 contour lines characterizing the coherence range in terms of ζ and D tot . Using ideal Bose gas theory, we calculated the one-body coherence function G 1 (r) and determined the distance r f over which it decreases by a given factor f with respect to G 1 (0). We choose the value f = 20 to explore the long tail that develops in G 1 when phase coherence emerges. The contour lines shown in Fig. 4 correspond to given values of r 20 /λ T ; they should not be considered as fits to the data, but as an indication of a coherence significantly larger than the one obtained from Boltzmann statistics (for which r 20 ≈ λ T ). The fact that the threshold phase-space densities D tot, c follow quite accurately these contour lines validates the choice of tools (non Gaussian velocity distributions, matter-wave interferences) to characterize the onset of coherence.
Observation of topological defects. From now on we use the weak trap along z (ν z = 365 Hz) so that the onset of extended coherence is obtained thanks to the transverse condensation phenomenon. We are interested in the regime of strongly degenerate, interacting gases, which is obtained by pushing the evaporation down to a point where the residual thermal energy k B T becomes lower than the chemical potential µ (see methods for the calculation of µ in this regime). The final box potential is ∼ k B × 40 nK, leading to an estimated temperature of ∼ 10 nK, whereas the final density (∼ 50 µm −2 ) leads to µ ≈ k B × 14 nK. In these conditions, for most realizations of the experiment, defects are present in the gas. They appear as randomly located density holes after a short 3D ToF (Fig. 5 a,b), with a number fluctuating between 0 and 5. To identify the nature of these defects, we have performed a statistical analysis of their size and contrast, as a function of their location and of the ToF duration τ (Fig 5 c,d). For a given τ , all observed holes have similar sizes and contrasts. The core size increases approximately linearly with τ , with a nearly 100 % contrast. This favors the interpretation of these density holes as single vortices, for which the 2π phase winding around the core provides a topological protection during the ToF. This would be the case neither for vortex-antivortex pairs nor phonons, for which one would expect large fluctuations in the defect sizes and lower contrasts.
Dynamical origin of the topological defects. In principle the vortices observed in the gas could be due to steady-state thermal fluctuations. BKT theory indeed predicts that vortices should be present in an interacting 2D Bose gas around the superfluid transition point [11]. Such "thermal" vortices have been observed in nonhomogeneous atomic gases, either interferometrically [14] or as density holes in the trap region corresponding to the critical region [20]. However, for the large and uniform phase-space densities that we obtain at the end of the cooling process (ρλ 2 T ≥ 100), Ref. [44] predicts a vanishingly small probability of occurrence for such thermal excitations. This supports a dynamical origin for the observed defects.
To investigate further this interpretation, we can vary the two times that characterize the evolution of the gas, the duration of evaporation t evap and the hold duration after evaporation t hold (see Fig. 1 a). For the results presented in this section, we fixed t hold = 500 ms and studied the evolution of the average vortex number N v as a function of t evap . The corresponding data, given in Fig. 6 a, show a decrease of N v with t evap , passing from N v ≈ 1 for t evap = 50 ms to N v ≈ 0.3 for t evap = 250 ms. For longer evaporation times, N v remains approximately constant around 0.35 (5).
The decrease of N v with t evap suggests that the observed vortices are nucleated via a Kibble-Zurek (KZ) type mechanism [22,23,45], occurring when the transition to the phase coherent regime is crossed. However applying the KZ formalism to our setup is not straightforward. In a weakly interacting, homogeneous 3D Bose gas, BEC occurs when the 3D phase-space-density reaches the critical value g 3/2 (1). For our quasi-2D geometry, transverse condensation occurs when the 3D phase-spacedensity in the central plane z = 0 reaches this value. At the transition point, the KZ formalism relates the size of phase-coherent domains to the cooling speedṪ . For fast cooling, KZ theory predicts domain sizes for a 3D fluid that are smaller than or comparable to the thickness a z of the lowest vibrational state along z; it can thus provide a good description of our system. For a slower cooling, coherent domains much larger than a z would be expected in 3D at the transition point. The 2D nature of our gas leads in this case to a reduction of the in-plane correlation length. In the slow cooling regime, we thus expect to find an excess of topological defects with respect to the KZ prediction for standard 3D BEC.
More explicitly we expect for fast cooling, hence short t evap , a power-law decay N v ∝ t −d evap with an exponent d given by the KZ formalism for 3D BEC. The fit of this function to the measured variation of N v for t evap ≤ 250 ms leads to d = 0.69 (17) (see Fig. 6 a). This is in good agreement with the prediction d = 2/3 obtained from the critical exponents of the so-called "F model" [21], which is believed to describe the universality class of the 3D BEC phenomenon. For comparison, the prediction for a pure mean-field transition, d = 1/2, is notably lower than our result.
For longer t evap , the above described excess of vortices due to the quasi-2D geometry should translate in a weakening of the decrease of N v with t evap . The nonzero plateau observed in Fig. 6 a for t evap ≥ 250 ms may be the signature of such a weakening. Other mechanisms could also play a role in the nucleation of vortices for slow cooling. For example due to the box potential residual rugosity, the gas could condense into several independent patches of fixed geometry, which would merge later during the evaporation ramp and stochastically form vortices with a constant probability.
Lifetime of the topological defects. The variation of the number of vortices N v with the hold time t hold allows one to study the fate of vortices that have been nucleated during the evaporation. We show in Fig. 6 b the results obtained when fixing the evaporation to a short value t evap = 50 ms. We observe a decay of N v with the hold time, from N v = 2.3 initially to 0.3 at long t hold (2 s). To interpret this decay, we modeled the dynamics of the vortices in the gas with two ingredients: (i) the conservative motion of a vortex in the velocity field created by the other vortices, including the vortex images from the boundaries of the box potential [48], (ii) the dissipation induced by the scattering of thermal excitations by the vortices, which we describe phenomenologically by a friction force that is proportional to the nonsuperfluid fraction of atoms in the gas [47]. During this motion, a vortex annihilates when it reaches the edge of the trap or encounters another vortex of opposite charge. The numerical solution of this model leads to a non-exponential decay of the average number of vortices, with details that depend on the initial number of vortices and their locations.
Assuming a uniform random distribution of vortices at the end of the evaporation, we have compared the predictions of this model to our data. It gives the following values of the two adjustable parameters of the model, the initial number of vortices N v,0 = 2.5 (2) and the superfluid fraction 0.94 (2); the corresponding prediction is plotted as a continuous line in Fig. 6 b. We note that at short t hold , the images of the clouds are quite fuzzy, probably because of non-thermal phononic excitations produced (in addition to vortices) by the evaporation ramp. The difficulty to precisely count vortices in this case leads to fluctuations of N v at short t hold as visible in Fig. 6 b. The choice t hold = 500 ms in Fig. 6 a was made accordingly.
The finite lifetime of the vortices in our sample points to a general issue that  (2) in the presence of a phenomenological damping coefficient [47]. The inferred superfluid fraction is 0.94 (2). Confidence ranges on these parameters are obtained from a χ 2 -analysis.
one faces in the experimental studies on the KZ mechanism. In principle the KZ formalism gives a prediction on the state of the system just after crossing the critical point. Experimentally we observe the system at a later stage, at a moment when the various domains have merged, and we detect the topological defects formed from this merging. In spite of their robustness, the number of vortices is not strictly conserved after the crossing of the transition and its decrease depends on their initial positions. A precise comparison between our results and KZ theory should take this evolution into account, for example using stochastic mean-field methods [49][50][51][52].

Discussion
Using a box-like potential created by light, we developed a setup that allowed us to investigate the quantum properties of atomic gases in a uniform quasi-2D configuration. Thanks to the precise control of atom number and temperature, we characterized the regime for which phase coherence emerges in the fluid. The uniform character of the gas allowed us to disentangle the effects of ideal gas statistics for in-plane motion, the notion of transverse condensation along the strongly confined direction, and the role of interactions. This is to be contrasted with previous studies that were performed in the presence of a harmonic confinement in the plane, where these different phenomena could be simultaneously present in the non-homogeneous atomic cloud. For the case of a weakly interacting gas considered here, our observations highlight the importance of Bose statistics in the emergence of extended phase coherence. This coherence is already significant for phase-space densities D 0 ∼ 3 − 4, well below the values required for (i) the superfluid BKT transition and (ii) the full Bose-Einstein condensation in the ground state of the box. For our parameters, the latter transitions are expected around the same phase-space-density (∼ 8 − 10) meaning that when the superfluid criterion is met, the coherence length set by Bose statistics is comparable to the box size.
By cooling the gas further, we entered the regime where interactions dominate over thermal fluctuations. This allowed us to visualize with a very good contrast the topological defects (vortices) that are created during the formation of the macroscopic matter-wave, as a result of a Kibble-Zurek type mechanism. Here we focused on the relation between the vortex number and the cooling rate. Further investigations could include correlation studies on vortex positions, which can shed light on their nucleation process and their subsequent evolution [53].
Our work motivates future research in the direction of strongly interacting 2D gases [19], for which the order of the various transitions could be interchanged. In particular the critical D for the BKT transition should decrease, and reach ultimately the universal value of the "superfluid jump", D = 4 [54]. In this case, the emergence of extended coherence in the 2D gas would be essentially driven by the interactions. Indeed once the superfluid transition is crossed, the one-body correlation function is expected to decay very slowly, G 1 (r) ∝ r −α , with α < 1/4. It would be interesting to revisit the statistics of formation of quench-induced topological defects in this case, for which significant deviations to the KZ power-law scaling have been predicted [55,56].

Methods
Characterization of the box-like potential. We create the box-like potential in the xy plane using a laser beam that is blue-detuned with respect to the 87 Rb resonance. At the position of the atomic sample, we image a dark mask placed on the path of the laser beam. This mask is realized by a metallic deposit on a wedged, anti-reflective coated glass plate. We characterize the box-like character of the resulting trap in two ways. (i) The flatness of the domain where the atoms are confined is characterized by the root mean square intensity fluctuations of the inner dark region of the beam profile. The resulting variations of the dipolar potential are δU/U box ∼ 3%, where U box is the potential height on the edges of the box. The ratio δU/k B T varies from ∼ 40% at the loading temperature to ∼ 10% at the end of the evaporative cooling (ramp of U box , see below). In particular, it is of ∼ 20% at the transverse condensation point for the configuration in which the vortex data have been acquired. (ii) The sharp spatial variation of the potential at the edges of the box-like trapping region is characterized by the exponent α of a power-law fit U (r) ∝ r α along a radial cut. We restrict the fitting domain to the central region where U (r) < U box /4 and find α ∼ 10-15, depending on the size and the shape of the box.
Imaging of the atomic density distribution. We measure the atomic density distribution in the xy plane using resonant absorption imaging along z. We use two complementary values for the probe beam intensity I. First we use a conventional low intensity technique with I/I sat ≈ 0.7, where I sat is the saturation intensity of the Rb resonance line, with a probe pulse duration of 20 µs. This procedure enables a reliable detection of low density atomic clouds, but it is unfaithful for high density ones, especially in the 2D geometry due to multiple scattering effects between neighboring atoms [57]. We thus complement it by a high intensity technique inspired from [58], in which we apply a short pulse of 4 µs of an intense probe beam with I/I sat ≈ 40. ToF bimodality measurements (where the cloud is essentially 3D at the moment of detection) were performed with the low intensity procedure. This was also the case for the matter-wave interference measurements, for which we reached a better fringe visibility in this case. In situ images in Fig. 1 and all data related to vortices (e.g. Fig. 5 a,b) in the strongly degenerate gases were taken with high intensity imaging. We estimate the uncertainty on the atom number to be of 20%.
Ideal gas description of trapped atomic samples. We consider a gas of N non-interacting bosonic particles confined in a square box of size L in the xy plane, and in a harmonic potential well of frequency ν z along z. The eigenstates of the single-particle Hamiltonian are labelled by three integers j x , j y ≥ 1, j z ≥ 0: where a z = (h/mν z ) 1/2 /(2π) and χ j is the normalized j-th Hermite function. Their energies and occupation factors are where µ < 0 is the chemical potential of the gas and N = j n j . The average value of any one-body observableÂ can then be calculated: Estimation of the interaction energy for weakly interacting gases. We estimate the local value of the interaction energy per particle int = (2π 2 a/m)ρ (3D) (r), where a = 5.1 nm is the 3D scattering length characterizing s-wave interactions for 87 Rb atoms and ρ (3D) (r) the spatial 3D density estimated using the ideal gas description. It is maximal at trap center r = 0. For example, using a typical experimental condition with N = 40 000 atoms in a square box of size L = 24 µm at T = 200 nK, we find a maximal 3D density of ρ (3D) (0) = 13.8 µm −3 . The mean-field interaction energy for an atom localized at the center of cloud is then int = k B ×2.1 nK. We note that int is negligible compared to k B T and hν z for all atomic configurations corresponding to the onset of an extended phase coherence. In this case the interactions play a negligible role in the 2D ToF expansion that we use to reveal matter-wave interferences.
Temperature calibration. All temperatures indicated in the paper are deduced from the value of the box potential, assuming that the evaporation barrier provided by U box sets the thermal equilibrium state of the gas. This hypothesis was tested, and the relation between T and U box calibrated, using atomic assemblies with a negligible interaction energy. For these assemblies, we compared the variance of their velocity distribution ∆v 2 obtained from a ToF measurement to the prediction of Eq. (8).
The calibration obtained from this set of measurements can be empirically written as where the values of the dimensionless parameter η and of the reference temperature T 0 slightly depend on the precise shape of the trap. For the square trap of side 24 µm we obtain T 0 = 191 (6) nK and η = 3.5 (3). The reason for which T saturates when the box potential increases to infinity is due to the residual evaporation along the vertical direction, above the barrier created by the horizontal Hermite-Gauss beam.
Power exponent for fitting N c . We estimate the behavior of ∆(N ) at fixed T using Bose law for a ideal gas. We compute from (8) the equilibrium velocity distributionρ(v). Then we estimate the spatial density after a ToF of duration τ (for a disk trap of radius R) via ρ(r) ∝ρ(r/τ ) * Θ(r R) where * stands for the convolution operator and Θ for the Heaviside function. We fit ρ(r) to a double Gaussian and compute the atom fraction in the sharpest Gaussian ∆, similarly to the processing of experimental data. To simulate our experimental results, we consider ν z = 350 Hz, R = 12 µm, τ = 14 ms and T varying from 100 to 250 nK. For a given T , we record ∆ while varying the total atom number N from 0.06 to 4 times the theoretical critical number for BEC ⊥ N c,th = ζ(π 2 /6)A/λ 2 T (see Sup. Mat.). We fit ∆(N ) between N min = 0.06 N c,th and a varying N max in 1.1 -4 N c,th , to f (N ) = (1 − (N c /N ) α ) with N c as a free parameter and a fixed α. For all considered T and N max , choosing α = 0.6 provides both a good estimate of N c (between 0.93 and 0.99 N c,th ) and a satisfactory fit (average coefficient of determination 0.94).
Chemical potential in the degenerate interacting regime. To compute the chemical potential µ of highly degenerate interacting gases, we perform a T = 0 mean-field analysis. We solve numerically the 3D Gross-Pitaevskii equation in imaginary time using a split-step method, and we obtain the macroscopic ground state wave-function ψ(r). Then we calculate the different energy contributions at T = 0 -namely the potential energy E pot , the kinetic energy E kin and interaction energy E int -by integrating over space : with ω z = 2πν z and = h/(2π). We obtain the value of the chemical potential µ by taking the derivative of the total energy with respect to N and subtracting the single-particle ground state energy: In the numerical calculation, we typically use time steps of 10 −4 ms and compute the evolution for 10 ms. The 3D grid contains 152 × 152 × 32 voxels, with a voxel size 0.52 × 0.52 × 0.26 µm 3 .
Analysis of the density holes created by the vortices. We first calculate the normalized density profile ρ/ρ where the averageρ is taken over the set of images with the same ToF duration τ . Then we look for density minima with a significant contrast and size. Finally for each significant density hole, we select a square region centered on it with a size that is ∼ 3 times larger than the average hole size for this τ . In this region, we fit the function to the normalized density profile, where A 0 accounts for density fluctuations. We also correct for imaging imperfections (finite imaging resolution and finite depth of field) by performing a convolution of the function defined in Eq.(14) by a Gaussian of width 1 µm, which we determined from a preliminary analysis.
Supplementary material:

Transverse condensation and 2D coherence
Most of the experimental data have been taken with a confinement frequency along the z axis ν z = 365 Hz = (k B /h) 18 nK. These data thus lie in the regime ζ = k B T /hν z > 1. However, thanks to Bose statistics, one can still reach a situation where the surface density ρ 2D 0 associated to the ground state of the z motion |j z = 0 is comparable to the total surface density ρ 2D tot . Indeed the surface density ρ 2D exc = ρ 2D tot −ρ 2D 0 that can be accumulated in the excited states of the z motion is bounded [1]. When ρ 2D tot largely exceeds this bound, essentially every additional atom accumulates in the ground state of the z motion |j z = 0 . In this supplementary material, we analyse this transverse condensation phenomenon (BEC ⊥ ) and show that the coherence length in the xy plane is also significantly affected when BEC ⊥ occurs.

Thermal equilibrium of an ideal Bose gas in a quasi-2D geometry
We consider a gas of N non-interacting bosonic particles confined in a square box of size L in the xy plane and in a harmonic potential well of frequency ν z along z.
We use Dirichlet boundary conditions in the xy plane so that the eigenstates of the single-particle Hamiltonian are labeled by two quantum numbers j x , j y describing the state in the xy plane: 1 L sin(πj x x/L) sin(πj y y/L), j x , j y strictly positive integers, (15) and the quantum number j z ≥ 0 describing the vibrational state along z. Putting by convention the energy of the ground state at zero, the energies and occupation factors of these energy levels are where j = (j x , j y , j z ), = h/(2π), µ < 0 is the chemical potential of the gas and the total atom number in the gas is The number of atoms in the j z vibrational state is so that the 2D phase-space-density associated to that state is Turning the discrete sum over j x , j y into an integral in the limit where L/λ T is very large, we find where Z = µ/k B T and ζ = k B T /hν z . The total 2D phase-space-density is 2 Transverse condensation in a quasi-2D geometry In order to investigate the transverse condensation phenomenon, we use a treatment very similar to that of usual 3D Bose-Einstein condensation. We focus on the case ζ 1 where many vibrational states along the z direction are populated. Using the semi-classical approximation that consists in replacing the discrete sum in (22) by an integral over j z , we calculate the total phase-space-density : Since Z < 1 and g 2 (Z) remains finite when Z → 1 [g 2 (1) = π 2 /6], this semi-classical approximation leads to the paradoxical result that for a given ζ, the total 2D phasespace-density is bounded from above by ζπ 2 /6. The paradox is lifted by noticing that when the fugacity Z approaches 1, the population of the lowest vibrational state |j z = 0 is not properly accounted for when one replace the discrete sum in (23) by the integral (24). More precisely within this semi-classical approximation, when the total phase-space-density approaches the value ζπ 2 /6, the result above must be replaced by with for Z ≈ 1 : Therefore, when the total phase-space-density D 2D tot is significantly larger than ζ π 2 /6, the phase-space-density D 2D exc saturates and additional atoms accumulate essentially in the |j z = 0 state (see Fig. 7) [1]. tot . The calculation is made for an ideal gas confined in a square box in the limit L/λ T → ∞. The values of ζ are: 0.1 (olive), 1 (magenta), 5 (cyan), 10 (black), 15 (red), 20 (blue). The dotted lines indicate the critical phase-space-density ζ π 2 /6 (same color code).
Let us estimate the 2D phase-space-density associated to the ground state of the z motion when D 2D tot reaches the threshold for transverse condensation D 2D tot,c = ζ π 2 /6.
At this point, the population of the state |j z = 0 is significantly different from that of |j z = 1 , so that it cannot be accounted for properly by the integral (24). This implies that the chemical potential is on the order of −hν z , i.e., Z ∼ e −1/ζ and we thus predict More precisely, a numerical calculation (see Table 1) gives at the condensation point for ζ between 5 and 20.
In the limit ζ → ∞, this transverse condensation phenomenon (BEC ⊥ ) constitutes a phase transition (see Fig. 8): denoting the transversely condensed fraction as f 0 = D 2D 0 /D 2D tot and the reduced total phase-space-density as x = D 2D tot /D 2D tot, c , we find in this limit:

Central 3D density and transverse condensation
At thermal equilibrium the 3D density at a point z is given by where χ j is the j-th normalized Hermite function and a z = (1/2π) h/mν z . In the limit λ T L, we replace again the sum over j x , j y by an integral and get  In the following we are interested in the value of the 3D density in the central plane z = 0 where it is maximum. We recall that The value of the 3D phase-space-density in the plane z = 0 is thus: which is plotted in Fig.9, as a function of D 2D tot . As above we calculate the contribution of the excited states of the z motion by replacing the sum over j z = 2n by an integral. Using for n 1 Above the transverse condensation threshold, the contribution of the j z = 0 must be calculated separately and the total phase-space-density reads with for Z ≈ 1 : : Total 3D phase-space-density in the plane z = 0, as a function of the total 2D phase-space-density (same color code as in Fig. 7). The horizontal line corresponds to D 3D tot (z = 0) = 2.612, which corresponds to the threshold for BEC in a uniform 3D gas.
In other words, transverse condensation occurs when the central 3D phase-space density approaches the value 2.612, which would corresponds to a "true" Bose-Einstein condensation for a 3D gas (Fig. 9).

Coherence length in the xy plane
So far we have only addressed the thermodynamics of the atoms along the z axis. However when BEC ⊥ occurs for ζ 1, the coherence length in the xy plane is also affected and it can become significantly larger than λ T . In order to characterize this in-plane coherence, we consider the one-body correlation function for r in the xy plane. At thermal equilibrium, this is equal to Note that in this section we turn to periodic boundary conditions in order to simplify the calculations. Thus we label the single particle states by their in-plane momentum k = (2π/L)(j x , j y ), with j x , j y ∈ Z. The integral over k can be calculated by expanding n k,jz as a power series in the fugacity Z which is an infinite sum of Gaussian functions with a width increasing with n. For large r, it is useful to determine the dominant term in this sum over n, which for a given j z , is obtained approximately for and takes a value approximately proportional to Let us focus on the ground state of the z motion j z = 0 and consider the regime where Z ≈ 1 so that D 2D 0 = − ln(1 − Z) 1 and ln(1/Z) ≈ e −D 2D 0 . We then find the approximate value contribution of j z = 0 : which coincides with the result derived by another method in [2]. A rough estimate for at the transverse condensation point can be obtained by plugging the approximate value D 2D 0 ∼ ln ζ obtained at Eq. (29) into this expression: where a ho = (2π) −1 h/mν z . A more precise estimate is given in Table 1 for the relevant range of values for ζ. To obtain these results, we computed numerically the variations of G 1 using (42), and we looked for the wings of this function. More precisely we considered the point r 20 where the function G 1 is divided by 20 with respect to its value in r = 0. At this point we define as If G 1 had an exponential variation for all r, this quantity would take the same value independently of the location r where it is calculated. For a non-strictly exponential G 1 , the present definition is a good compromise between considering the far wings of G 1 in order to monitor the appearance of an extended coherence, and restricting to sufficiently small values of r so that the values of G 1 are still significant. The variation of /λ T with D 2D tot for various values of ζ is shown in Fig. 10. For Boltzmann statistics G 1 (r) = exp(−πr 2 /λ 2 T ), so that the expected value of at r = r 20 is Bolt. = 0.16 λ T . We see in table 2 that at BEC ⊥ , the value of is  (47) with the total 2D phase-space-density (same color code as in Fig. 7).
increased by a factor 6 − 9 for ζ = 10 − 20 with respect to Bolt. . More precisely the result (46) states that at the threshold for BEC ⊥ , the coherence length in the xy plane is comparable to the size of the ground state along the z direction. The physical meaning of (46) is related to the fact that the 3D phase-space-density in the plane z = 0 reaches the value 2.612 when BEC ⊥ occurs. For a uniform, infinite 3D gas with this spatial density, the coherence length would diverge, signaling the occurrence of a true Bose-Einstein condensation. Here the z confinement limits the extension of the coherent part of the gas along z to a ho , which prevents the divergence of and limits its value also to a ho in the xy plane. In the regime ζ 1, the appearance of a large coherence length in the xy plane and the occurrence of transverse condensation for the z degree of freedom are thus linked.

Finite size effects
So far we assumed that the size of the box in the xy plane was arbitrarily large compared to the other length scales of the problem, λ T in particular, so that the sum over j x , j y could safely be replaced by an integral. For a finite-size box this approximation ceases to be valid when the phase-space-density becomes large enough. A full 3D Bose-Einstein condensation (BEC full ) can then take place, with a macroscopic accumulation of particles in the single-particle ground state j x = j y = 1, j z = 0. The existence of two successive condensations when the phase-space-density increases, first BEC ⊥ and then BEC full , was highlighted in [1].
BEC full is expected to occur when the chemical potential is chosen smaller (in absolute value) than the gap 3h 2 /(8mL 2 ) between the true ground state of the box and the first excited states. At the point where BEC full occurs, the 2D phase-space-density associated to the |j z = 0 level is This value is notably larger than the value D 2D 0 = ln(ζ) when BEC ⊥ occurs if meaning simply that the gas must have a flat shape. An example is given in Fig. 11 where we plot the result of a calculation summing the populations of the individual quantum states given in (16). The calculation is performed for a typical box size L = 200 λ T and for ζ = 5. The two successive transitions are clearly visible on this figure. First when D 2D tot ≈ D 2D tot, c with D 2D tot, c = ζπ 2 /6 ≈ 8.2, the phase-space-density associated to the excited state of the z motion saturates at a value close 1 to D 2D tot, c , which is the signature of BEC ⊥ . Then when D 2D tot reaches a value around D 2D tot = ζπ 2 /6 + ln (4/3π) L 2 /λ 2 T ≈ 18, the fraction of atoms f BEC occupying the overall ground state starts to be significant. It is clear on this figure that there exists a domain of values of D 2D tot for which the population of the excited states of the z motion is saturated, with no macroscopic occupation of the single particle ground state in the box.
For this typical box size L = 200λ T , we infer from Eq. (48) that BEC full occurs when D 2D 0 reaches the value ∼ 9.7. For Rb atoms and a trapping frequency ν z in the range 300 − 1500 Hz, this phase-space-density is similar to the one at which the (interaction-induced) BKT superfluid transition occurs (D 2D BKT = 8 − 10) [2,3]. Therefore the relevance of the previous sections of this supplementary material, based on calculations for a non-interacting gas in the limit L λ T (no BEC full ), is limited to the region of parameters where D 2D 0 8. This is precisely the region where the emergence of coherence studied in the paper occurs ( D 2D 0 ∼ 4).