Energy-space random walk in a driven disordered Bose gas

Motivated by the experimental observation [1] that driving a non-interacting Bose gas in a 3D box with weak disorder leads to power-law energy growth, $E \propto t^{\eta}$ with $\eta=0.46(2)$, and compressed-exponential momentum distributions that show dynamic scaling, we perform systematic numerical and analytical studies of this system. Schr\"odinger-equation simulations reveal a crossover from $\eta \approx 0.5$ to $\eta \approx 0.4$ with increasing disorder strength, hinting at the existence of two different dynamical regimes. We present a semi-classical model that captures the simulation results and allows an understanding of the dynamics in terms of an energy-space random walk, from which a crossover from $E \propto t^{1/2}$ to $E \propto t^{2/5}$ scaling is analytically obtained. The two limits correspond to the random walk being limited by the rate of the elastic disorder-induced scattering or the rate at which the drive can change the system's energy. Our results provide the theoretical foundation for further experiments.


I. INTRODUCTION
The emergence of simple and universal behaviors insensitive to system parameters and past trajectories is one of the most fascinating aspects of the physics of complex systems.Although the theory of universal behaviors was traditionally developed for equilibrium critical phenomena [2], recent experimental and theoretical studies have extended these ideas to a wide range of far-from-equilibrium systems [3][4][5][6][7][8][9][10][11].
While interactions are usually central to the universal dynamics, in a recent experiment [1], we demonstrate that in absence of interactions, an interplay between drive and disorder can also lead to universal behavior.This system, with a power-law energy growth [E ∝ t η with η = 0.46 (2)] and self-similar momentum distributions well characterised by a compressed exponential, shows qualitatively different behavior from its interacting counterpart.
In Ref. [1], these observations are reproduced with Schrödinger-equation simulations and qualitatively explained by a semi-classical model.In this paper, we formalize our theoretical results.First, we extend the Schrödingerequation simulations to a wider parameter range and observe a crossover from η ≈ 0.5 to η ≈ 0.4 with increasing disorder strength (Section II), which hints at the existence of two distinct dynamical regimes.We then present the semi-classical model (Section III) that captures the simulation results and allows an understanding of the dynamics in terms of an energyspace random walk.This in turn leads to a simple energyspace drift-diffusion equation (Section IV) that reproduces the crossover between the two regimes, and analytic predictions of E ∝ t 1/2 and E ∝ t 2/5 that emerge in the limits where the random walk is limited by the rate of disorder-induced scattering or the rate at which the drive can change the system's energy.Our results offer a new example of a dynamical system undergoing energy-space drift-diffusion [27][28][29][30][31] and provide the theoretical foundation for further experimental studies.

II. SCHRÖDINGER-EQUATION SIMULATIONS
The non-interacting dynamics in Ref. [1] can be described by the Schrödinger equation where m is the particle mass, V box is the clean trapping potential, V D is the disorder, L is the box length along the drivingforce direction z, and U /L is the amplitude of the driving force.Here, we model the trap as a cubic box [Fig.1(a)] of infinite depth [32], and the disorder V D is chosen to be an uncorrelated (zero-mean) Gaussian random potential.The choice of an uncorrelated potential is sensible because the correlation length of V D in an optical trap, which is on the order of the laser wavelength λ, is small compared to the atomic de-Broglie wavelength in the experiment.The strength of the random potential V D is characterized by its r.m.s.value σ.
In Fig. 1(b), we illustrate the evolution of the momentum distribution, n k (k) = |ψ(k)| 2 , for one choice of parameters {U , ω, σ}.As in the experiment [1], the drive rapidly increases the momentum spread along z, and cross-dimensional coupling due to V D causes energy to leak into the transverse directions.At long times, n k (k) is nearly isotropic and gradually broadens [Fig.1(c)].
In agreement with the observations in Ref. [1], the energy growth is well described by a power-law, E (t ) ∝ t η , as shown in Fig. 1(d).The (nearly-)isotropic momentum distributions n k (k, t ) at different t are self-similar, with where t ref is an arbitrary reference time, β = −η/2, and α = 3β corresponding to particle-conserving transport [33].This self-similarity is illustrated in Fig. 1(e) for different parameters {U , ω, σ}; for each simulation, the distributions at different t collapse onto a single curve when rescaled according to Eq. (2).The collapsed curves are well described by compressed exponentials [black lines in Fig. 1

(e)] of the form
with exponent κ and momentum scale k s ∝ E .
In Ref. [1], only a relatively narrow range of η and κ was observed.Here, by extending the range of disorder 10 -8 ) t/t 0 1 3 10 30 strengths σ, we observe a crossover from η ≈ 0.5 and κ ≈ 4 to η ≈ 0.4 and κ ≈ 5 [Fig.1(f)].This hints at the existence of two distinct dynamical regimes, corresponding to weak and strong disorder.Analytically understanding the emergence of these two regimes is the goal of the subsequent sections.

III. SEMI-CLASSICAL MODEL
The key ideas used to develop our model for the interplay of the drive and disorder are illustrated in Fig. 2. First, we note that in the absence of disorder, strongly driving the gas along a separable axis of the trap leads to 1D chaotic dynamics with bounded energy growth [1,34,35].This is illustrated by the (disorder-free) 1D Schrödinger-equation simulations shown in Fig. 2(a), where we initialize the system in different sine-basis states |ψ(t =0)〉 = |k z,0 〉 of the box (red dots).For small k z,0 , the strong drive mixes the k z states only up to a cutoff k c (horizontal dashed line), while for large k z,0 , the drive only weakly perturbs the system.However, the presence of disorder significantly modifies the picture in 3D.While the drive can only increase k z up to about k c , the disorder can scatter particles to equal-energy states with lower k z , where the drive can again increase their energy [see Fig. 2(b)].This cooperative process is the key to the unbounded energy growth.
To model this process, we propose the following semi-classical kinetic equation where Θ(k) is the Heaviside function, and s and f , respectively, characterize the rate of the elastic disorder-induced scattering and the rate at which the drive can change the system's energy.The first line of Eq. ( 4) describes the elastic disorderinduced scattering.A perturbative treatment using Fermi's golden rule gives the scattering rate from a state |k〉 as For uncorrelated V D , after ensemble-averaging where s ∝ σ 2 (see Appendix A 1), and the factor of k arises from the 3D density of states.This leads to the −s |k|n k (k, t ) term in Eq. ( 4) for the population out-flux from state |k〉.The integral term in the first line of Eq. ( 4) describes the in-flux to state |k〉; since the out-flux from each state contributes an equal in-flux to every other state on the same k-shell, the total in-flux to state |k〉 due to scattering is given by the out-flux averaged over the shell.c /(2m), and the time t by 1/ f .The colored lines are power-law fits used to extract η.The gray dashed lines show analytic predictions for sk c / f = 10 2 and 10 −3 using Eq. ( 18) and Eq. ( 23), respectively.(b) Momentum distributions for the simulations shown in (a), dynamically scaled according to Eq. ( 2), using the extracted η and an arbitrary reference time t ref = 5 × 10 3 / f .The solid lines show compressed-exponential fits used to extract κ.(c) Extracted η and κ as a function sk c / f , showing a similar crossover between the two dynamical regimes as seen in Fig. 1(f).
The second line of Eq. ( 4) heuristically models the driving process.While the chaotic 1D dynamics is not amenable to an exact treatment, the simulation results in Fig. 2(a) inspire a simple model, where the drive randomly mixes k z states up to k c at a phenomenological rate f , without affecting states with k z > k c .While we also treat k c phenomenologically, its value can be estimated from the time-averaged energy of the driven 1D system (see Appendix A 2).
Before analytically studying this model, we validate it through stochastic numerical simulations of Eq. ( 4) for different values of the dimensionless parameter sk c / f , which sets the ratio of the elastic scattering rate to the rate at which the drive can change the system's energy.As shown in Fig. 3, our model, despite its simplicity, captures all the key features seen in the Schrödinger-simulation results in Fig. 1.

IV. ANALYTIC ANALYSIS A. Qualitative ideas
To solve our model analytically, we switch from momentum space, where Eq. ( 4) describes a highly non-local jump process, to energy space, where the process is quasi-local.In energy-space, the trajectory of a particle can be described by a sequence of events of the form ...
where S and D refer to individual scattering and driving events, and Z S (E ) and Z B ( Ē ) label the state of the particle, where the subscripts stand for shell and band, respectively [see Fig. 4 (a)].Immediately after a scattering event S, the particle is randomly distributed on a k-shell of energy E , so its state can be labeled Z S (E ).Similarly, immediately after a driving event D, the particle is randomly distributed on a cylindrical band with radius is known; we label such a state Z B ( Ē ).
Note that successive S or D events do not change the state Z S or Z B , but the state changes when S and D alternate, as illustrated in Figs.space random walk with d dt where r (E ) is the (generally energy-dependent) rate at which S and D alternate.
In the strong-scattering regime, sk ≫ f , the rate r (E ) is limited by the occurrence of D events, so r (E ) ∝ (k c /k) f , where the k c /k factor arises because the particle takes part in the random walk only when k z < k c .Since k ∝ 〈E 2 〉 1/4 , Eq. ( 7) implies in agreement with both the Schrödinger-equation simulations and the stochastic simulations of our model.On the other hand, in the strong-driving regime, f ≫ sk, the rate r (E ) is limited by the occurrence of S events and given by r (E ) ∝ (k c /k) sk = sk c .As the suppression from k c /k is cancelled by the density of states factor k in the scattering rate, Eq. ( 7) implies which is also consistent with both the Schrödinger-equation simulations and the stochastic simulations of our model.

B. Energy-space drift-diffusion equation
We now formalize the ideas from Eq. ( 7) and derive an energy-space drift-diffusion equation valid for all values of sk c / f .Since the particle changes state only when S and D alternate, we can more succinctly describe its trajectory as where S and D, respectively, stand for D...DS and S...SD.Heuristically, the energy change may be described by where the 'velocity' v E and 'random force' ζ(t ), respectively, lead to energy drift and diffusion.Denoting 〈T S,D 〉 as the mean waiting time for S and D to happen, and µ S,D and σ 2 S,D the energy drift and the energy-variance production in each step, we have As derived in Appendix B, we have The expressions for 〈T S + T D 〉, µ D , and σ S,D agree with the qualitative discussion in Section IV A. The non-zero µ S arises because a particle in Z B ( Ē ) is more likely to be scattered out of the band at higher k z , due to the k dependence of the scattering rate.
Combining Eqs. ( 10) -( 12), we can write a drift-diffusion equation [36] for the energy distribution P (E , t ): The formal derivation of Eq. ( 13), using the theory of continuous-time random walks [37], is given in Appendix B. This equation satisfies the non-equilibrium fluctuationdissipation relation proposed in Ref. [29], as discussed in Appendix C. The solid line shows the analytic prediction K s = (πk c /k 0 ) 5 /45 [calculated from Eq. ( 20)] with no free parameters.

C. Limiting regimes
For sk ≫ f , Eq. ( 13) reduces to with diffusion constant Following Ref. [38], Eq. ( 14) can be shown to support selfsimilar solutions of the form The corresponding momentum distribution, is a compressed exponential [see Eq. ( 3)] with κ = 5, and the energy growth is a power law with η = 2/5 in agreement with Eq. ( 8).
For sk ≪ f , Eq. ( 13) reduces to with diffusion constant The self-similar solution supported by Eq. ( 19) is The corresponding momentum distribution, is a compressed exponential with κ = 4, and the energy growth is a power law with η = 1/2 in agreement with Eq. ( 9).
To quantitatively verify Eqs. ( 18) and ( 23), in Fig. 3(a), we compare them to our stochastic-simulation results for sk c / f = 10 2 and 10 −3 and observe good agreement.
In the low-disorder limit (sk ≪ f ), we can also directly compare our analytically predicted energy-diffusion coefficient D s [Eq.(20)] with the Schrödinger-equation simulations, because D s depends only on s and k c , both of which can be obtained from the input parameters {U , ω, σ} (such a comparison is not possible for D d because we cannot calculate f ).For each simulation in the low-disorder limit [39], we fit the E (t ) curve [such as shown in Fig. 1(d)] to Eq. ( 23) and extract D s .In Fig. 5(a), for fixed ω and various U , we plot the extracted D s versus s calculated from σ and observe the linear behavior predicted in Eq. (20).Then, in Fig. 5(b), for several ω and a range of U , we show that the fitted constants of proportionality between D s and s [the slopes of lines in Fig. 5(a)] agree with Eq. (20).

D. General solution
After analysing the two limits, we now examine the general solution to Eq. ( 13).We can remove the parameters {s, k c , f } from the equation by introducing the dimensionless quantities with This transforms Eq. ( 13) to  E/E sys FIG. 6. Results of extended stochastic simulations of the semi-classical model with different sk c / f .We show the evolution of the energy E (left) and of the compressed-exponential exponent κ for the instantaneous shape of n k (right).We perform simulations with different sk c / f , indicated by the color, which merge into common curves when t and E are scaled using units given in Eq. ( 25).In both plots, the dashed lines show the analytic predictions in the two limiting regimes, and the solid black lines show numerical predictions using Eq. ( 13).
The small deviations of stochastic simulations from the solid lines (not visible in the E plot) arise due to initial transients in each simulation.
with P ′ (E ′ , t ′ ) = E sys P (E , t ).This shows that, under appropriate scaling, solutions to Eq. ( 13) follow a universal E − t trajectory.We illustrate this in Fig. 6.In principle, at very long times, one should always observe E ∝ t 2/5 and κ = 5.In terms of system parameters, we classify the system as low-disorder if sk c ≪ f and highdisorder if sk c ≫ f , but the dynamics is actually controlled by the ratio sk/ f , which increases as the energy grows.Thus, a low-disorder system at long times is mathematically identical to a high-disorder one at short times.However, note that the crossover between the two regimes occurs over an enormous timescale, so any realistic experiment will sample a small region of the universal trajectory, with the energy growth well fitted by a power law and an essentially constant κ.Also note that for a strongly disordered system, t sys is very short (sk c / f ≫ 1, so t sys ≪ 1/ f ), so by the time any significant energy is pumped into the system, t /t sys is already large.

V. CONCLUSION & OUTLOOK
In conclusion, we have developed a semi-classical model for a driven non-interacting box-trapped Bose gas in the presence of uncorrelated disorder.The dynamics at the heart of this model can be understood in terms of an energy-space random walk, and the resulting analytic predictions reproduce the key features seen both in the experiment of Ref. [1] and in Schrödinger-equation simulations.
Our work points to several future directions.First, it would be interesting to experimentally explore the dynamics beyond the weak-disorder regime.This could also lead to further theoretical questions, as the scattering in our model is treated within first-order perturbation theory, which does not hold for arbitrarily strong disorder; for example, the onset of Anderson localization [40] may lead to additional dynamical regimes.Second, an analogous study in 2D may reveal even richer physics.Taking our model at face value, we would expect η = 2/5 across all parameter regimes in 2D, because there is no density-of-states enhancement factor k in the scattering rate, so we always have r ∝ 1/k and E ∝ t 2/5 .However, this may be inaccurate due to the more prominent role of fluctuations in 2D.For example, our treatment of the scattering rate relies on ensemble-averaging.While this is a good approximation in 3D, its validity in 2D is not obvious, as far fewer states are involved in the scattering process.This poses interesting questions both experimentally and theoretically.In this section, we present the calculation for Γ s (k) in Eq. ( 5) and derive an explicit expression for s in Eq. ( 6) for a cubic box of volume L 3 in the presence of a disorder potential V D (r).The unperturbed basis states of the box are |k〉 states of the form |k〉 = 8 The ensemble-averaged matrix element |〈k|V D (r)|k ′ 〉| 2 is explicitly where By substituting the Fourier representation into Eq.(A2), we get The {x, y, z} integrals are separable, with 1D integrals of the form For large momenta, the sinc functions in the above equation have almost no overlap, so we can drop the interference terms: and using lim L→∞ sinc 2 (qL/2) = (2π/L) δ(q), we get Substituting this expression back into Eq.(A5), we get where λ i , λ ′ i can take on values ±1.Intuitively, {k ′′ } is the set of 64 k-vectors that connect the plane-wave components in |k〉 to those in |k ′ 〉.
In numerical simulations, we discretize real space into (N − 1) 3 grid points, which requires some changes to the above equations.First, the decomposition of the C (r 1 − r 2 ) is written in terms of a Fourier series rather than a Fourier transform, where the summation is performed over the first Brillouin zone of the grid.Second, the delta functions δ(k ′′ − q) become G δ k ′′ −q,G where {G} is the set of reciprocal lattice vectors.This change introduces unphysical Umklapp scattering processes.However, for the uncorrelated disorder potential V D (r) with zero mean and variance σ 2 , the correlation function is C (r) = δ r,0 σ 2 , so C (k ′′ ) is constant and the Umklapp scattering processes do not affect the physics.Eq. (A9) then gives Note that the factor of 1/64 in Eq. (A9) disappears because all 64 k ′′ contribute to the sum.Finally, substituting Eq. (A11) into Eq.( 5), we get where the sum has been approximated by an integral and we can read off the scattering parameter s in Eq. ( 6) as (A13)

The cutoff momentum k c
Here we detail our method to determine k c .In our semi-classical model, we have assumed that the drive randomly mixes states with k z < k c at a rate f .Therefore, in 1D, a state initialized with k z,0 < k c will reach a momentum distribution n 1D k that is (on average) uniform below k c and zero above it.For k z,0 < k c , the mean energy of the driven system is thus Therefore, k c may be estimated by computing 〈E 〉 from 1D Schrödinger-equation simulations for a driven particle in a disorderfree box.In Fig. A1(a), we show 〈E 〉 for U = 1500 E 0 and ω = 75 E 0 /ℏ, starting from different k z,0 .For low k z,0 , 〈E 〉 is indeed essentially independent of k z,0 [see also Fig. 2  The waiting time T D is a continuous random variable.We begin by calculating its mean, 〈T D 〉, before calculating its full distribution.Suppose that at t = 0, the particle has just been scattered into Z S (E ) with unknown position on the k-shell.Then, with probability (1 − k c /k), it sits in the upper part of the shell ('upper shell'), where k z > k c .In this case, the next event has to be S, and the waiting time for it is Exp[sk] [an exponential distribution with time constant 1/(sk)].After the S event, the clock is reset because the particle is still in the same statistical state Z S .Therefore, the additional waiting time before D happens is again T D .Hence, where T D|u is the conditional waiting time till D if the particle is initially in the upper shell.On the other hand, if the particle is initially in the lower shell, where k z < k c , the next event could be either S or D, and the waiting time for this event is Exp[sk + f ].With probability f /(sk + f ), the event is D, and the particle is driven out of the shell.Otherwise, the event is S, the clock is reset, and the additional waiting time before D happens is again T D .This can be written as where T D|ℓ,S and T D|ℓ,D are the conditional waiting times till D if the particle is initially in the lower shell and the first event after t = 0 is S or D, respectively.We can thus write an equation for 〈T D 〉: The solution for 〈T D 〉 is For k ≫ k c , to leading order in k c , Generalizing the analysis above, we can also write an equation for the distribution function φ D (t ) for T D : where g (t ; λ) denotes the exponential distribution function with time constant 1/λ.By taking the Laplace transform the integral equation Eq. (B7) can be reduced to an algebraic one, Its solution is The exact expression for φ D (t ), obtained from the inverse Laplace transform of φD (u), is complicated, but it can be well approximated by an exponential distribution.This can be seen from the above equation, where for small u (corresponding to large t and large k), we have [using 〈T D 〉 in Eq. (B6)] When Z S (E ) D − → Z B ( Ē ) happens, the new energy Ē is a random variable with distribution G( Ē |E ).Irrespective of T D , it is equally likely for the particle to be driven from any point on the lower shell.This means that c , k because the particle is in the lower shell.After some manipulation, we get c /6 .From this distribution, we calculate (B14) 2. Calculations for S = D...DS

a. Calculation of T S
The waiting time T S is calculated along the same lines as T D above, but we need to also average over the initial k z in the band (with fixed k ⊥ ) since the scattering rate sk depends on k z via k = k 2 ⊥ + k 2 z .This gives where the first term in the integral comes from the mean waiting time till the first event (either D or S) after t = 0, and the second term comes from the additional time needed if the first event is D. Solving the equation, we get (B16) Note that the integrals in the above equation cannot be evaluated analytically, but by treating k c and k z as small parameters, we can perform Taylor expansions and obtain the simple result above.It is also possible to obtain the distribution of T S using the Laplace transform, Therefore, the distribution for T S is also exponential to leading order.When Z B ( Ē ) S − → Z S (E ) happens, the new energy E is a random variable.The exact distribution for E is tricky to calculate because it is correlated with T S .In particular, if the particle is scattered when it has a higher k z (and hence a higher k), E will be larger, and T S is likely to have been shorter due to the k-dependence of the scattering rate.However, for an approximate calculation, we ignore this correlation and calculate the distribution of E irrespective of T S .The error of this approximation is a higher-order term.
First, let us calculate the probability p(k z | Ē ) that the particle is scattered out at k z , and hence into the shell with E = (k 2 ⊥ + k 2 z )/2.If the particle is at k ′ z at t = 0, then, with probability sk ′ /(sk ′ + f ), the first event after t = 0 is S, and we get a contribution to p(k z | Ē ) only if k ′ z = k z ; otherwise, with probability f /(sk ′ + f ), the first event after t = 0 is D, and the probability that the particle leaves at k z in some future S event is p(k z | Ē ).Averaging over k ′ z , we get ) . (B20)

The drift-diffusion equation
We can now derive Eq. ( 13) by generalizing the approach of Ref. [37].First, we denote the probabilities of the particle being in states Z S (E ) and Z B (E ) by P 1 (E , t ) and P 2 (E , t ), respectively.We express the rate of change of P 1 (E , t ) and P 2 (E , t ) as where J + 1,2 and J − 1,2 are the in-and out-fluxes, respectively.Since the state of the particle has to alternate between Z S and Z B between steps of the energy-space random walk, we have where G 1,2 (E |E ′ ) are the energy transition probabilities for D and S, respectively.The out-fluxes J − 1,2 can be written as where φ 1,2 (E , t ) = φ D,S (E , t ) are the waiting-time distributions.The first terms in Eqs.(B23) represent the fluxes contributed by particles originally in the states 1, 2 at t = 0, and the second terms represent the fluxes contributed by particles that enter the states at t ′ and leave at t .Using Eq. (B21), we eliminate J + 1,2 from Eq. (B23) and get These integral equations can be solved using the Laplace transform, which gives with the memory kernels given by M1,2 (E , u) = φ1,2 (E , u) 1 − φ1,2 (E , u) .
Appendix C: Non-equilibrium fluctuation-dissipation relation In cases where the dynamics of a periodically driven and thermally isolated system obeys an energy-space drift-diffusion equation of the form (C1)

)FIG. 1 .
FIG. 1. Schrödinger-equation simulations of a driven non-interacting Bose gas in a box with disorder.(a) Illustration of the simulation geometry.For the box of size L, the natural units of momentum, energy, and time are, respectively, ℏk 0 = ℏπ/L, E 0 = ℏ 2 / mL 2 , and t 0 = ℏ/E 0 .(b) Snapshots of the projected momentum distribution ñk (k x , k z ) for drive parameters U = 1500 E 0 and ω = 75 E 0 /ℏ, and disorder strength σ = 〈V 2 D 〉 1/2 = 750 E 0 (the simulation grid is of size 127×127×127, which leads to a UV cutoff of 127k 0 ).For comparison to Ref. [1], using the 39 K atom mass m = 6.5 × 10 −26 kg and L = 50 µm, the simulation parameters here correspond to U /k B = 7.4 nK, ω/(2π) = 9.5 Hz, and σ/k B = 3.7 nK.The scale bar corresponds to 20 k 0 .(c) Evolution of the (spherically averaged) momentum distribution n k for parameters as in (b) and t /t 0 ∈ [0.85, 55.8].(d) Energy-growth dynamics for U = 1500 E 0 , ω = 75 E 0 /ℏ, and various σ.The solid lines show power-law fits used to extract the energy-time scaling exponent η.(e) For each σ value in (d), momentum distributions at different t (such as shown in (c) for σ = 750 E 0 ) collapse onto a single curve when dynamically scaled according to Eq. (2), with arbitrarily chosen t ref = 10 t 0 .The solid lines show fits according to Eq. (3), used to extract the compressed-exponential exponents κ.(f) Extracted η and κ as a function of σ for fixed U = 1500 E 0 and ω = 75 E 0 /ℏ.

FIG. 2 .
FIG. 2. Key ideas underpinning our semi-classical model.(a) Numerical simulations of the (disorder-free) 1D Schrödinger equation for U = 1500 E 0 and ω = 75 E 0 /ℏ, starting from different initial sine-basis states (red dots).The density plots show the 1D momentum distribution n 1D k (k z , t ).The horizontal dashed line indicates the cutoff momentum k c (see text and Appendix A 2).(b) The unbounded energy growth process in our model.The dots indicate the sine-basis eigenstates |k〉 of the disorder-free box.At t = 0, the particles start in the ground state (red dot), and their k z may be increased by the drive (blue arrows) up to k c .The disorder-induced scattering (orange arrows) moves the particles along equal-energy shells and provides opportunities for the drive to further pump energy into the system.

FIG. 3 .
FIG. 3. Stochastic-simulation results for our semi-classical model.(a) Energy growth over time for different sk c / f .We normalize the energy E by E c = ℏ 2 k 2c /(2m), and the time t by 1/ f .The colored lines are power-law fits used to extract η.The gray dashed lines show analytic predictions for sk c / f = 10 2 and 10 −3 using Eq.(18) and Eq.(23), respectively.(b) Momentum distributions for the simulations shown in (a), dynamically scaled according to Eq. (2), using the extracted η and an arbitrary reference time t ref = 5 × 10 3 / f .The solid lines show compressed-exponential fits used to extract κ.(c) Extracted η and κ as a function sk c / f , showing a similar crossover between the two dynamical regimes as seen in Fig.1(f).

FIG. 4 .
FIG.4.Qualitative ideas underpinning the energy-space random-walk picture.(a) Statistical states of the particle in energy space: Z S (E ) (orange) labels the state of a particle randomly distributed on a spherical shell with definite k and E , while Z B ( Ē ) (blue) labels the state of a particle randomly distributed in a cylindrical band with definite k ⊥ , random k z ∈ [0, k c ], and mean energy Ē .(b) Effects of a driving event after scattering events.When the driving event happens, the particle originally in Z S (E ) (orange line) may be driven into a range of Z B ( Ē ) states (blue shaded area) with a range (∆k ⊥ ) of k ⊥ , and hence a range of Ē .(c) Effects of a scattering event after driving events.When the scattering event happens, the particle originally in state Z B ( Ē ) (vertical blue line) may be scattered into Z S (E ) states (orange shaded area) with a range (∆k) of k, and hence a range of E .
Appendix A: Calculation of the semi-classical model parameters 1.The scattering rate FIG. A1.Extraction of k c from 1D Schrödinger-equation simulations without disorder.(a) Time-averaged energy, 〈E 〉, for the system initialized in different states k z,0 ; see also Fig. 2(a).We normalise k z,0 by k 0 = π/L and 〈E 〉 by E 0 = ℏ 2 / mL 2 .The solid line shows the prediction of our model: For k z,0 < k c , the particles evenly probe momentum states up to k c , so 〈E 〉 = E c /3 independently of k z,0 (blue dashed line), while for k z,0 > k c , we have 〈E 〉 = ℏ 2 k 2 z,0 /(2m) (red dotted line), so at k z,0 = k c the energy jumps by a factor of 3. (b) Extracted values of k c for all the simulation parameters used in Fig. 5(b).The solid lines are guides to the eye.
(a)].To estimate k c and its error, we use the mean and the standard deviation of 〈E 〉 for k z,0 < 10 k 0 .Fig. A1(b) shows the values of k c calculated for the values of U and ω used in Fig. 5(b).

1 .
Calculations for D = S...SD a. Calculation of T D