1. Introduction
A magneto-optical trap (MOT) is known worldwide as the staple technology for producing cold and confined ensembles of neutral atoms. Conceptually simple, this device is essential in the preparation stages of highly involved tasks, e.g., Rydberg atom quantum computing [
1,
2] and Bose-Einstein condensation [
3,
4]. Additionally, the limit of a large atom number
N presents a wide array of interesting physics, such as random lasing [
5,
6], Anderson localization [
7,
8], subradiance [
9,
10,
11], superradiance [
12,
13,
14] as well as spontaneous oscillations that bear similarities to plasma and stellar phenomena [
15,
16]. The last example falls into the category of MOT instabilities, which are, by their dynamic nature, out of thermal equilibrium. In the stable MOT regime, however, the thermodynamical considerations can be complex, requiring handling with particular care.
Stable MOTs governed by single-atom physics normally possess isotropic Gaussian velocity distributions and, thus, are deemed in thermal equilibrium. Exceptions occur in cases with Sisyphus cooling, where the distributions resemble double Gaussians [
17]. For conservative traps (such as magnetic traps or optical dipole traps, used in, e.g., Bose-Einstein condensation), the evaporative cooling processes depend on short-range collisional interactions (in particular, van der Waals interactions); after few collisions, the atoms thermalize and acquire isotropic Gaussian velocity distributions, while anisotropically trapped [
18]. In the case studied in this paper, light-mediated long-range interactions are considered in a stable MOT containing many atoms. Compared to the case with short-range interactions, where detailed analytical results of the velocity distribution can, in principle, be obtained [
18], the collective forces stemming from long-range interactions (see below) pose a different challenge, making it unclear whether analytical modelling is viable. We invoke here numerical simulations and predict a velocity distribution typical of a nonequilbrium steady state (NESS), which we observe to act as a precursor to the unstable (self-oscillating) regime. In general, such a state can appear in a non-periodically driven system with broken reversibility and is characterized by the existence of nonzero net currents and a nonequilibrium probability distribution [
19]. States of this kind are abundant in nature, existing in processes as distant as biological ones [
20,
21].
The many-atom physics in a MOT become important typically for
atoms [
22]. The collective forces that appear include, e.g., the shadow force [
23] and the rescattering force [
22]. The former one is compressive and is due to an imbalance of the beam intensities in the cloud, caused by their attenuation due to the light’s scattering. The latter force is, on the opposite hand, repulsive and appears as the atoms rescatter the scattered photons. These antagonistic forces are critical for understanding different features induced by multiple-scattering, e.g., the stable-cloud size increase with
N [
22,
24] as well as the occurrence of spontaneous oscillations when the beam detuning is brought close to the atomic resonance [
15,
25].
This article’s main objective is to numerically study the three-dimensional (3D) particle velocity distributions of the stable regime, including how they are impacted as the unstable regime is approached. At large N, velocity anisotropy is observed, which is a signature of a NESS in a MOT. The NESS behavior is studied with respect to all of the main MOT parameters—the atom number N, the beam detuning , the magnetic field gradient along the strong axis, and the intensity I of a single beam. Moreover, given the 3D nature of the numerical simulations, we are permitted to investigate the full spatiotemporal characteristics of the NESS.
This article continues as follows. In
Section 2, the theoretical model employed in the 3D simulations is briefly described as well as the details surrounding its numerical implementation. Then,
Section 3 presents the main results, including the identification of the NESS, its behavior versus the MOT parameters and its spatiotemporal characteristics. Finally, in
Section 4, we conclude and discuss the future perspectives.
2. Theoretical Approach
Our employed theoretical model and its numerical implementation have been detailed in Ref. [
26], and here we provide a brief description. We consider the so-called balanced MOT configuration, where the laser beams are independent and of the same, constant intensities before entering the cloud. The model is based on the hyperfine transition
, with each of the three Zeeman transitions
treated as an independent 2-level system and driven by, respectively,
,
,
polarized light (see Figure 1b in Ref. [
26]). The magnetic field (which splits the Zeeman levels) is linearized under the assumption of the atom position being much smaller than the radius of the MOT coils and the separation between them. The main physical effects of the model are (i) the (Doppler) trapping force; (ii) the diffusion stemming from its fluctuations; (iii) the beam intensity attenuation due to the light’s scattering by the atoms; and (iv) the rescattering caused by the scattered photon exchange between the atoms. Working with a
system (as opposed to a 2-level system) permits a proper description of the features related to the magnetic field and light polarization.
An important feature of our model is it being spatially nonlocal. Particularly, this is due to the atomic cross sections (in the expressions of the main effects) depending on the beam intensities, which are subject to attenuation. Note that the attenuation affects all of the main effects, including itself. In the case of the trapping force, the attenuation’s inclusion results in an additional compression, i.e., the shadow force (see Figure 2a in Ref. [
26]), which is antagonistic to the rescattering (see Figure 2b in Ref. [
26]). The attenuation naturally introduces the beam cross-saturation effect, which is computed using a devised numerical scheme.
Moreover, the modelling is done in 3D, without any assumption of spherical symmetry. As shown and discussed more in
Section 3, this allows us to observe stable clouds in different non-centrosymmetric configurations.
The complete system dynamics in our model are described by the following collisionless Vlasov-type kinetic equation for the atomic phase-space density
(refer to Equation (26) in Ref. [
26]):
where
M is the atomic mass,
is the trapping force,
is the rescattering force, and
D is the momentum diffusion coefficient; the effects depend implicitly on time as well as on the MOT parameters mentioned earlier—
N,
,
,
I. The dependence on
N arises from their dependence on
I that contains information on the cloud density (for
, each particle also contributes to greater repulsion), whereas the dependencies on
and
arise as these are crucial for a proper description of cooling and confining processes. Finding direct numerical solutions to this equation can be deemed impractical considering both local and nonlocal spatial dependencies. Thus, a more sensible yet an equivalent way of solving the system’s dynamics is used by implementing a numerical produce based on a superparticle [
27,
28] treatment. As an important remark, note that the Vlasov equation is known to contain NESS solutions [
29,
30], meaning a MOT can potentially be found in such a state. In
Section 3, we argue that a MOT indeed contains a NESS.
Our numerical implementation proceeds as follows. We start by generating a Gaussian cloud composed of
superparticles representing a larger amount of real particles (here,
to
). The cloud then evolves under the action of the effects of the model. The analysis is done by randomly selecting cloud images (here, 200, unless specified otherwise) after the initial transient time (see Figure 3 in Ref. [
25]).
Note that we here concentrate on the predictions of the stable MOT regime. The unstable (self-oscillating) one has been studied with our numerical model in previous works [
16,
25] showing qualitative agreements with the experimental observations. Nevertheless, the rationale behind picking particular MOT simulation-parameters is based on our observation, to be elaborated in
Section 3, that the MOT contains a unique stable regime (NESS phase) before entering the unstable regime.
3. Main Results
In this section, the main results of our simulations are presented. The 3D velocity distributions of stable MOT clouds are studied, including how they are impacted as the unstable regime is approached. The study is done with respect to all of the main MOT parameters—N, , , I—in this particular order. The spatiotemporal profiles of the 3D velocity distributions are also discussed. The analysis of our observations leads us to argue for the existence of a NESS in a MOT.
We begin by observing
Figure 1, showing distributions for the three velocity dimensions at different values of
N (the dotted lines are Gaussian fits to the histograms). The shape of the histograms is Gaussian due to the stochastic force (proportional to the diffusion coefficient) being of this kind. When above
atoms, the distributions broaden and, more remarkably, the distribution of the
axis (strong magnetic field axis) narrows compared to the remaining ones. We relate these observations to the growing antagonism between the shadow force and the rescattering force (many-atom effects) as
N is increased. Note that different parts of the large clouds (
) possess varying degrees of velocity anisotropy and even have its sign switched when the atoms near the cloud edge are considered. This is illustrated in
Figure 2, for
: The cloud core (
of all particles;
) contributes positively to the anisotropy and is responsible for the enhanced peak in
Figure 1, whereas the parts further away either do not contribute (
) or contribute negatively (
) and thus diminish the overall degree of anisotropy. The change observed in the distribution widths in all the velocity dimensions (
Figure 2) is indicative of spatially inhomogeneous flows, as in, e.g., 3D vortexes. As discussed more when presenting the spatiotemporal velocity profiles, the appearance of a nonequilibrium/anisotropic velocity distribution is a signature of a NESS in our system and vortex-like flows can indeed form.
To proceed quantifying the velocity anisotropy (in whole clouds, as in
Figure 1), we use the normalized difference of the root-mean-square (RMS) widths of the fitted Gaussians,
, where
,
are the corresponding Gaussian RMS widths.
Figure 3 shows the velocity anisotropy plotted versus
N for two different
values. Above
atoms, a transition occurs from isotropic to anisotropic velocity, and the anisotropy grows as
N is increased but at a slower rate when
is greater (at
versus
). The transition is correlated with the unstable regime occurring above
atoms, and the slower growth rate—that for a greater
the instability threshold is farther [
25]. Below we present a natural way of extracting the NESS threshold at a fixed
N.
We continue by studying the NESS velocity anisotropy as the unstable regime is approached. In
Figure 4, we display the behavior of the RMS widths
versus
. As can be observed, the widths first decrease and then increase as
gets lower, with the velocity remaining isotropic up to around the point, where they reach their minima. The NESS threshold occurs roughly at this point. We naturally mark it at the detuning beyond which the cloud shape diverges from ellipsoidal, to be discussed when presenting the NESS spatiotemporal characteristics. The decrease of the widths occurs due to the strengthening of the trap’s friction coefficient, whereas the increase is explained by the influence of the many-atom physics. Indeed, in the single-atom limit, we find a continuous decrease of the widths up to a detuning of
, beyond which an increase is observed, in agreement with the Doppler theory [
31]. In
Figure 4, the anisotropy is observed to grow larger as the instability threshold is approached, which is correlated with the cloud shape changes discussed below. Note that we have added the result for the widths for an unstable cloud close to the instability threshold. We observe in this case the velocity distributions to be Gaussian and accompanying greater anisotropy than in the stable regime. We report that the distributions in the unstable regime are generally not Gaussian (close to the threshold or beyond), as different instability regimes exist [
16]. This is, however, beyond the scope of this paper and not discussed further.
In
Figure 5, we plot the NESS threshold versus
(dashed line), superimposed with the instability threshold (solid line; taken from Figure 6 in Ref. [
25]). The NESS threshold is seen to follow roughly the same qualitative behavior as the instability threshold, which decreases linearly with
, with the slope of
per G/cm. The negative
offset from the instability threshold coincides with the cloud changing its shape from ellipsoidal (far-detuned case) to a different non-centrosymmetric configuration. Approaching the instability threshold closer leads to greater deformations, which we discuss when presenting the spatiotemporal profiles of the 3D velocity distributions.
We consider next examining the impact of
I on the NESS velocity anisotropy. In
Figure 6, the results versus
I are displayed for the NESS near the instability threshold (
away), i.e., where the anisotropy is close to being maximal (see the example
Figure 4). The behavior is observed to be highly non-monotonous and separated into three parts: At low
I (
mW/cm
) the anisotropy decreases, then increases at intermediate
I (between 1 and 10 mW/cm
) and, finally, decreases again at large
I (≥10 mW/cm
). This behavior bears some resemblance to the corresponding instability threshold behavior (see the inset), which we note has been qualitatively verified by experimental results (see Figure 3.34 in Ref. [
32]). This resemblance indicates that the anisotropy becomes enhanced the closer the instability threshold is to the resonance. The behaviors are, however, not exactly identical: At small
I the anisotropy is small (
%) compared to the intermediate region anisotropy (growing from
to
%), whereas the thresholds mirror equidistantly in these regions (logarithmic scale), and at large
I the drop in the anisotropy (toward
%) corresponds to the threshold getting plateaued. Due to the NESS phase’s presence being necessary for the MOT’s transition into the unstable regime, this downward trend in the anisotropy size may indicate that the instability mechanism eventually gets broken for larger intensities.
Finally, we discuss the spatiotemporal properties of the 3D velocity distributions of the stable clouds and visually distinguish between the
regular and NESS clouds. In
Figure 7, the velocity arrow plots are displayed for the clouds at different detunings from the instability threshold detuning
(rows), viewed from the side (left column) and the top (right column). Far from the instability threshold (
), the cloud takes on a flattened ellipsoid shape, as expected for a MOT having the magnetic field stronger along one of the axes (the
axis); the velocities are clearly disordered as the particles are subject to the effect of diffusion. Whereas this cloud is outside the NESS region (see
Figure 5), the other two are inside and possess noticeably different non-centrosymmetric configurations. The first one (
) appears triangular (side view) due to being pinched along the beam axes (top view; diagonal beam directions); the velocities start to acquire order as the particles move in a vortex-line manner around the lobes separated by the regions where the pinches occur. Closer to the instability threshold (
), the pinches become more pronounced and so does the order of the velocities.
Figure 8 displays how the superparticle number in each octant is affected for the NESS region: Increasingly unequal amounts are found in adjacent octants when approaching the unstable regime, with the NESS threshold naturally set at the point before this deviation (at
for the considered gradient). We observe the shape changing to be universal as the instability threshold is approached and accompanying growing velocity anisotropy (see the example
Figure 4). This change is indeed due to the growing antagonism between the shadow effect and the rescattering, as the unstable regime is not reached when either of the two effects are removed from the simulations [
26]. On an equivalent note, the stable cloud takes a flattened ellipsoid shape when attenuation (causing the shadow effect and thus fuelling the antagonism) is removed (see Figure 3.45 in Ref. [
32]).
The observation of the inhomogeneous particle-movement is important and supports the argument that a NESS is present. Indeed, as mentioned in
Section 1, such a state can appear in non-periodically driven systems (e.g., our balanced MOT) and is characterized by the existence of nonzero net currents and a nonequilibrium probability distribution (which corresponds to velocity anisotropy [see
Figure 1]). Note that although the particles move inhomogeneously, the envelopes of the NESS region clouds nevertheless exhibit as little fluctuation as the regular stable clouds (see Figure 3.27 in Ref. [
32]), proving that the definition of cloud stability is respected.
As the last remark, there is a possibility of the NESS behavior being quasi-stationary, which otherwise has been observed in different many-body systems with long-range interactions [
33,
34]. In such a case, the NESS clouds would eventually evolve into other stable configurations. The timescales for that would, however, need to exceed 5 s, given this is the longest duration of our simulations. For comparison, such timescales are much longer than the damping period of MOT clouds, taking some tens of milliseconds.