Observation of subdiffusive dynamic scaling in a driven and disordered Bose gas

We explore the dynamics of a tuneable box-trapped Bose gas under strong periodic forcing in the presence of weak disorder. In absence of interparticle interactions, the interplay of the drive and disorder results in an isotropic nonthermal momentum distribution that shows subdiffusive dynamic scaling, with sublinear energy growth and the universal scaling function captured well by a compressed exponential. We explain that this subdiffusion in momentum space can naturally be understood as a random walk in energy space. We also experimentally show that for increasing interaction strength, the gas behavior smoothly crosses over to wave turbulence characterized by a power-law momentum distribution, which opens new possibilities for systematic studies of the interplay of disorder and interactions in driven quantum systems.

Usually interactions are at the heart of the emergent dynamics, but naturally present disorder can also play a crucial role.The study of disorder is a vast field, with highlights including localization and quantum-Hall phenomena in 2D electron gases and quantum wires [26][27][28][29][30], coherent backscattering of acoustic and electromagnetic waves [31,32], and Anderson localization of cold atoms [33,34].Moreover, the interplay of disorder and interactions can result in new phenomena such as many-body localization in lattice systems [35][36][37] and time crystals [38][39][40][41].
In this Letter, we explore the dynamics of a 3D boxtrapped Bose gas strongly driven in presence of weak disorder.In absence of interatomic interactions the gas shows subdiffusive dynamic scaling: its energy grows sublinearly with the drive time t s , approximately as t 0.5 s , and its momentum distribution at different t s is described by a scaling function that is captured well by an isotropic compressed exponential.This behavior is in stark contrast to that expected for a disorder-free noninteracting gas in our cylindrical geometry with forcing along the box axis [see Fig. 1(a)]; in that case the system is effectively 1D and one expects chaotic dynamics with bounded energy growth [42][43][44].Our observations can be explained in terms of a random walk in energy space (see [45] for our detailed theoretical study) and give credence to the proposals that such random walks are a generic feature of thermally isolated driven systems [46][47][48].We also experimentally show, by tuning the interaction strength, that the energy-space random walk observed in the noninteracting limit is continuously connected to wave turbulence characterized by a power-law scaling function [17,[49][50][51].This points to interesting future studies in the regime where the drive, the disorder, and the interactions all play a significant role.
We start with a quasi-pure 39 K Bose-Einstein condensate (BEC) in the lowest hyperfine state, trapped in a cylindrical optical box [52][53][54].The condensate is prepared at a scattering length a = 200 a 0 (where a 0 is the Bohr radius), and we slowly (in 5 s) tune a to zero by tuning the bias magnetic field to 350.4(1) G [55,56].For a noninteracting BEC in our box of length L ≈ 50 µm and radius R ≈ 15 µm, the frequency of the lowest-lying axial excitation is (ignoring disorder) ω z = 3π 2 ℏ/(2mL 2 ) ≈ 2π × 1.5 Hz, where m is the atom mass, while our variable trap depth U D is always larger than 2πℏ × 400 Hz.Weak optical disorder, proportional to the trapping laser power and hence U D , is al-  E x,z , and 1D momentum distributions, ñk (k x,z ), obtained by integrating 2D distributions over k z or k x ; note that E x ≈ E z at t s = 0.The axial E z initially (in < 0.1 s) rises far above E x , but at long t s the two energies are almost equal and grow in unison.The long-time momentum distribution is essentially isotropic [ ñk (k x ) ≈ ñk (k z )], but highly nonthermal.We show ñk (k x,z ) for t s = 3 s, together with the calculated equilibrium distribution for the same energy per particle, E /k B ≈ 23 nK (shaded curve); the experimental distributions show no condensate peak, even though the condensation temperature is 180 nK and the equilibrium distribution has 66% condensed fraction (gray) [57].arXiv:2304.06697v2[cond-mat.quant-gas]17 Dec 2023 Semiclassical picture of a driven and disordered noninteracting Bose gas.(a) The drive mixes only states with k z values up to some k c [42], so in absence of disorder the radial momentum k r (in the k x − k y plane) remains zero, and the growth of a particle's energy is bounded by E c = ℏ 2 k 2 c /(2m).However, disorderinduced elastic scattering distributes energy into the radial modes, and allows the drive to populate states above E c (outside the blue shaded area).(b) This picture implies subdiffusive long-time dynamics, with a sublinear energy growth.Here we consider a particle (red dot) that already has an energy above E c , and k z > k c , so it does not interact with the drive.If elastic scattering (solid circle) reduces its k z to below k c , the particle temporarily interacts with the drive, until another scattering event increases its k z above k c .The interaction with the drive can either increase or decrease the particle's energy, as exemplified by the red arrows.The alternation of scattering and driving events thus results in an energy-space random walk, with a characteristic step size E c .ways present in our holographically created trap [52,59,60], but is typically irrelevant in interacting-gas experiments [61].We inject energy into the system along the box axis z, using a spatially uniform time-varying force [49] of magnitude [42]).After driving the gas for a variable time t s , we probe its momentum distribution using absorption imaging after 50 ms of time-of-flight expansion [58]; this gives line-of-sight integrated 2D distributions n k (k), which we normalize such that n k (k) dk = 1.
In Figs.1(b,c) we illustrate our qualitative observations; here E x,z = ℏ 2 k 2 x,z /(2m)n k (k x , k z ) dk and 1D momentum distributions, ñk (k x,z ), are obtained by integrating n k (k) over k z or k x .Initially the dynamics is essentially 1D, with the drive rapidly (in < 0.1 s) increasing only E z .This is what is expected in absence of disorder, in which case the growth of E z would be bounded [42].However, at long times E x ≈ E z and the energy keeps growing.The long-time momentum distributions along k x and k z are nearly identical, but far from thermal; they show no BEC peak even though the energy per particle is far below the equilibrium condensation value.
In Fig. 2(a) we outline a semiclassical picture of how the interplay of the drive and disorder can lead to isotropic dynamics with unbounded energy growth.The drive mixes only axial modes, with k z up to some k c [42].The disorderinduced elastic scattering transfers energy into the radial modes and allows the drive to increase a particle's energy above E c = ℏ 2 k 2 c /(2m).Alternating scattering and driving events then lead to an unbounded energy growth.
In Fig. 2(b) we explain why at long times this growth is sublinear, corresponding to subdiffusion in momentum space.Once the average energy per particle, E , is significantly larger than E c , most particles have k z > k c and do not interact with the drive.When a particle is occasionally scattered in and out of the k z < k c space, its temporary interaction with the drive can either increase or decrease its energy (red ) ) arrows).This results in an energy-space random walk with a characteristic step size E c [45]: The rate r is energy-dependent, because at any time only a fraction of particles ∝ ℏk c / 2mE interacts with the drive, and because the density of states for elastic scattering is ∝ E .As shown in Ref. [45], this model predicts E ∝ t η s , with 2/5 ≤ η ≤ 1/2 depending on the ratio of the elastic scattering rate and the rate at which the drive mixes k z states.
To isolate and quantify the disorder-induced scattering in our system, we prepare an anisotropic n k with a short t s , stop the drive, and study the subsequent cross-dimensional relaxation for different disorder strengths (∝ U D ).In Fig. 3(a) we show an example of how E z and E x both approach E /3, and in Fig. 3(b) we show how the anisotropy (E z − E x )/E decays for different U D .The initial decay is captured well by an exponential (solid lines), and we use the decay constant Γ d as a measure of the typical scattering rate [62].At long times, the anisotropy decay slows down, which we also observe in simulations of the Schrödinger equation with disorder [42], and attribute to the quantization of states in a finite-size box [63].In Fig. 3(c) we show that, as expected for single-particle scattering, Γ d is independent of the particle density.
In Fig. 3(d) we show the dependence of Γ d on U D .The solid line is based on our numerical simulations [42], which give Γ d ∝ U 2 D , as expected from perturbation theory.We match the data well by setting the rms disorder strength to 2% of U D , and adding a small offset to Γ d .The 2% disorder is compatible with ≈ 1% observed in the bench tests of the optical potentials used for our box trap [60] (see also [59,64]), and also with the measurements with atoms in Ref. [61], where the uniformity of the gas density was confirmed down to energies of a few % of U D .In simulations we assume that (a) Evolution of the atom number N and the radial per-particle energy E r .For t s ≲ 15 s (vertical dashed line), N is essentially constant [65], while E r ∝ t η s with η = 0.46 (solid line); at longer t s some particles have enough energy to leave the trap, so N drops and E r saturates.(b) Evolution of the radial momentum distribution n k (k r ); each curve corresponds to a point in (a), with the same color coding.The dashed line, ∝ k −2 r , is tangent to all the self-similar curves, and k D is the momentum-space trap depth.(c) For scaling exponents α = −0.45 and β = −0.23 (and arbitrarily chosen reference time t 0 = 3 s), the distributions for t s ∈ [1.1, 14.4] s and all k r collapse onto a universal curve.The dotted line shows a compressedexponential fit, f ce ∝ exp − (k r /k 0 ) κ , with κ = 3.0 (see text).The inset shows the results of numerical simulations with the same drive, disorder, and scaling parameters (see text and [42]); for comparison with the experiments, the dotted line is the same as in the main panel.
the disorder is uncorrelated down to a lengthscale of 800 nm, set by the simulation grid, which is comparable with the expected correlation length of experimental disorder, set by the trap-laser wavelength, λ = 532 nm.The small offset in Γ d could arise from trap-shape imperfections or magnetic-field inhomogeneity.
We now turn to the study of the long-time isotropic dynamics for continuous driving (Fig. 4).Here we take images along the drive axis z, which avoids the small effects of the center-of-mass oscillation [42].The distribution in the x + k 2 y ) 1/2 .We normalize 2πk r n k (k r ) dk r = 1 and define E r = E x + E y , so E r ≈ 2E /3 for long t s .
In Fig. 4(a) we plot N (t s ) and E r (t s ) for U s /k B = 7.0 nK, ω s /(2π) = 10 Hz, and Γ d = 8.0 s −1 .At t s ≈ 15 s some atoms reach the momentum-space trap depth k D = 2mU D /ℏ 2 = 3.8 µm −1 , at which point N starts dropping and E r saturates.Until then, N is essentially constant [65] and E r shows power-law growth, E r ∝ t η s , with η = 0.46 (2).In Fig. 4(b) we show the evolution of n k (k r ).For t s ≳ 1 s the distributions are self-similar, with a well-defined front moving towards the UV until it reaches k D ; for t s ≳ 15 s the distribution is essentially stationary.
The fact that n k is self-similar for 1 s ≲ t s ≲ 15 s, while E r grows algebraically, implies dynamic scaling: with β = −η/2, reflecting the subdiffusive energy growth, and α = 2β, reflecting a particle-conserving transport; t 0 is an arbitrary reference time.
In Fig. 4(c) we show that, for t s ∈ [1.1, 14.4] s and all k r , the distributions can be collapsed onto a universal curve, with α = −0.45(2)and β = −0.23(1)[42].The calculations in [45] predict such scaling with the 3D momentum distribution captured by a compressed exponential, ∝ exp − k/k ′ 0 κ 3D , with κ 3D varying between 4 (for η = 1/2) and 5 (for η = 2/5), and k ′ 0 ∝ t 1/κ 3D s [66].We empirically find that n k (k r ), obtained by integrating the 3D distribution along k z , is captured by a normalized compressed exponential [dotted line in Fig. 4(c)] [67].As shown in the inset of Fig. 4(c), we reproduce our observations in numerical simulations; here we show the results of simulations for t s ∈ {2.9 − 18} s, obtained with the same U s , ω s , and Γ d , and collapsed with the same α, β, and t 0 as in the experiments.
We conclude by pointing to an interesting question for future study -what happens if the drive, the disorder, and the interactions all play a significant role?In the noninteracting (a = 0) dynamics observed here, the rate at which energy is absorbed from the drive decays as t η-1 s .On the other hand, in interacting wave-turbulent cascades [17,[49][50][51] this rate is constant and the turbulent steady state is characterized by , with γ = 3.3 (3).The two types of dynamics are qualitatively different, but are continuously connected by tuning a, as we illustrate in Fig. 5(a) for one set of drive and disorder parameters, and fixed gas density (∝ N ).This crossover should be controlled by some dimensionless parameter(s), but from the drive, disorder, and interaction properties, one can construct many candidates.Moreover, further qualitatively different outcomes are possible: for strongenough interactions, thermalization should be the fastest process, while for strong-enough disorder, localization effects should prevail.Constructing the full dynamical phase diagram for the driven and disordered interacting gas is thus a fascinating challenge.As a first step in this direction, in Fig. 5(b) we show that for our parameters the crossover from an energy-space random walk to turbulence depends on the product N a, which suggests that it can be captured within the mean-field Gross-Pitaevskii framework.
In summary, we have observed subdiffusive dynamic scaling in a noninteracting Bose gas driven far from equilibrium in the presence of weak disorder, which we explain in terms of an energy-space random walk.The tuneability of our system opens the possibility to study the interplay of (with γ = 3.4) for 0.6 µm −1 < k r < 2.5 µm −1 , with c and n 0 as free parameters.For various N and a, the crossover from c = 1 (subdiffusion) to c = 0 (turbulence) depends only on the product N a.
the drive, disorder, and interactions in regimes where they all play a significant role, which we illustrate by showing how the energy-space random walk crosses over to turbulentcascade dynamics.Our far-from-equilibrium states with low and tuneable energy per particle also provide a novel starting point for studies of equilibration in closed quantum systems [24,25,68].

SUPPLEMENTAL MATERIAL I. NUMERICAL SIMULATIONS
We describe our system using the time-dependent Schrödinger equation for a particle trapped in a cylindrical box of radius R = 15 µm and length L = 50 µm and driven with a potential V s (r, t ) = U s z/L × cos(ω s t s ).We model the disorder present in our system as V d (r i ) = 3V 0 X i , where r i are the lattice points in our simulation grid, X i are independent random numbers sampled uniformly between −1 and 1, and V 0 is the rms disorder strength.Our simulations always begin with the system in the ground state at t s = 0.
For V 0 = 0, the problem is separable into radial and axial directions, since the drive is oriented along the box axis z.We solve the resultant 1D problem by expanding the wavefunction in the box eigenbasis (using the lowest 80 states) and solving the coupled differential equations for the expansion coefficients.In Fig. S1 we plot the calculated momentum distribution ñk (k z ) and energy E z versus excitation time t s for U s /k B = 10.5 nK and ω s /(2π) = 10 Hz (cf.

II. FREQUENCY RESPONSE
Here we experimentally study the frequency response of our noninteracting box-trapped gas.In Fig. S2(a), we show for a fixed Γ d and three different U s how E r depends on ω s for a t s = 4 s sinusoidal drive.At low frequency E r increases roughly linearly with ω s but above some U s -dependent ω peak the system stops responding strongly.
As illustrated in Fig. S2(b), we observe qualitatively different behavior below and above ω peak .Here we plot the evolution of ñk (k z ) around t * s = 4 s, for U s /k B = 7.0 nK.For ω s /(2π) = 10 Hz (below ω peak ) the distribution shows no condensate peak, while for ω s /(2π) = 25 Hz (above ω peak ) the lowest mode remains macroscopically occupied.
In the main paper we always use ω s < ω peak , which ensures that our drive uniformly mixes states with k z < k c .

III. EXTRACTION OF THE SCALING EXPONENTS
Here we detail our extraction of the scaling exponents (α, β) from n k (k r , t s ).We start by fitting n k (k r ) for different t s with f ce = A 0 exp[−(k r /k 0 ) κ ], with A 0 , k 0 and κ as free parameters, as shown in Fig. S3(a) for three experimental curves from Fig. 4(b) in the main paper.As shown in Fig. S3(b), we get essentially constant κ for some t s range (1.1 − 14.4 s); this is the scaling range in which the n k curves are self-similar.We then fix κ to its average value in the scaling range and refit the curves to obtain k 0 and A 0 .Finally, as shown in Fig. S3(c), we extract α and β using power-law fits of the form k 0 ∝ t −β s and A 0 ∝ t α s (solid lines).

FIG. 1 .
FIG. 1. Noninteracting box-trapped Bose gas driven far from equilibrium.(a) Box geometry and driving force, F = (U s /L) cos(ω s t s ) ẑ, where L ≈ 50 µm is the box length.(b) Time-of-flight images, giving 2D momentum distributions n k (k x , k z ), for N = 3.3 × 10 5 atoms in a box of depth U D /k B = 90 nK, driven at ω s /(2π) = 10 Hz with U s /k B = 10.5 nK.The scale bar shows 1 µm −1 and the optical density (OD) saturates at 3. (c) Energies,Ex,z , and 1D momentum distributions, ñk (k x,z ), obtained by integrating 2D distributions over k z or k x ; note that E x ≈ E z at t s = 0.The axial E z initially (in < 0.1 s) rises far above E x , but at long t s the two energies are almost equal and grow in unison.The long-time momentum distribution is essentially isotropic [ ñk (k x ) ≈ ñk (k z )], but highly nonthermal.We show ñk (k x,z ) for t s = 3 s, together with the calculated equilibrium distribution for the same energy per particle, E /k B ≈ 23 nK (shaded curve); the experimental distributions show no condensate peak, even though the condensation temperature is 180 nK and the equilibrium distribution has 66% condensed fraction (gray)[57].
FIG. 3. Disorder-induced cross-dimensional coupling.(a) Here we stop the drive with U s /k B = 10.5 nK and ω s /(2π) = 10 Hz after t s = 0.1 s [see Fig. 1(b)] and then hold the gas for a variable time t hold in a trap with U D /k B = 90 nK.The energies E z and E x both relax towards E /3 = (E z + 2E x )/3 (by symmetry E x = E y ).(b) The decay of the anisotropy (E z − E x )/E is initially exponential, ∝ exp(−Γ d t hold ) (solid lines); Γ d grows with disorder strength (∝ U D ) and is independent of the initial anisotropy.(c) As expected for single-particle scattering, Γ d is independent of the gas density (∝ N ); here U D /k B = 90 nK.(d) Dependence of Γ d on U D .Up to a small offset of 1 s −1 , the data are captured by our numerical simulations with rms disorder equal to 2% of U D (solid line).

FIG. 5 .
FIG. 5. Crossover from subdiffusion to turbulence, for U s /k B = 7.0 nK, ω s /(2π) = 10 Hz, Γ d = 15 s −1 , and t s = 1 s.(a) n k for N = 10 5 atoms and different scattering lengths a.For subdiffusion, n k is captured by a compressed exponential f ce .For a turbulent cascade, n k ∝ k −γ+1 r (dashed line) for k r ≳ 1/ξ, where ξ is the healing length; here ξ = 0.6 µm −1 for 100 a 0 .(b) We quantify the crossover between the two regimes by fitting n k = c f ce (k r ) + n 0 k −γ+1 r Fig. 1 in the main text).The drive only couples k z modes up to a characteristic momentum k c [Fig.S1(a)], and the average E z saturates (with large fluctuations) after a rapid initial increase [Fig.S1(b)].To model the 3D dynamics observed in experiments, we solve the full 3D Schrödinger equation using a pseudospectral method with fourth-order Runge-Kutta time evolution with a 10 µs timestep and a numerical grid of size 100 × 100 × 100 µm 3 discretized with ζ = 0.8 µm resolution, which is on the order of the trapping laser wavelength, λ = 532 nm.This results in Fourier components of the disorder up to π/ζ ≈ 4 µm −1 , similar to the largest k in the experiment, set by k D .Simulating the experimental protocol from Fig. 3 of the main text, we find that V 0 /k B ≈ 2 nK reproduces the experimental Γ d = 8.0 s −1 for U D /k B = 90 nK.In these simulations we also observe, as in the experiments, that the decay of anisotropy slows down at long times.For the simulations in the inset of Fig. 4(c) of the main text, we have also extracted the scaling exponents as in the experiment (see Section III) and found η = 0.46(1), κ = 3.2(1) [κ 3D ≈ 4], α = −0.46(1),and β = −0.24(1).

6 FIG
FIG. S1.Simulations of a strongly driven Bose gas in a disorderfree box.(a) The 1D ñk (k z ) versus t s for U s /k B = 10.5 nK and ω s /(2π) = 10 Hz.The dashed lines indicate k c .(b) The corresponding bounded energy growth.The energy E z shows large fluctuations around an average value that quickly saturates (see inset).

ω
FIG. S2.Frequency response of a noninteracting Bose gas.(a) E r versus ω s for a t s = 4 s sinusoidal drive with different U s (see legend) and Γ d = 8.0 s −1 .We show data taken with two different initial atom numbers (circles: 2.8 × 10 5 , diamonds: 1.1 × 10 5 ).The inset shows a typical time-of-flight image taken along ẑ for ω s < ω peak .(b) Evolution of the axial momentum distribution around t * s = 4 s for ω s < ω peak (left) and ω s > ω peak (right).