Thermalization and Nonequilibrium Steady States in a Few-Atom System

1 Introduction

The description and investigation of nonequilibrium systems offer a vast field of open questions for both theory and experiment [1], [2]. Specifically, the thermalization of few-particle systems out of equilibrium and the properties of the steady state that is approached are studied intensely [3], [4]. Experimentally, when employing large atom numbers, atomic interaction typically drives the thermalization of the sample. A paradigmatic application of this thermalization is the widely used technique of evaporative cooling [5]. For few-particle samples, by contrast, interaction between the particles is usually negligible because the extremely low particle number results in a low density of the sample and, hence, negligible scattering cross sections. Nevertheless, they provide an exciting area of research because, for such few-particle systems, the thermodynamic limit eventually breaks down, thus, opening a way to investigate fundamental questions like the definition of entropy [6]. A great advantage of ultracold atom systems is that, even in the few-particle limit, interaction between the particles can be tuned to a significant level by, e.g. employing a strong confinement of the sample [7].

Here, we investigate an isolated system of few ultracold cesium (Cs) atoms, forming a nonthermal steady state. Based on experimental data, we devise a numerical model indicating a parameter regime, where thermalization dynamics of a few-particle system can be resolved. Key ingredients are, first, the large s-wave scattering length for Cs–Cs collisions and, second, low temperatures created by sympathetic cooling of the Cs sample with an ultracold cloud of Rb atoms.

2 Microscopic Model

The interaction between ultracold atoms is characterized by the scattering cross-section σ for atomic collisions. In the effective-range approximation, this scattering cross section reads

for bosons [8], [9]. There, a is the Cs–Cs s-wave-scattering length, which, for Cs, takes relatively large values in the range of 103a0 with a0≈0.5 Å being Bohr’s radius, k is the collision wave-vector, and r_eff is the effective range characterizing the scattering potential. The effective range is accessible by the approximation [10]

where Γ is the Gamma function and

is the asymptotic mean scattering length, which is a property of the scattering potential [11]. Here, μ is the reduced mass of the colliding atoms, and C6=6890.48a.u. is the van der Waals constant of the scattering potential [10]. These values yield an asymptotic mean scattering length of a¯c=97 a0 and an effective range of reff=300 a0. The dependence of the scattering cross section on the collision energy Ec=ℏ2k2/2μ, which results from (1) is shown in Figure 1. The dependence features two limiting cases: For small energies and, hence, small wave vectors k, the cross section approaches σ(Ec=0)=8πa2. For large energies Ec≫0, the cross section asymptotically approaches σ=8π/(k2+reff2k4/4) from below, even for very large scattering lengths a. As a result, in order to harness a large s-wave-scattering length, sufficiently low collision energies and hence temperatures need to be realized experimentally. This provides a pathway towards large scattering cross sections and, thereby, significant interaction even in the few-particle regime.

Figure 1:

The Cs–Cs scattering cross-section σ (blue line) depends strongly on the collision energy E_c in the relevant energy range. The dashed green line indicates the scattering cross-section σ given by the scattering length a(B=0)=−3068 a0 for Ec→0. For E_c > 0, an upper bound arises from the wave-vector (solid green line). The vertical red lines indicate the relevant energy scales for the simulation with mean collision energy at temperature T = 0.2 μk and T = 1.6 μk.

For a thermal sample at temperature T, the mean collision energy may be estimated by ⟨Ec⟩=μv¯2/2=4kBT/π, with the mean relative velocity v¯=(16kBT/πmCs)1/2 of the gas [12]. In the following, we assume Cs atoms to be prepared initially at a temperature of T=200 nK and subjected to a quench of the dipole-trap potential, which increases the total energy of the sample corresponding to an effective temperature of T = 1.6 μK. The mean collision energies for these situations are indicated in Figure 1 by vertical lines.

Figure 2:

Dependence of the Cs–Cs s-wave-scattering length (blue line) on the magnetic background field B, reproduced using the values given in [13]. The Feshbach resonances are indicated by red vertical lines at their corresponding resonance field B0,i*. The dashed blue line indicates the background scattering length a_bg. The shaded gray area corresponds to a hypothetical negative magnetic field strength, which is not accessible.

The above considerations apply for Cs atoms at small magnetic fields, which are in the order of the earth’s magnetic field of a few hundred Milligauss. However, by employing magnetic Feshbach resonances, the s-wave-scattering length a can be tuned within a broad range around the background scattering length abg=1875a0 [13]. At every Feshbach resonance i, which is characterized by its resonance field B0,i*, the scattering length features a pole and, therefore, allows tuning of the scattering length in the vicinity of the resonance. The analytical formula describing the dependence of the scattering length on the magnetic field reads [13]

where Bi* is the corresponding zero crossing fields that quantify the width of every resonance. The Cs–Cs Feshbach resonances in the magnetic field range up to 100 G were investigated experimentally in [13]. We use the values for Bi* and B0,i*, which were found there for the three dominant resonances to reconstruct the Cs–Cs scattering length a(B), which is shown in Figure 2. This reproduction of data also yields the scattering length for small magnetic fields a(0)=−3068 a0. This is a large absolute value compared to other scattering lengths, which are relevant for, e.g. two-species Rb–Cs mixtures with aRbCs=645 a0 and aRbRb=101 a0 [14], [15] and illustrates the comparably strong Cs–Cs interaction in the experiment. The large negative value stems from a broad Feshbach resonance at a resonance field of −11.1 G [13]. While this negative magnetic field strength cannot be realized physically, the influence of this resonance dominates the characteristics of the scattering length in the field range from zero to 100 G. This includes a zero-crossing around B=18 G. Overall, the broad Feshbach resonance offers the opportunity to experimentally tune the Cs–Cs interaction from large negative to positive values and also facilitates the suppression of interaction by setting the magnetic field to the zero crossing field. Thereby, the interaction within few-particle samples of Cs atoms can be tightly controlled.

In order to simulate the long-time evolution based on microscopic dynamics, we solve the microscopic equations of motion for every Cs atom numerically by a third-order Runge–Kutta method [16]. The effect of s-wave scattering is included in the simulation by a simple hard-sphere model: The relative distance rij=|r→i−r→j| between all pairs of particles i,j=0…N is evaluated after every simulation time step. A scattering event is triggered, when the inter-atomic distance is smaller than the scattering cross-section rij<(σij/π)1/2. There, the cross-section σij(kij) is calculated from (1) for every atomic pair, taking into account the scattering wave-vector kij=μvij/ℏ, which is determined by the relative velocity vij=|v→i−v→j| of the particles. When a scattering event is triggered, the relative velocity in the center-of-mass frame of the particles is redistributed to a randomly drawn angle. This simulates the momentum exchange and takes into account the isotropic nature of the s-wave scattering. This microscopic treatment of the scattering processes implies a quadratically growing complexity of the problem 𝒪(N2) with the particle number of the sample N. While this renders the method in its present implementation unfeasible for typical many-particle systems, the computational effort for few-particle systems is small: For a typical atom number of N = 10, the computation time used for the implementation of the scattering events is smaller than the time that is needed for finding the solutions to the equations of motion. Typical simulations of several thousand trajectories require computation times around 1 day on a standard office desktop computer. Therefore, the simple numerical implementation of s-wave scattering presented here is a well-performing approach to model the microscopic dynamics of interacting few-particle systems.

Figure 3:

(a) Illustration of the dipole trapping beams, which create the confinement of the atoms (yellow) and the optical lattice that can be employed for axial localization of the atoms (blue). (b) Potential cut along the z-axis. The solid line indicates the combined potential from horizontal and vertical beam. The dashed line is the corresponding harmonic approximation, and the dotted line is the harmonic approximation of only the horizontal dipole-trap laser beam. The green shaded distribution indicates the initial Cs distribution in the simulation.

Figure 4:

Numerically simulated dynamics of the axial position distribution f(z) for an initially displaced thermal Cs sample having (a) N = 1 and (c) N = 10 atoms. The atomic density is represented by the colormap and plotted versus the time and the position z in the trap. The regions indicated by the vertical red, orange, and blue lines starting at 0.0 s, 0.5 s, and 2.0 s, respectively, are used to compute the averaged axial density distributions, which are shown in (b) and (d). The insets illustrate the axial phase space distribution at the start of every region. The horizontal axis of the phase space illustration displays the z-axis in the range of ±200 μm, while the vertical axis corresponds to the velocity v_z in the range of ±50 mm/s. (b, d) The mean value is shown as solid lines, while the area between the 10^th and 90^th percentile is indicated by colored shading illustrating the time-averaged fluctuations around the mean. Red, orange, and blue indicates the corresponding starting times of 0.0 s, 0.5 s, and 2.0 s, respectively. In both cases, the solid black line indicates the distribution of a thermal state in harmonic approximation with a temperature of T = 1.6 μk, corresponding to the total energy in the system.

3 Steady-States Versus Thermalization Dynamics

For the few-particle dynamics considered here, the optical dipole trapping setup is a crucial component, as it provides the conservative potential, which confines the particles in vacuum [5] and allows to tune the effective scattering rate independently from the scattering length. A schematic view of the most important constituents of the setup considered here is shown in Figure 3a. The horizontal and vertical dipole-trap laser beams are generated by the same laser source at a wavelength of 1064 nm and induce an attractive potential for Cs atoms. The main confinement is created by the horizontal beam with a total trap depth of kB×365 μK at a power of 1 W and a Gaussian beam waist of 21 μm. The crossed horizontal beam with a beam waist of 165 μm and power of 1 W adds a contribution of kB×15 μK, as indicated in Figure 3b. The purpose of this beam is to create a small anharmonicity, which facilitates the dephasing of atomic trajectories, which will be discussed later. This induced anharmonicity is illustrated by the dashed and dotted lines in Figure 3b, which represent the harmonic approximation at the trap center, as well as the harmonic approximation for the potential created exclusively by the horizontal beam, respectively. While the former is relevant for the dynamics at the trap center, the latter characterizes the trapping potential outside of the center. An additional optical lattice potential is created by two counter-propagating laser beams along the horizontal (z-) direction at a wavelength of 790 nm, which creates a repulsive potential for Cs atoms. It, hence, tightly confines Cs atoms axially in the horizontal dipole-trap beam. By applying a relative detuning between the two counter-propagating beams, the interference pattern can be set into motion, and thereby, precise transportation of Cs atoms along the z-axis can be achieved, realizing an atomic conveyor belt [17], [18].

We first numerically investigate a thermal distribution of Cs atoms at a temperature of 200 nK. This temperature can be achieved, for example, by employing sympathetic cooling within a cold cloud of Rb atoms, which can be prepared experimentally at this temperature. We already demonstrated the successful immersion of Cs atoms into cold clouds of Rb atoms in [19], [20], and [21]. By removing the cold Rb bath from the dipole trap by means of a resonant pulse of laser light, a pure Cs sample at such low temperatures can be prepared. Subsequently, all Cs atoms are transported out of the trap center by means of the conveyor-belt optical lattice, preparing them at a distance of 100 μm from the potential minimum. Dynamics is initiated by releasing these Cs atoms, which creates a nonequilibrium state featuring an oscillatory behavior. The Cs mean atom number can be controlled experimentally by varying the Cs MOT loading time, allowing for typical mean atom numbers between 1 and 10 atoms.

Using the numerical simulations described previously, we compute in total 5×103 trajectories. We extract the dynamics of the axial position distribution f(z), which is shown in Figure 4 for two different settings: a non-interacting case is realized by setting the Cs atom number per sample to N = 1, which nullifies the occurrence of scattering events (see Fig. 4a,b). For the interacting case, a sample size of N = 10 atoms is considered, allowing for atomic collisions (see Fig. 4c,d). This range of atom numbers corresponds to values that are routinely accessible in the experiment.

During the first oscillation periods, the evolution of the non-interacting and interacting samples is almost identical, exhibiting the expected oscillation at an amplitude of 100 μm. In the non-interacting case, these oscillations dephase more and more because different parts of the atomic cloud oscillate at different frequencies due to the anharmonicity of the trapping potential. This effect can be observed nicely in the full atomic phase space ρ(z,vz) along the z-axis. The phase-space representation reveals the creation of a spiral-like feature, which is winding up more and more as the evolution time increases (see insets of Fig. 4a). After 2 s of evolution, the system approximately reaches a steady state, where the structure of the wound-up spiral is blurred, and the steady-state phase-space distribution features a donutlike shape. The approach of this steady state can be quantified by analyzing the evolution of the position distribution during a full oscillation period as illustrated by the colored vertical lines in Figure 4. Extracting the mean as well as the 10^th and the 90^th percentiles from the simulation data characterizes the distribution and its fluctuations. We clearly observe a reduction in fluctuations over time and, thereby, the approach to a double-peak position distribution. The mean of the distribution, however, does not change significantly during the evolution. The approached steady state is nonthermal, illustrated by the difference to the black, thermal distribution at the corresponding total energy of the system in Figure 4b. For the interacting case, however, the steady state, which is approached, changes drastically: While also, here, the effect of dephasing and the corresponding reduction of fluctuations are visible, additionally, the scattering events between the atoms lead to a redistribution of the phase-space density and particularly populate the center of phase space (see Fig. 4c). This redistribution creates a steady-state phase-space distribution with Gaussian shape. Quantitatively, the scattering events manifest themselves in a dynamic of the mean position distribution, which is absent in the non-interacting case. As seen in Figure 4d, the mean position distribution approaches the thermal distribution

for a temperature T and a harmonic trap with ωz=2π×39 Hz, which correspond to the values of the harmonic approximation at the trap centre. The temperature T=⟨Etot⟩/3kB=1.6 μK is derived from the mean total energy ⟨Etot⟩ of the nonequilibrium state. The numerical considerations, hence, demonstrate how an initially nonthermal sample of few Cs atoms relaxes to a thermal state in the presence of interaction. The proposed interaction time of 2 s is experimentally achievable because typical lifetimes for Cs atoms trapped in the optical dipole trap exceed 10 s in our experiment.

3.1 Accessing the Atomic Total Energy Distribution

Besides the atomic density distribution along the z-axis, which is measured directly in the experiment, also the total energy distribution of the atoms is potentially accessible. Employing the method of adiabatic lowering, the cumulated total energy distribution can be measured by means of counting the remaining atoms in the trap after the potential was adiabatically lowered to various values [17], [22]. From the numerical data, the total energy distribution can be analyzed quantitatively, similar to the approach employed for the position distribution presented previously. The corresponding distributions for the non-interacting and interacting cases are shown in Figure 5. When atomic interaction is absent, the total energy of every particle is a constant of motion in the conservative potential. As a result, the fluctuations of each distribution vanish, and furthermore, the distributions at the three evolution times of t = 0 s, 0.5 s, and 2 s are identical. The presence of interaction, however, allows for redistribution of energy among the atoms and causes nonvanishing fluctuations reflecting the effect of the scattering events. Under the influence of interaction, the atomic total energy distribution converges to the distribution

(6)f(Etot)=12Etot2(kBT)3exp(−EtotkBT),

Figure 5:

Quantitative analysis of the total energy distribution for slices around t = 0 s, 0.5 s, and 2 s, displayed in red, yellow, and blue, respectively. Evaluation is analogous to Figure 4b and d. (a) Non-interacting setting. The total energy distribution shows no time dependence, and all three curves are identical. The yellow and blue curves were shifted upwards by 0.01 each in order to make the curves distinguishable. (b) Interacting setting with N = 10 atoms per experiment. The total energy distribution relaxes towards the distribution expected for a thermal state at temperature T = 1.6 μk (black line).

with the temperature T = 1.6 μk, which is expected for a thermal state. In contrast to the position distribution, the total energy distribution is not influenced by the initially excited oscillation along the z-axis, as the total energy is conserved for this oscillation. The fluctuations, which are observed for the interacting case in Figure 5, rather directly reflect the influence of the scattering events on the energy distribution. The convergence of the energy distribution towards a thermal distribution is similar to the convergence of the position distribution presented in Figure 4. Already after t = 0.5 s, a significant convergence is observed, and at t = 2 s, the system has almost completely approached the thermal state. The measurement of the atomic total energy distribution, hence, may be used as a complementary method to observe the thermalization of the sample.

Figure 6:

Experimental analysis. (a) Experimental setup and imaging data. The Cs atoms (blue circles) are illuminated in the dipole trapping potential (orange) by an optical molasses (red). The axial position of the atoms is pinned by an additional optical lattice dipole trap potential (blue). The fluorescence light of the atoms is collected by a microscope objective and magnified onto an EMCCD camera chip. This yields the atomic in situ fluorescence (colormap), which is employed to extract the atomic positions; data shown are a randomly selected data set. Moving the atomic distribution by means of the optical conveyor-belt lattice, subsequent images (1–4) are used to extend the field of view. This facilitates the analysis of broad Cs distributions. (b) The atomic positions from the stitched images are binned in a histogram for many experimental repetitions to yield the in-trap density distribution (red bars). The experimental data shown here corresponds to the starting point of the dynamics presented in Figure 7, which is a peaked distribution, displaced from the trap center.

Figure 7:

Dynamic of a few Cs atom sample. (a) Timeline and measured Cs distributions resolving the in-trap dynamics within one oscillation period for starting times of 0, 50, 100, and 200 ms. (b–e) The mean (solid line) as well as the minimum and maximum (shaded area) over the three distributions around 0, 50, 100, and 200 ms, respectively, analogous to the analysis of the slices presented in Figure 4b and d. The dephasing of the initial distribution is nicely visible, and after 200 ms, the residual dynamics on the timescale of the axial trapping period is in the order of the experimental uncertainty. This indicates the approach of the nonequilibrium steady state. This interpretation is supported by strongly reduced fluctuations for dephased Cs distributions.

4 Comparison to Experimental Results

Atomic position distributions are measured experimentally by means of fluorescence imaging as illustrated in Figure 6. To this end, the atoms are trapped in the horizontal dipole trap with the overlapped optical lattice, which provides axial confinement and, thereby, resolution of the atomic position along the z-axis. The atomic fluorescence is excited by applying an optical molasses, consisting of cooling light at a detuning of Δcool=−12MHz to the Cs F=4⇒F′=5 transition of the D2 line and repumping light, which drives resonantly the F=3⇒F′=4 transition. The total intensities of cooling and repumping light are Icool=4.7 mW/cm2≈1.7Isat and Irep=1.1 mW/cm2≈0.4Isat. The fluorescence of the atoms is collected with a microscope objective and imaged onto an EMCCD camera with a resolution of 2.2 μm. Single Cs atoms appear on the fluorescence image as bright peaks and, therefore, allow to extract the position of every single atom from the images. For many repetitions of the experiment, the acquired atomic positions are binned into a histogram, thereby yielding the atomic density distribution along the z-axis as shown in Figure 6b.

Figure 7 shows the dynamics of such density distributions for Cs atoms with a temperature of approximately 10 μk caused by the magneto-optical trap that is employed to produce the atomic sample. In order to account for this relatively large temperature, a horizontal dipole trap power of 3 W is used instead of the simulated value of 1 W. While at these temperatures, no significant effect of atomic collisions is expected due to the relation of collision energy and scattering cross section shown in Figure 1, the measurements illustrate the dephasing of the position distribution in the non-interacting case as well as the reduction in fluctuations as the steady state is approached. Here, the Cs atoms are prepared at a distance of roughly 300 μm from the trap center, and the oscillations are observed experimentally by recording the atomic position distribution for three different times with 4 ms delay between subsequent distributions. These delays capture slightly less than half a trap period and are illustrated in the insets of Figure 7. These three distributions are taken at starting times of 0, 50, 100, and 200 ms. By evaluating the mean as well as the minimum and maximum for every set of distributions, experimental results are extracted, which qualitatively match the predictions of the simulations presented in Figure 4. The reduction in fluctuations due to the dephasing of the atomic distribution is clearly visible as well as the emergence of the characteristic double-peak structure, which is a hallmark of the nonequilibrium steady state that is approached in the non-interacting case.

5 Conclusion

In conclusion, we investigated the dynamics of a few-particle system in the non-interacting and interacting cases in order to study the thermalization of nonequilibrium systems far below the thermodynamic limit. The possibility to experimentally tune the atom number allows us to observe the transition from the single-particle regime to the onset of a many-body system. We showed the experimental feasibility of the method by measurement results, which confirm the emergence of a nonthermal steady state for a non-interacting atomic ensemble. While for the simulations presented here classical trajectories of the particles were considered, the cold atom system features an intrinsic quantum nature of the particles that could be exploited to investigate the thermalization of quantum systems [1], [2]. Recent experimental progress on investigating these nonequilibrium systems with ultracold gas experiments includes, for example, the observation of universal dynamics during the thermalization process [4] or the possibility to experimentally resolve the impact of single s-wave-scattering events on the dynamic after a quench [23]. By taking advantage of the excellent control in our cold gas experiment, also the realization of driven systems is within reach [24], [25]. The proposed system, hence, offers the opportunity to study open questions ranging from classical thermodynamics to driven quantum systems in a versatile and well-controlled experimental setting.