Теоретическая модель описала трехмерный переход Андерсона в атомных волнах

Abstract In a recent experiment (Yu et al 2026 arXiv: 2602.07654), energy-resolved measurements of an atomic matter wave spreading in a speckle potential enabled the direct observation of the three-dimensional Anderson transition. In this work, we present a quantitative theoretical description of the matter-wave dynamics based on a tailored implementation of the self-consistent theory of localization, which incorporates both the spectral and spatial properties of the state prepared in the disorder. We benchmark this theoretical approach against ab initio numerical simulations, and use it to analyze the atom density profiles observed experimentally in the localized, diffusive, and critical regimes. Particular emphasis is placed on the key role of the atomic energy distribution, especially on the distinct contributions of Bose-condensed and thermal atoms to interpret the experimental profiles. Our framework provides a versatile and efficient theoretical toolbox for quantitatively describing wave-packet dynamics in three-dimensional disordered quantum systems, which remains challenging for state-of-the-art large-scale numerical simulations. Original content from this work may be used under the terms of the Creative Commons Attribution 4.0 license. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI. 1. Introduction Understanding quantum transport in disordered media is a long-standing challenge, central to condensed-matter [1], wave physics [2], and cold-atom experiments [3]. A key phenomenon in this context is Anderson localization, whereby interference effects suppress diffusion and can ultimately lead to the absence of transport [4–6]. Ultracold atoms provide a particularly powerful platform to investigate this physics, as disorder, dimensionality, and interactions can be controlled with high precision, while observables such as density profiles are directly accessible [7–11]. These capabilities have stimulated an intense theoretical activity [12–26], and enabled experimental demonstrations of Anderson localization in one-dimensional [27, 28], two-dimensional [29] and three-dimensional (3D) [30–33] disordered potentials, as well as precise studies of analogous dynamical localization phenomena in kicked-rotor systems [34–39]. In a recent experiment [40], atoms were transferred using a radio-frequency (rf) scheme into a 3D disordered potential with a well-controlled and narrow energy distribution [41, 42]. This protocol enables the exploration of wave-packet dynamics at selected energies across the localized, diffusive, and critical regimes of the Anderson transition. In particular, it overcomes difficulties encountered in previous setups [30–33], which were limited by a broad atomic energy distribution that prevented the direct observation of the transition and the study of the critical regime (see, e.g. [23]). While the measurements reported in [40] provided a direct observation of the Anderson transition, a fully quantitative connection to microscopic theory remains to be established, in particular for an accurate description of the experimentally observed spatial density profiles. This calls for the development of a dedicated theoretical framework capable of reliably describing wave propagation over large distances and long time scales across the mobility edge in three dimensions. Such a framework must incorporate the key experimental ingredients, notably the spatial and spectral properties of the initial states prepared by the rf transfer in the disorder. In the present work, we undertake this task by developing a tailored implementation of the 3D self-consistent theory (SCT) of localization [43–45], designed to capture the above features of wave-packet dynamics. We first benchmark this approach against ab initio numerical simulations of 3D wave-packet spreading, which explicitly implement the rf-loading scheme and reproduce the statistical properties of the experimental speckle disorder. We then confront the theoretical predictions with experimental measurements of the density profiles across the mobility edge. This analysis reveals the central role of the atomic energy distribution in shaping the spatial structure of wave packets in the disordered potential. The article is organized as follows. In section 2, we recall the main elements of the experiment reported in [40], which enabled the realization of energy-resolved propagation of atomic wave packets in a disordered potential. We then develop a theoretical framework adapted to the experimental scenario, based on the SCT of localization. In section 3, we benchmark this approach against ab initio numerical simulations performed under conditions close to those of the experiment. Section 4 is devoted to the comparison with the experimental measurements, including a characterization of the spatial and spectral distributions of the atomic cloud transferred into the disorder. In particular, we discuss the respective contributions of condensed and thermal atoms to the energy distribution. Finally, we summarize our findings and outline perspectives in section 5. 2. Energy-resolved wave packet propagation across the Anderson transition: theory 2.1. Experimental context In the experiment of [40], the dynamical spreading of ultracold atomic clouds prepared with well-defined energies in a three-dimensional random speckle potential was investigated in the vicinity of the Anderson mobility edge. Such energy-resolved transport was achieved using the scheme illustrated in figure 1: atoms from a Bose–Einstein condensate (BEC) are first prepared in a disorder-free hyperfine state , then transferred to a disorder-sensitive state at a well-defined target energy . The transfer is driven by a radio-frequency pulse whose frequency sets the detuning from the bare resonance of state . As shown in [40, 41], in the weak-coupling regime where only a small fraction of the BEC is transferred, and provided that the condensate is large enough to be approximated by the uniform plane-wave state , the atoms loaded into the disorder acquire the following energy distribution: Here is the zero-momentum spectral function in the disordered potential, i.e. the probability density that an initial plane-wave state has energy in the random potential [41], and is an energy-filter function determined by the temporal shape of the rf pulse. The filter is peaked around the target energy , with a typical width . Note that the spectral function alone would correspond to the energy distribution obtained after an abrupt switch-on of the disorder applied to the BEC (see, e.g. [31]). This protocol thus enables the preparation of atomic wave packets with a much narrower energy distribution than in previous experiments. Moreover, the target energy can be finely tuned across the mobility edge , allowing for a precise investigation of the Anderson transition. 2.2. Energy-resolved state: model Having introduced the experimental protocol, we now develop a simple model describing the preparation and the subsequent evolution of an energy-resolved state in the random potential. We denote by the condensate wave function prior to the transfer, and by the state of the atoms loaded into the disorder. For a coherent transfer, as realized in [40], we model the relation between these two states by where is the Hamiltonian in the disorder-sensitive state, with the speckle potential, and is an amplitude filter function such that . In the second equality, we have introduced a complete set of eigenstates of with energies , which will prove convenient in the following. Experimentally, the filter function is directly linked to the Fourier transform of the temporal envelope of the rf pulse, as given by first order perturbation theory for an initial state transferred to a continuum. Equation (2) is therefore valid in the weak-coupling regime and provided interactions can be neglected. One can easily verify that the definition (2) correctly yields the energy distribution (1) of the atoms in the disorder. Indeed, we have by definition where equation (2) has been used in the second equality, together with the approximation for the condensate wave function in state . Here and in the following, the overline denotes averaging over disorder realizations. 2.3. Disorder-averaged density We now derive an expression for the disorder-averaged density , suitable for the application of the SCT of localization, starting from the filtered initial state (2). The wave function of the cloud in the random potential at time is related to the initial state by where is the Heaviside step function (unless otherwise stated, from now on we set ). By inserting the closure relation and using the spectral representation the wave function can be rewritten as where we have introduced the retarded Green’s function From equation (6), we obtain a general expression for the disorder-averaged density: where is the advanced Green’s function, and we have changed the energy variables to and , with . Although equation (8) is fully general, it is not convenient for practical use, as it involves a correlator of Green’s functions evaluated at three distinct spatial points. However, in the so-called hydrodynamic regime, where fast spatial and temporal variations of the profile (typically occurring on the scale of the mean free path and the mean free time, respectively) are neglected, this expression can be significantly simplified. The derivation is presented in appendix In this relation, is the energy distribution (1), and is the disorder-averaged propagator, defined as with the density of states per unit volume. Physically, gives the probability density for a particle with energy to propagate from at time to at time in the random potential [47]. The averaged density (9) then follows by summing over all possible initial positions and all energies, weighted by the energy distribution of the atoms. In the ideal limit of a very sharp filter function , equation (9) further simplifies to so that the density becomes a direct probe of transport in the disorder at energy . In section 3, equation (11) will be compared with numerical simulations performed in the ideal sharp-filter limit. By contrast, when comparing with the experimental data in section 4, it will be necessary to account for the finite width of the energy filter. We will therefore rely instead on equation (9). 2.4. Propagator: SCT We now derive a general expression for the propagator capturing the wave-packet dynamics across the mobility edge of the Anderson transition. To this end, we use the SCT of localization, a simple model in which the quantum interference responsible for Anderson localization are modeled by a re-summation of Cooper diagrams [43–46]. Within this approach, the Fourier transform of the propagator (10) is given by involving a generalized diffusion coefficient renormalized by quantum interference according to: where is the so-called Boltzmann or classical diffusion coefficient and the density of states. An exact calculation of the frequency-dependent diffusion coefficient and the corresponding density was previously performed in [48] in one and two dimensions. Here, we extend this analysis to the 3D case. To this end, we first note that the integral on the right-hand side of equation (13) is divergent and must be regularized by introducing a UV cutoff , typically of the order of the inverse mean free path. This can be conveniently achieved by isolating the divergence in the momentum integral as follows [49]: Carrying out the two integrals yields a third-order algebraic equation for , whose solution is where and (restoring ). In this formulation, the three physical quantities , and of the model are recast in terms of the three composite parameters , and . The latter, in particular, represents the mobility edge. While these parameters can be approximately estimated using perturbation theory [47, 50, 51], they are difficult to determine accurately in the strong-disorder regime where the Anderson transition occurs. Here we will therefore simply take them as constant free parameters of the theory, exploiting that they vary smoothly with energy in the vicinity of the mobility edge we are focusing on. Together, equations (1), (9), (12) and (15) provide a complete description of the wave-packet density at any position and time across the critical point. It is nevertheless instructive, at this stage, to recall the predictions of the theory in the limit of asymptotically long times, corresponding to [52]: In terms of wave-packet spreading, these frequency scaling laws imply fundamentally different dynamics around the mobility edge, which can be traced in the time evolution of the mean squared width of the wave packet, . Assuming a sharp filter function, so that the wave packet has a well defined energy , readily follows from equations (11) and (12): Inserting the asymptotic laws (16), we infer: which are the usual diffusive, anomalous and localized behavior expected above, at and below the mobility edge [5, 53]. 3. Benchmarking the theory against numerical simulations Before analyzing the experimental results of [40] within the theoretical framework introduced in the previous section, we first assess its validity through comparison with ab initio numerical simulations. These simulations are performed in an idealized setting assuming sharp energy filtering and a spatially narrow initial state, yet under conditions close to those of the experiment. To this end, we numerically prepare, in a blue-detuned speckle potential , an initial filtered state of the form given in equation (2), and then compute its time evolution in the random potential. The density is monitored as a function of time and subsequently averaged over disorder realizations, yielding the average density . 3.1. Speckle potential The speckle potential is generated numerically on a regular rectangular grid of volume , using a Fourier filtering method similar to the one described in [41]. The disorder corresponds to a blue-detuned laser speckle, as in [40], implying that the potential is strictly positive. Its statistical distribution follows the exponential law [7]. Here, denotes the root-mean-square (rms) value of the disorder and is referred to as the disorder amplitude. For laser speckle potentials, also coincides with the mean value of the disorder, i.e. . In the following we present numerical results for , as in [40]. The speckle potential is generated to accurately reproduce the experimental configuration, where the speckle grains are elongated along the axis (see figure 2). More precisely, the parameters are chosen such that the autocorrelation function (where ) matches the experimental one both along the longitudinal axis and in the transverse – plane. As shown in figure 2(c), the agreement between the numerical computation, the theoretical expression given by equation (42) in the appendix C.1, and the experimental measurements [40, 41] of the autocorrelation, is excellent. Consistently with [23], the correlation lengths are defined as the half width at half maximum divided by a factor 1.39156. It yields and on the transverse plane. These correlation lengths set the correlation energy , where is the geometric mean over the different directions. Download figure: Standard image High-resolution image In the simulations, the discretization steps are typically and , yielding , . These ratios are sufficient to properly resolve both the speckle grains (see figure 2(c)) and the de Broglie wavelength at the considered energies. Details of the numerical generation of the speckle potential and the theoretical form of the autocorrelation function are provided in appendix C.1. 3.2. Numerical propagation of an energy-resolved wave-packet To numerically prepare the initial state defined by equation (2), we start from a Gaussian wave packet , and consider a narrow Gaussian filter , with and Hz. The target energy is varied around the expected mobility edge, typically from to Hz. To apply the filtering operation involved in equation (2), we do not compute the eigenstates of the random Hamiltonian, which would be numerically prohibitive for the 3D system considered here. Instead, we use an expansion of the operator on a basis of Chebyshev polynomials. This strategy was previously employed in [54–56] in the context of Anderson localization in momentum space, and is briefly described in appendices C.2 and C.3. Once the energy-resolved state is computed following equation (2), we numerically propagate it in time by applying the time evolution operator . The time propagation is performed with high accuracy and efficiency using a Chebyshev expansion of the evolution operator, which is applied iteratively to obtain the state at successive time steps . This method, briefly reviewed in appendices C.2 and C.5, is more stable and efficient than Crank–Nicolson schemes, as it allows for relatively large time steps [57, 58]. Within this approach, we obtain the spatial density at time , which we average over realizations of the disorder to get the average density . To characterize the dynamics, we follow the methodology of [40] and we focus on the disorder-averaged linear density profile along , defined as . Figures 3(a)–(c) show such linear density profiles at different evolution times obtained from the numerical simulations (colored points) for three target energies , corresponding respectively to the diffusive (), critical () and localized regime (). For the chosen numerical parameters, we estimate the mobility edge to be . This value slightly exceeds the numerical estimate reported in [23], which is based on extensive transfer-matrix simulations. We attribute this discrepancy to the finite size of the numerical grid, which restricts the accessible expansion times, and to the relatively coarse spatial discretization required to keep the computational cost of long-time expansion in a 3D disordered potential tractable. Despite this, the simulation results clearly reveal the characteristic Gaussian spreading and the time-independent exponential profile expected on the diffusive (figure 3(a)) and localized (figure 3(c)) sides of the Anderson transition, respectively. Near the mobility edge (figure 3(b)), the profile exhibits an intermediate shape together with a slow expansion. Download figure: Standard image High-resolution image 3.3. Comparison with the 3D SCT To compare the numerical results of figure 3 with the SCT, we start from the narrow-filter approximation (11) applied to the linear density, where . The 1D propagator is deduced from equation (12): with the diffusion coefficient determined from equation (15). Note that the propagator is normalized, , and so is the linear density . Because the profiles represented in figure 3 typically correspond to times well beyond to the mean free time 4, equation (20) can be further simplified by using the small-frequency asymptotic relations (16) for . In the diffusive regime , they lead to the well-known Gaussian distribution where , while in the localized regime one has where . At the mobility edge , finally, the linear density is an Airy function [59, 60]: with . Figures 3(a)–(c) show as solid curves the SCT prediction (19) for the linear density, using the asymptotic expressions (21)–(23). For the comparison, we use the average numerical profiles immediately after loading into the disorder as the initial density distribution, and take and as fit parameters. These are calibrated from a global comparison between the SCT and the numerical profiles. To ensure consistency when comparing the SCT to both the numerical and experimental profiles (see section 4.3), we use the central density as the fitting observable. This is motivated experimentally since the central density is less noisy than, for example, the tails of the profiles (which are affected by finite signal-to-noise ratio), while still fully capturing the dynamics of the cloud expansion. The parameters and are determined by minimizing the least-squares cost function where denotes here the numerical central density at energy , its associated uncertainty and the corresponding SCT prediction for a given parameter pair . The double sum runs over target energies spanning the mobility edge and over evolution times up to five seconds. The minimization yields the values , where the quoted uncertainties correspond to the fitting errors. Overall, figures 3(a)–(c) show very good agreement between the SCT and the ab initio simulations over the range of times and energies explored (typically , for evolution times up to s). In particular, the tails of the density profiles are very well reproduced by the theory, both in the localized regime and at the critical point—an encouraging result given the approximate nature of the SCT. Note that despite the relatively narrow initial density distribution, the profiles remain close to the initial shape in the localized regime for the chosen disorder amplitude. This behavior is also observed experimentally (see section 4). In the diffusive regime, by contrast, we observe significant deviations between theory and simulations in the tails. These discrepancies originate from the finite size of the numerical domain, which starts to affect the profiles once the wave packet reaches the domain boundaries. This interpretation is supported by theory: replacing the unbounded expression (21) with the diffusive solution in a finite box of size (with periodic boundary conditions, as in the simulations), we find [47] which yields an excellent agreement with the simulations, including in the tails (see figure 3(a)). Figure 3(d) finally shows the inverse central density of the wave packet, , as a function of time for three target energies in the diffusive, localized, and critical regimes. The plot reveals the three characteristic long-time scaling laws (18) for the packet width, , with again very good agreement between the SCT and the ab initio simulations. This analysis confirms that the theoretical framework reliably captures the dynamics across the mobility edge, even in a regime of relatively strong disorder, large spatial scales, and long times, where direct numerical simulations become highly demanding. 4. Theory versus experiment We now present comparisons between the SCT and the experimental observation of the density-profile dynamics obtained with the setup of [40]. We begin with a brief description of the experimental conditions (section 4.1) and of the corresponding initial density profiles, which serve as the input to the SCT. We then characterize the atomic energy distribution realized in the experiment, accounting for the residual thermal background cloud (section 4.2). Finally, we provide a quantitative comparison between the measured density profiles and the theoretical predictions over a range of target energies spanning the mobility edge. 4.1. Energy-resolved state preparation The crucial feature of the experiments reported in [40] is the energy-resolved state preparation, as sketched in figure 1. The sequence starts with the preparation of a BEC in state , containing around atoms in a nearly isotropic optical trap with frequencies Hz. The temperature is and the condensate fraction is around . The BEC is characterized by the chemical potential , yielding the Thomas–Fermi radii . The state-dependent disorder is created using a bichromatic optical speckle, following the principle described in [42]. A small fraction of atoms () is then transferred to the state by applying a radio-frequency (rf) pulse. The rf frequency sets the target energy to , where is the detuning from the bare resonance in absence of disorder (see figure 1). The rf pulse duration, ms, defines the energy width of the energy-resolved wave packet (see below). Immediately after the transfer, at , the optical trap is switched off and the remaining atoms in state are removed. The subsequent evolution of the wave packet in the disordered potential is then monitored up to 5 s by in-situ fluorescence imaging along the axis. During the whole experimental sequence, the atoms are suspended against gravity using a magnetic levitation scheme. Such levitation yields a shallow residual harmonic confinement5 in the horizontal - plane, such that the expansion occurs essentially along the vertical axis. Following the analysis done in [40], we then integrate further the recorded column density , and we study the one-dimensional profile along . For each energy , we record the initial profile immediately after transfer. Examples are shown in figure 4 for , and Hz (solid curves with orange shading). The dashed lines indicate the estimated initial density profile of the BEC in state before transfer, including the thermal component. The experimental profiles exhibit weak tails, which we attribute to a combination of short-time interaction effects and possible dynamics in state during the rf-coupling stage. They also display a slight asymmetry, likely caused by small fluctuations and inhomogeneities in the levitation magnetic field. Download figure: Standard image High-resolution image 4.2. Energy distribution The key objective of the rf transfer scheme is to prepare an atomic cloud with a narrow energy distribution in state . In the weak-coupling regime in which the experiment operates [40], the energy spread of the transferred atoms is Fourier-limited, . In this case, and assuming that the atoms in state form a pure BEC, the transferred state is also pure and has an energy distribution given by equation (1). The filter function is determined by the Fourier transform of the temporal envelope of the rf pulse, which has the shape of a Kaiser function [40]. This leads to a narrow rms energy width . In the experiment, however, the initial atomic cloud in state is not a fully pure BEC but it contains a finite thermal fraction, . This thermal component gives rise to an additional contribution in the energy distribution, which generalizes to: where (see appendix The thermal contribution describes the rf transfer into the disorder of the thermal atoms initially present in state . As shown in appendix Figure 5 shows energy distributions calculated numerically from equation (26) for three target energies, , and (indicated by vertical dashed lines), for a disorder amplitude of and a thermal fraction of . The shaded area represents the zero-momentum spectral function that is used to calculate the energy distribution from equation (1) (see appendix C.4 for details on the numerical computation of the spectral function). Although the thermal contribution yields a broad pedestal (see inset of figure 5), it remains small enough such that the energy distribution remains essentially very narrow, one [31, 32] or even two [30, 33] orders of magnitude smaller than in previous experiments. Download figure: Standard image High-resolution image 4.3. Comparison with 3D SCT We now compare the experimental density profiles after evolution in the speckle disorder with the SCT predictions. For this comparison, the SCT is supplied with three quantities: the initial experimental density profiles (figure 4), the energy distribution computed numerically from equations (26) and (27), and the experimental value of the mobility edge, . The latter was measured independently in [40] through a study of the inverse central-density dynamics, , thereby providing an unambiguous identification of the localized and diffusive phases. Within our theoretical framework, the 1D integrated density profile after an evolution time follows from equation (9) as where is the disorder-averaged propagator defined in equation (20). To compute this integral, we use for the initial distribution fits of the experimental density profiles measured immediately after the rf transfer to an asymmetric pseudo-Voigt function (see [61]). Examples of these fits are shown as solid lines in figure 4. As explained in section 2.4, the SCT depends on three effective parameters: the mobility edge , and the coefficients and . Here the mobility edge is fixed independently using the experimental determination reported in [40] (and in agreement with transfer-matrix computations [23]). The parameters and are calibrated according to the same method explained in section 3.3, where the sum in equation (24) runs over target energies spanning the mobility edge and over evolution times up to three seconds. This procedure is applied both to the ideal energy distribution —corresponding to a pure BEC in state —and to the full distribution that accounts for the 25% thermal fraction (see equation (26)). We first show in figure 6 the experimental density profiles of atoms in state along the -axis at the maximum evolution time , for the three energies , and Hz considered so far. The pronounced change in the profile shapes illustrates the expected behavior across the mobility edge : from a sharply peaked profile with exponentially decaying tails in the localized regime, to a broad profile in the diffusive regime. The experimental data are compared with the theoretical prediction (28) using either (, dashed curves) or (with , solid curves) for the energy distribution. In the diffusive regime, both choices provide a good description of the data. However, at the mobility edge and in the localized regime, significantly better agreement—particularly in the tails of the density profiles—is obtained when the thermal fraction is properly taken into account. This confirms the crucial role of the atomic energy distribution and highlights the need for an accurate characterization of the initial cloud in order to reproduce the observed spatial profiles (see also [31–33, 62]). Download figure: Standard image High-resolution image We note that even when thermal atoms are included, a residual discrepancy remains in the density tails in the localized regime. This may originate from effects present in the experiment but not included in the theory, such as interactions during the rf transfer or the residual shallow harmonic confinement in the – plane. These mechanisms may also account for the values and , which are slightly smaller, though of the same order of magnitude, than those obtained under ideal conditions in the numerical simulations of section 3. We further show in figure 7 the time evolution of 1D density profiles for the same three energies (points, experimental data), together with the prediction of the SCT taking into account the finite thermal fraction of the cloud (solid curves). The plots clearly reveal atomic spreading in the diffusive regime, while the profiles remain nearly time-independent in the localized regime. Figure 7(a) also suggests that, in the localized regime, the core of the distribution appears slightly narrower than the initial density profile shown in figure 4. We tentatively attribute this feature to the intrinsically strong spatial fluctuations of the localized eigenstates onto which the density profile decomposes. In the experiment, the averaging over disorder realizations is only partial, so that residual signatures of these fluctuations may persist in the measured profiles. Such effect lies beyond the scope of the SCT, which essentially provides a mean-field description of the Anderson transition, but would be worthwhile to investigate in future work. Download figure: Standard image High-resolution image 5. Conclusion In this work, we have developed a flexible theoretical framework for quantitatively describing the dynamics of matter-wave packets in three-dimensional random potentials. The approach combines an accurate characterization of the atomic energy distribution with a three-dimensional implementation of the SCT of localization, providing access to the dynamics across the mobility edge. Although the SCT is intrinsically approximate, we have shown that it accurately reproduces exact 3D numerical simulations in a nontrivial anisotropic speckle potential, while being significantly less computationally demanding. Within this framework, we have further achieved quantitative agreement with the density profiles reported in the recent experiment of [40], where the dynamics near the Anderson transition could be probed using well-controlled atomic energy distributions. Beyond the present study, our methodology can be readily extended to other platforms, such as fermionic systems [30, 33]. It would also be interesting to generalize the approach to more complex scenarios, for instance by including a shallow transverse confinement, as in the experiment of [40], or by incorporating weak atom-atom interactions within the SCT, following the lines of previous works [46, 63]. On the experimental side, it would be highly appealing to apply the rf-transfer scheme to other physical problems. Relevant examples include measurements of critical exponents or multifractality [60, 64, 65], investigations of percolation transitions in random potentials [66] and their connection to Anderson physics [26], exploration of localization phenomena in momentum space [17, 67, 68], and characterization of mobility edges and critical phases in one-dimensional quasi-periodic potentials [69–72]. Acknowledgment The authors thank V Denechaud, B Lecoutre, A Signoles for early work and D Clément, D Delande and M Filoche for fruitful discussions. Data availability statement The data cannot be made publicly available upon publication because the cost of preparing, depositing and hosting the data would be prohibitive within the terms of this research project. The data that support the findings of this study are available upon reasonable request from the authors. Funding This work benefited financial support from the project Localization of Waves of the Simons Foundation (Grant No. 601 939), the French National Research Agency under Grant No. ANR-24-CE30-6695 (FUSIoN), Grant No. ANR-25-CPJ1-0059-01 (CPJ LUMIS), Grant No. ANR-22-CMAS-0001(QuanTEdu-France) and Region Île-de-France in the framework of DIM QuanTiP. Author contribution J-P B, S B and N C conducted the numerical simulations and developed the theoretical predictions. K X, H-M Q, X Y, M N and Y G carried out the experiments. All authors contributed to data analysis, project development, and the writing of the paper. Appendix A: Disorder-average density In this appendix, we provide details on the derivation of the disorder-average density, equation (9). Our starting point is the general relation (8), in which we express the Green’s function correlator in the hydrodynamic regime as [47] Here () is the averaged retarded (advanced) Green’s function, is the quantum propagator, which gives the probability density for an atom of energy to propagate from to , is the scattering mean free path and the density of states per unit volume. Next, we insert this decomposition into equation (8) and introduce the Wigner change of variables . This leads to Here, we have further approximated , using that in the long-time dynamics, , the dominant contribution to the integral over arises from the vicinity of owing to the rapidly oscillating factor , so that . To simplify this expression, we proceed in two steps. First, instead of directly working with the product of initial wave functions in real space, it is convenient to introduce the Wigner distribution of the initial state, defined by the relation: The main interest of the Wigner distribution is to simultaneously encode the spatial and momentum structures of the initial state. In particular, it gives access to the spatial density after integration over the momentum variable: Second, we use that at scales larger than the mean free path (hydrodynamic regime), is a smooth function of and , so that . With this simplification, the three integrals over , and in equation (30) can be easily performed, yielding: In the last step, we note that the momentum width of is of order , where is the scattering mean free path, whereas the Wigner distribution has a momentum width , with the spatial extent of the initial wave packet. In the experiment [40], one typically has , so that in equation (33) one can approximate . This enables us to perform the integral over using equation (32). Noting further that , where is the spectral function, we recover equation (9) of the main text. Appendix B: Energy distribution of the atomic cloud In this appendix, we present a derivation of the energy distribution (26), which accounts for a finite thermal contribution given by equation (27). The starting point is the density matrix of the Bose gas in state prior to transfer: where is the free-space Hamiltonian, a normalization factor, and we have approximated for the BEC wave function. The rf pulse transfers an initial energy in state to a final energy in the disorder-sensitive state . For an initial state, , and the transfer in the disorder corresponds to the transformation , which generalizes equation (2), with the transferred energy. The energy distribution in the disorder after transfer is defined as The initial BEC state has energy , so its energy distribution after transfer in the disorder reads which immediately leads to equation (1) of the main text. Consider now the energy distribution of the thermal component. Since the thermal factor in equation (34) involves the free Hamiltonian, it is convenient to introduce the basis states to evaluate the trace: We then approximate , and change variable in the momentum integral from to . This leads to: with the free-space density of states in three dimensions. We finally perform the integral over using that the filter function is much narrower than the spectral function, the thermal factor and the density of states: Using that in three dimensions (with the Heaviside step function), and imposing the normalizing condition , we finally obtain equations (27) of the main text. Appendix C: Numerical methods C.1. Numerical generation of the speckle potential Experimentally, the speckle potential is generated by illuminating a diffusive lens with a laser beam of vacuum wavelength nm (wavenumber ). The incident beam is Gaussian, with a waist of mm at the lens plane, and is truncated by a diaphragm of diameter mm. The lens has a focal length mm. The beam propagates along the -axis and is linearly polarized in the -plane. Numerically, the speckle potential is computed using a Fourier-filtering method similar to that in [41]. The electric field along the polarization direction is expressed in the Fourier representation as Here, is an arbitrary amplitude, and is a uncorrelated, complex Gaussian stationary random process with zero mean. The wave vector denotes the projection of onto the -plane, and is the cutoff set by the diaphragm. The Dirac distribution enforces the light’s dispersion relation in vacuum. Numerically, the integral (40) is evaluated using a 3D fast Fourier transform on a rectangular grid, approximating the Dirac distribution by a narrow Gaussian of width : . Equation (40) implies that the field is a sum of independent Gaussian random variables . Consequently, the probability density of the potential , where is a constant characterizing the light-atom interaction, follows the exponential distribution [73]. The average potential is obtained from , using the uncorrelated nature of : . This yields , where To numerically generate the speckle potential shown in figures 2(a) and (b), we proceed as follows. We evaluate equation (40) with using a numerically generated random noise (the corresponding numerical parameters are given in section 3.1). The resulting field is normalized by , and the potential is then obtained as , where is the experimental rms disorder strength. From the numerically generated speckle potential, the autocorrelation function shown in figure 2(c) is computed as . In the same figure, this result is compared with the theoretical prediction, which follows directly from equation (40) using the Gaussian character of : C.2. Chebyshev polynomial expansion of operators The numerical computation of the initial energy-resolved state, the spectral function, and the time evolution of the wave packet presented in the main text rely on a unified framework based on Chebyshev polynomial expansions of operator functions of the Hamiltonian. The calculations are performed on the spatial grid used to generate the speckle potential, with the Laplacian term of the Hamiltonian discretized using a three-point finite-difference stencil in each spatial direction. In the three cases, the problem reduces to applying a given operator function of the Hamiltonian, , to an initial state. Since full diagonalization of the 3D Hamiltonian would be computationally prohibitive, we approximate by a truncated Chebyshev expansion. To this end, the discrete Hamiltonian is first shifted and rescaled so that its spectrum lies in the interval , where Chebyshev polynomials are defined. Given bounds and of the spectrum, we define and introduce the rescaled operator , whose eigenvalues lie in . Using together with the Chebyshev expansion for the Dirac distribution [74], one obtains the Chebyshev expansion where denotes the Chebyshev polynomial of the first kind of order . Applying to a state thus amounts to computing the vectors , generated recursively via Each new vector requires only a sparse matrix-vector multiplication and a vector subtraction, making the method both efficient and numerically stable. In practice, the series is truncated at finite order to reach the desired accuracy and evaluated using the Clenshaw algorithm [75]. C.3. Numerical generation of the energy-resolved state We first show how the Chebyshev method is used to numerically generate the energy-resolved state , defined by equation (2). In this problem, on applies the operator to the BEC state prior to transfer, where is the filter function discussed in the main text. In this case, we use , where the coefficients are computed from equation (45). The integral defining in equation (45) is conveniently evaluated by introducing the change of variable . Using the identity , the coefficients can be written as the cosine transform Numerically, the variable is sampled uniformly in the interval with points, (), and the integral is approximated by the discrete cosine transform where and . C.4. Numerical generation of the spectral function The computation of the energy distribution appearing in equation (26) requires the numerical evaluation of the spectral function . This is performed using a Chebyshev expansion with . To this end, we write: where is the spectral function for the rescaled Hamiltonian, and we have introduced . Using the Chebyshev expansion of the Dirac distribution [equation (44)], we obtain To suppress Gibbs oscillations arising from truncation at finite order , the coefficients are multiplied by Jackson weights. This procedure—equivalent to convolving the spectral function with an optimal smoothing kernel—is known as the Kernel Polynomial Method [74]. For the energy distributions shown in figure 5, we use , ensuring that the Jackson broadening remains well below the intrinsic width of the filtered spectral functions (below for the narrowest spectral function of 5 Hz). C.5. Time evolution The initial state is propagated in time with step according to . Numerically, the evolution operator is applied efficiently using a Chebyshev expansion with [57]. To this end, the evolution operator is first rewritten as where . The first factor in equation (52) is a global phase, which leads to an irrelevant energy shift and has no impact on the dynamics. The second factor, denoted , can be expanded as where equation (44) has been used in the second equality. The coefficients are given in closed form as where denotes the Bessel function of the first kind of order . Applying the evolution operator thus amounts to computing where the vectors are generated recursively using equation (46). Owing to the rapid decay of the Bessel functions with , the series can be safely truncated without additional regularization. Footnotes - 4 - 5Although the precise influence of the weak confinement remains to be precisely investigated, the transport dynamics remains fully 3D since the transverse size of the cloud is much larger than the expected scattering mean path , see [40].