{'fraction_non_alphanumeric': 0.0360340541610753, 'fraction_numerical': 0.014233231405504103, 'mean_word_length': 4.660316274436559, 'pattern_counts': {'":': 0, '<': 7, '<?xml version=': 0, '>': 14, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 1, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 6, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Long-time expansion of a Bose-Einstein condensate: Can Anderson localization be observed?

We numerically explore the long-time expansion of a one-dimensional Bose-Einstein condensate in a disorder potential employing the Gross-Pitaevskii equation. The goal is to search for unique signatures of Anderson localization in the presence of particle-particle interactions. Using typical experimental parameters we show that the time scale for which the non-equilibrium dynamics of the interacting system begins to diverge from the non-interacting system exceeds the observation times up to now accessible in the experiment. We find evidence that the long-time evolution of the wavepacket is characterized by (sub)diffusive spreading and a growing effective localization length suggesting that interactions destroy Anderson localization.

We numerically explore the long-time expansion of a one-dimensional Bose-Einstein condensate in a disorder potential employing the Gross-Pitaevskii equation. The goal is to search for unique signatures of Anderson localization in the presence of particle-particle interactions. Using typical experimental parameters we show that the time scale for which the non-equilibrium dynamics of the interacting system begins to diverge from the non-interacting system exceeds the observation times up to now accessible in the experiment. We find evidence that the long-time evolution of the wavepacket is characterized by (sub)diffusive spreading and a growing effective localization length suggesting that interactions destroy Anderson localization.

I. INTRODUCTION
Anderson localization in disordered systems has been predicted almost 60 years ago. Anderson could show that quantum transport is suppressed in the presence of disorder. The particle remains trapped with exponentially decaying probability density [1]. This phenomenon turned out to be quite universal as it occurs in a wide variety of systems whenever a multitude of scattered paths can destructively interfere with each other. Localization effects have been observed in systems as diverse as electrons in solids (for a review see [2]), light in disordered media [3][4][5], microwaves [6], sound waves [7] and in driven timedependent systems where it is referred to as dynamical localization [8,9]. Many open questions have remained: one of the most important is the role of particle-particle interaction. In its original formulation Anderson localization has been identified within the one-particle Schrödinger equation. Going beyond the one-particle picture, does the many-body wavefunction of the interacting system support Anderson localization, i.e. the exponential decay of the reduced one-body density? The study of the influence of interactions on Anderson localization has a long history (e.g. [10,11], for a review of early results see [12]), but has mostly been confined to time-independent condensed matter electronic systems near their ground state. The underlying manybody localization can be viewed in terms of localization in Fock space close to the Hartree-Fock state of localized single particle states [13]. For one-dimensional (1D) disordered systems a transition to delocalization has been found for increasing attractive interactions for fermions and increasing repulsive interactions for bosons [14,15]. The ground state of the disordered Fermi-Hubbard model has been found to show different phases, a metallic phase, a Mott insulator phase, and an Anderson localized phase, * iva.brezinova@tuwien.ac.at depending on the interplay between the strength of disorder and on-site interactions [16]. Similar quantum phase diagrams have been investigated for the ground state of ultracold Bose gases in 1D [17]. Several theoretical (see e.g. [18][19][20]) and experimental studies (see e.g. [21]) have investigated this transition from localized to delocalized many-body states in the framework of spin chains and ultracold atoms. For finite temperature (T > 0) ensembles, it has been shown that at sufficiently high temperatures interactions lead to dephasing and finite conductivity [22]. It was further found that electron-electron interactions yield a finite hopping above a critical temperature [13,23]. A similar result, i.e. a metal-insulator transition, has been found recently for weakly interacting bosons in one dimension [24]. Recently, Bose-Einstein condensates (BECs) have offered a new platform to investigate not only Anderson localization in its pure form, i.e. without interactions, but most importantly, in the presence of interactions [25][26][27][28][29][30][31]. Atoms in a typical BEC interact via van der Waals interactions which can be approximated by contact interactions in the low temperature limit with the s-wave scattering length entering as the only parameter. Unless interactions are switched off by tuning Feshbach resonances in an external magnetic field, the interactions between ultracold atoms remain non-negligible [28,29]. While van der Waals or dipolar interactions are fundamentally different from the (screened) Coulomb electronelectron interaction in solids, investigation of interacting ultracold atoms in a BEC allows to probe fundamental questions about the role of interactions on Anderson localization with unprecedented experimental control with far-reaching implications also for condensedmatter physics. Numerical transport studies for 1D BECs through correlated disorder potentials showed a crossover between an Anderson localization like phase for weak interactions and a delocalized regime for larger interactions [32][33][34][35]. The localization dynamics of timedependent out-of-equilibrium systems has remained an open problem. Numerical simulations have been mostly performed within the framework of the discrete nonlinear Schrödinger equation where subdiffusive spreading of the mean width of an initially spatially localized wavepacket has been predicted and numerically observed [36][37][38][39][40]. Closely related studies focus on the subject of diffusion or suppression of transport in out of equilibrium disordered classical chains [41][42][43]. In the present paper, we aim for a numerical simulation following the experiment [27] in which first indicators of Anderson localization in an expanding BEC in one dimension were reported. The goal is to explore the long-time evolution of such a non-equilibrium system and to identify possible signatures of Anderson localization. This requires the solution of the continuous rather than the discrete non-linear Schrödinger equation whose longtime propagation poses a considerable numerical challenge. In the experiment [27], the BEC is initially confined in a cigar-shaped harmonic trap. The longitudinal harmonic potential is switched off after which expansion in a speckle-like disorder potential along the longitudinal axis sets in. Observation of an approximately exponentially localized particle density of almost macroscopic extent and over a period of more than one second after switching off the longitudinal trapping potential provided first indication of Anderson localization of an expanded BEC. We solve the Gross-Pitaevskii equation (GPE) on a large spatial grid and for long propagation times to explore the effect of particle-particle interactions on the expansion. In particular, we delineate the intrinsic difficulties in identifying unique signatures of Anderson localization within such a non-equilibrium scenario resulting from the simultaneous presence of multiple time and length scales in this problem. We relate our findings to the destruction of Anderson localization by subdiffusive spreading where the mean width increases as ∆x ∝ t α as a function of time t with α < 1/2 as observed in discrete systems [38,39]. We show that both a clear distinction between the expansion of an interacting and a non-interacting BEC as well as the approach to an asymptotically localized state, if occurring at all, manifests itself only on time scales exceeding present experiments. The paper is structured as follows: In Sec. II we introduce the system under consideration, the observables, and the numerical method used to propagate the system. In Sec. III we introduce the alternative expansion scenarios in order to disentangle the roles of particle-particle interactions and disorder. Results for observables signifying the degree of localization will be presented in Sec. IV followed by conclusions in Sec. V. Further technical details are given in the appendices A to C.

II. SYSTEM, AND NUMERICAL IMPLEMENTATION
We simulate a typical expansion scenario for a BEC following the experimental scenario of Ref [27]. Accordingly, we consider a BEC of 87 Rb atoms trapped in a cigar-shaped trap with transverse frequency ω ⊥ = 70 × 2π Hz and longitudinal frequency ω 0 = 5.4 × 2π Hz. Inter-atomic interactions are governed by the s-wave scattering length which is a s = 5.8nm [44]. The characteristic scales for length, time, and energy of the initial longitudinal harmonic trap potential serve as units throughout the paper, i.e. x 0 = /mω 0 ≈ 4.6µm, t 0 = 2π/ω 0 ≈ 29ms, and e 0 = ω 0 ( = m = ω 0 = 1). The BEC wavefunction along the longitudinal axis x is then given the solution of the 1D stationary GPE where N is the particle number and µ 0 is the chemical potential. For simplicity we abbreviate the strength of the inter-particle interaction by the nonlinearity-parameter g = 2 ω ⊥ a s N . For the numerical value of the nonlinearity we use g 0 = 400e 0 x 0 which corresponds to a particle number of N ≈ 1.2 × 10 4 . Calculating the ground state with Eq. 1 one obtains for the chemical potential µ 0 ≈ 35.6e 0 which is slightly larger than the energy of the first transverse excited state in the trap e 1 ⊥ = 2 ω ⊥ + ω 0 /2 = 26.4e 0 indicating that the radial degree of freedom might not be completely negligible in the experiment, at least initially. As discussed below, the BEC released from the longitudinal trap, while keeping the radial trap intact, rapidly expands such that the chemical potential falls below e 1 ⊥ within t 0 . We therefore assume (as has been suggested in [27]) that the ensuing dynamics can be regarded as effectively one-dimensional. The solution to Eq. 1 provides the initial condition at t = 0 for the subsequent expansion dynamics. For t > 0 we describe the dynamics within the 1D time-dependent GPE as Under the condition r s a s , where r s is the mean distance between the particles, the many-body ground state in the trap corresponds to a BEC well described by the GPE [45]. In Eq. 2, the harmonic trapping potential in Eq. 1 is switched off at t = 0 and the disorder (D) potential is switched on at the time t D > 0. We also allow for a time-dependent interaction strength whereby the particle-particle interactions, represented in the GPE by the non-linearity, are switched off at t I > 0. Choosing suitable t D and t I allows us to explore the competition between disorder and particle-particle interactions during the expansion process.
For the disorder potential we consider both Gaussian correlated disorder potentials and speckle potentials [46,47]. The Gaussian correlated disorder potential is generated by equidistantly placing Gaussian functions of width σ D and random amplitude A i The , the disorder potential becomes To be in the regime of Anderson localization [48] we The disorder is characterized by a Gaussian correlation function in real space (Fig. 1) and a Gaussian correlation functioñ also in Fourier space. For the speckle potential we follow the procedure described in Ref. [46,47]. The speckle potential is characterized by the real-space correlation function whose Fourier transform has a sharp high-momentum cut-off. We have chosen the parameters for the speckle potential such that the correlation functions in Fourier space have the same value at the maximal momentum k max = 1/ξ [48] within the Thomas-Fermi approximatioñ The numerical values used for the speckle potential are V D = 0.3[E int (t = 0) + E kin (t = 0)] and σ D = 0.57ξ. As will be discussed below, the choice of a (relatively) small value of σ D facilitates the observation of differences between the interacting and non-interacting scenarios already at (relatively) short propagation times. The choice of parameters [Eq. (11)] leads to somewhat different realspace correlation functions for the two disorder potentials (see Fig. 1). We have, however, verified that choosing different parameters for the speckle disorder such that the real-space correlation functions would more closely resemble each other does not significantly alter our results. The observables we will be focusing on are the disorderaveraged real-space particle density and momentum space density The disorder-averaged particle densities in either coordinate or momentum space displayed in the following are smoothed with Gaussians of width 6x 0 and 0.1/x 0 , respectively. For the analysis of (sub)diffusive spreading we investigate the mean width in real space with and the corresponding width in Fourier space As will be discussed below, due to the multi-scale nature of the expansion scenario, the mean width is insufficient to completely characterize the localization process. We therefore define, in addition, the inverse slope of the decaying density distribution L loc as a measure for the degree of localization (see below). For the numerical propagation we discretize space in an equidistant grid and use a spectral (Fourier) basis to represent spatial derivatives. A sufficiently dense grid was chosen such that the results have converged with respect to the grid spacing. This implies that we are in true continuum limit and are not treating a discrete non-linear Schrödinger equation as has been investigated in the past [36][37][38][39]. For all results presented here we use a box size x ∈ −2.4 × 10 4 , 2.4 × 10 4 x 0 with 1638400 grid points and periodic boundary conditions. For the time propagation we use a split-step method with adaptive step size. Details of the numerical method allowing to reach long propagation times are given in the Appendices A and B. The validity of the GPE for describing the nonequilibrium dynamics of this many-body body system is a priori not obvious. However, previous investigations focusing on disorder averaged observables such as ∆x have shown good agreement with advanced many-body treatment such as the multi-configurational time-dependent Hartree method for bosons (MCTDHB) [49][50][51] at least for short times.

III. SCENARIOS FOR THE EXPANSION
We consider now three different scenarios for the expansion of a BEC in a disorder potential (Fig. 2). They allow to disentangle the influence of particle-particle interactions and disorder on the non-equilibrium dynamics and localization of the wavepacket during the expansion.
In scenario A, we let the BEC freely expand up to t = t D . The switch-on time of the disorder t D is chosen such that the initial particle-particle interaction energy is almost completely converted into kinetic energy. After the time t D we also switch off interactions, i.e. we set t D = t I in Eq. (3) and Eq. (4). As numerical value we choose t D = 20t 0 . At this point the instantaneous inter- In scenario B, the disorder potential is present during the entire expansion time (i.e. t D = 0) but the particleparticle interactions (or nonlinearity) are switched off at t I = 20t 0 . In this case, E int (t) = 0.25E total , i.e. a significantly larger fraction of the total energy is neglected at the point of switch-off. Finally, in scenario C the disorder is also present from the start of the expansion (t D = 0) as in scenario B but the particle-particle interactions are never switched off during the entire simulation (i.e. t I → ∞). It is the latter scenario that most closely represents the experiment [27].

A. Scenario A
This scenario has been previously used to describe analytically [48] the experimental observation in Ref. [27]. The underlying assumption is that during the initial expansion stage (up to t = t D ) where most of the interaction energy is transferred into kinetic energy the disorder can be neglected while for later t > t D = t I the remaining interaction energy can be neglected. In the present case switching off the interaction, i.e. setting g = 0, at t D = t I amounts to neglecting approximately 4% of the total energy E total . Snapshots of the simulated evolved real-space and momentum density are shown in Fig. 3.
The analytical description [48,[52][53][54] of the emerging asymptotic real-space density for this scenario uses the fact that each momentum space component, i.e. plane wave of the wavepacket amplitudeψ(k, t I ) in the by now non-interacting system localizes with a k (or energy) dependent localization length with E = k 2 /2. Eq. 20 follows from the first Born approximation which applies to high energies and weak disorder, i.e. E V 2 D and V D σ D 2 k m (kσ D ) 1/2 . Consequently, the localization length diverges in the limit of large energies in this approximation. For a Gaussian disorder potential λ grows exponentially with E, for a speckle potential it diverges for E > 2 As the expanding BEC n(x, t) n(k, t) acquires a high-energy tail during the conversion of E int into kinetic energy, a fraction of the BEC executes a rapid quasi-ballistic expansion thereby counteracting localization. The quasi-ballistic expansion is clearly seen in the rapidly moving tails of the real space density [ Fig. 3 (a)]. The asymptotic density distribution can be derived from where n[λ(E), x] is the density of the Anderson localized state with localization length λ(E) [53,54]. A(k, E) denotes the spectral function (see Appendix C) and ψ (k, t I ) 2 is the initial momentum distribution. In this case the "initial" state of the localization scenario is the momentum distribution at the time the particle-particle interaction is switched off. The wavepacket at t = t I can be calculated analytically as [52,55] |ψ(k, t I )| 2 = 3ξ shown in Fig. 3 (b) where we have chosen t I such that the stationary momentum distribution [Eq. (22)] is reached to a high degree of accuracy. Eq. [21] predicts very flat tails reaching out to infinity [ Fig. 3 (a)]. Comparing with the numerically calculated densities as a function of time we observe very good agreement for t > t I . No fitting parameter is involved here. The numerically obtained non-stationary densities feature tails with local maxima above the analytical prediction. The local peaks move ballistically while at the same time successively spreading out. These tails contain high momentum components that only weakly scatter at the disorder potential. In fact, for the speckle potential, the scattering probability vanishes. Consequently, the local peaks are even more pronounced for the speckle potential than for the Gaussian correlated disorder (Fig. 4). Obviously, taking the width ∆x of the entire wavepacket as measure for the localization fails to capture the multi-scale nature of the expansion process that features simultaneously quasi-ballistic expansion of the outer tails and a fairly localized central region which may, possibly, display (sub)diffusive spreading.
It is worth noting that the "initial" moment distribu- at the moment of switching-off the interparticle interaction [Eq. (22)] is only weakly perturbed by the subsequent propagation in the disorder [ Fig. 3 (b)]. This, at first sight, surprising result can be explained as follows. For each disorder realization the wavepacket can be expanded in Anderson modes, the eigenstates of the linear Schrödinger equation with the disorder potential, asψ where φ n (k) are the Anderson modes in momentum space and a n are the expansion coefficients at time t I . Note that Eq. (23) is valid only in the non-interacting regime of scenario A and B [i.e. g(t) = 0 for t > t I ]. The momentum density is then given by For large differences in energy |E m − E n |, the Anderson modes will have a small overlap with each other in momentum space such that the corresponding off-diagonal contributions in the sum can be neglected. For small (but non-vanishing) differences |E m − E n | and large E m , E n the overlap still will be negligible since the width of the Anderson modes in momentum space ] for large momenta. The disorder average further reduces the offdiagonal terms such that we observe very good agreement between n(k, t I ) and n(k, t) for t = 100t 0 and t = 500t 0 , see Fig. 3 (b). The tails [see the log-plot in the inset of Fig. 3 (b)] represent the momentum distribution of the Anderson modes, i.e. the diagonal contributions to Eq. (24).

B. Scenario B
This scenario distinguishes itself from scenario A by the immediate presence of the disorder potential from t = 0 on (t D = 0). As in scenario A, the particle-particle interactions are switched off at t I , [g(t > t I ) = 0]. This scenario allows to investigate the influence of disorder on the initial expansion process, in particular, on the momentum distribution that provides the initial condition for the subsequent non-interacting Anderson scenario for t > t I . The presence of disorder during the time interval 0 < t < t I gives rise to new features not present in scenario A: In presence of the disorder potential the conversion of interaction energy into kinetic energy cannot proceed unhindered. Instead, the initial superfluid flow experiences friction in form of excitation of phonons as soon as momenta above the Landau velocity (or momentum) k L become available [56]. Those momenta then scatter inelastically at the disorder potential. Correspondingly, we observe for t < t I the emergence of minima in the momentum distribution around ±k L [50] where k L = µ(t I ), and an enhanced density for smaller momenta [ Fig. 5 (b)]. Once the pair interaction is switched off for t > t I , the disorder averaged n(x, t) n(k, t) momentum density remains essentially unchanged as was the case in scenario A. It should be noted that in scenario B the sudden switch-off of the interaction amounts to neglecting quite a large fraction (∼ 25%) of the initial total energy E tot (compared to 4 % in scenario A) reflecting the suppression of energy conversion by disorder. Unlike scenario A, scenario B cannot be treated fully analytically because the momentum distribution at t = t I is not known analytically. However, we observe only little deviations in the numerically obtained time-dependent densities as compared to scenario A confirming that the modified initial momentum distribution has indeed a minor effect on the shape of the densities at later times as assumed in [27,48,52] [ Fig. 5 (a)]. Likewise, the present results suggest that inserting the numerically found momentum distribution at t = t I , ψ (k, t I ) 2 , in Eq. 21 allows for a reasonably accurate description of the asymptotic real-space density n(x). n(x, t) n(k, t)

C. Scenario C
Within the mean-field approximation, scenario C is the one which should most closely describe the experimental situation: Disorder is present from the start and particleparticle interactions are present throughout the propagation. Moreover, comparing scenario C with B which differ only by the presence (or absence) of interactions beyond t = t I allows to identify the effect of interactions on Anderson localization for long times. At first glance, the time-dependent density n(x, t) for scenario C (Fig. 6) closely resembles that of the two previous scenarios A and B. The density features simultaneously a strongly peaked region and rapidly ballistically moving tails. There is, however, a unique feature only present in this case (see also [57] Fig. 5 (a)]. The effect can also be observed in the momentum distribution which continuously broadens as a function of time [ Fig. 6 (b)]. It is important to note that detecting this feature requires disorder averages with sufficiently large ensemble sizes to prevent blurring the signal by large fluctuations present. Observing them in the experiment may pose a considerable challenge as both a sufficient spatial resolution of the density as well as a large enough ensemble would be needed. In presence of particle-particle interactions Anderson modes are no longer eigenstates of the full Hamiltonian. Consequently, the expansion of the wavepacket in Anderson modes features now time-dependent expansion coefficients whose time evolution is governed by i∂ t a n (t) = g m1,m2,m3 where Eq. (27) describes the excitation of new modes not already excited initially. Taking a m (t = t I ) as initial conditions, it is the dynamics built into Eq. (27) that governs scenario C while for scenario B i∂ t a n = 0 holds for all t > t I . The redistribution of mode excitations is determined by a subtle interplay between spatial overlap of the modes and their energies. In a disordered system modes with near degenerate eigenenergies must have small (spatial) overlap in order to suppress level repulsion. Two mechanisms should play a major role: the near resonant excitation of modes with center of mass outside the initial spread of the wavepacket at t I [38][39][40], and a non-resonant excitation of modes with center of mass close to the origin but with larger energies and localization lengths. The spreading of the momentum distribution [ Fig. 6 (b)] indicates the importance of the latter mechanism.

IV. SUB-DIFFUSIVE SPREADING AND LOCALIZATION
In this section we focus on observables that have been used in the past to characterize spreading or localization of the expanding BEC. The behavior of the width ∆x [Eq. (15)] has gained considerable attention [30,[38][39][40]. For discrete systems described by a nonlinear Schrödinger equation, sub-diffusive spreading, i.e. ∆x(t) ∝ t α with α < 1/2 has been predicted and observed. The exponent is predicted to be α = 1/4 if the  ) for scenario A and C for the Gaussian correlated disorder. The labels "inner" and"outer" reflect the contributions to ∆x 2 (t) when only the contributions from the inner region x ∈ [−3000 : 3000] x0 or only from the outer region |x| > 3000x0 are included. These partial second moments are additive ∆x 2 inner + ∆x 2 outer = ∆x 2 total .
interaction energy is larger than the mean level spacing within the wavepacket, and α = 1/6 if the interaction energy is smaller than the mean level spacing [38,39]. These results for the discrete system cannot be directly transcribed to the continuum system investigated here. One important difference between the discrete system and the continuous system is that the energy spectrum of the discrete system is bounded both from below and from above while for the continuum system it is only bounded from below. Therefore, self-trapping, where the interaction energy is larger than the upper bound of the spectrum, cannot occur in the continuum system. As ∆x(t) measures the overall width of the wavepacket, it is particularly sensitive to the ballistically moving tails. Due to the presence of these tails in the asymptotic densities, even in scenario A without the interaction present (g = 0) the width as a function of time, ∆x(t), grows, and does not saturate (Fig. 7). To disentangle the contributions to the increase of ∆x 2 due to the spreading of the central part of the wavepacket from that for the quasi-ballistic growth of the tails we calculate the contributions to ∆x 2 (t) from both within the interval x ∈ [−3000 : 3000] x 0 and outside of it (Fig. 8). While in scenario A the central part of the wavepacket does not grow significantly after t ≈ 250t 0 , indicating the onset of localization, the wavepacket spreads continuously for scenario C showing no sign of localization. The outer part containing the tails originating from the high momentum components, however, shows a similar increase for both scenarios and strongly contributes to the total spread of the wavepacket for later times. Scenario B shows qualitatively the same behavior as scenario A. Interestingly, we observe in scenarios B and C for the speckle potential larger deviations of ∆x(t) for individual realizations of the disorder as indicated by the larger standard deviation [ Fig. 7 (b)]. This is due to the fact that the peak values of the disorder potential and, consequently, the total energy fluctuates more strongly between individual realizations of the disorder for the speckle than for the Gaussian correlated disorder. Since in scenario A the wavepacket has already substantially spread when disorder is switched on, the standard deviation of ∆x is reduced in this case compared to that for scenarios B and C. For the width in momentum space ∆k(t) we observe, in contrast to the behavior of ∆x(t), for g = 0 (scenarios A and B) a saturated momentum distribution which translates into ∆k(t) becoming stationary as a function of time (Fig. 9). In order to characterize the (sub)diffusive spreading and to compare with discrete models we extract estimates for the time-dependent exponents for all scenarios (Fig. 10). Within the time interval accessible by our simulation which exceeds the observation time in the experiment [27] by a factor 7, convergence to a well-defined exponent cannot (yet) be observed. For both scenarios A and B in which the long-time expansion proceeds in the absence of interactions, α(t) is still increasing for large t. Only for scenario C, α(t) is only weakly varying which may indicate an approach to a saturated value. If so, its value is, however, much higher than previously observed for discrete models and with α ≈ 0.4−0.5 close to the transition from the sub-diffusive to the diffusive regime.
The large values of α observed within all scenarios reflect the presence of fast quasi-ballistically moving tails absent in discrete models. It should be emphasized, however, that for an unambiguous determination of the longtime limit α(t → ∞) still much longer propagation times would be needed. These would, in turn, also require exceedingly large box sizes as the rapidly moving tails will eventually reach the wall of the box. A conclusive numerical test appears currently not yet computationally feasible. Generally, the extracted value of α cannot differentiate between the expansion of the outer region caused by the high-momentum tails and the expansion of the central region caused by particle-particle interactions (Fig. 8). As the ballistic expansion is generic for expanding BECs and not specific to the presence of particle-particle interactions, ∆x(t) appears not to be a suitable measure to probe the influence of interactions on Anderson localization of an expanding BEC. Since the long-time spread ∆x(t) is currently neither experimentally nor computationally accessible and, moreover, strongly influenced by the ballistic tails of the wavepacket, it is useful to focus on alternative and more local measures for localization that address the shape evolution on smaller length scales close to the center of n(x, t). While n(x, t) globally does not display a simple exponential shape near the center of the distribution one can extract an effective exponential slope within a spatially restricted interval. To this end, we fit the densities in Sec. III to Eq. (30) in different intervals. We address in the following only the scenarios B and C because the shape of the central peak of the density of these two is similar and considerably different from scenario A. Subsequently we focus on speckle disorder for comparison with the experimental results of Ref. [27]. The results for Gaussian disorder are, however, very similar. In previous works [27,52] it was estimated that the effective localization length L loc for scenario A is given by with λ given by Eq. (20) and the maximum energy E = E max = 1/2ξ 2 after the free expansion until t = t D = t I . Since within scenario B the momentum and energy distribution reached at t = t I differs from that of scenario A it is a priori not obvious that Eq. (31) should hold for scenario B as well. Nevertheless, since the high-momentum cut-off after the early expansion remains essentially unchanged by the friction of the wavepacket in the disorder potential, Eq. (31) provides still a useful estimate for L loc for scenario B. Since the shape of the numerically simulated central peak of the density is not characterized by a single exponential, the fitted value for L loc sensitively depends on the interval of x coordinates included in the fit. We probe this dependence for scenario B for which the density distribution near the maximum converges (see Fig. 5). The inverse slope L loc increases by more than a factor 3 when its interval is shifted from close to the central peak [35 : 100] in units of x 0 further out to the wings [650 : 1000]x 0 (Fig. 11). The approximately predicted value of L loc ≈ 242x 0 [Eq. (31)] is found near an intermediate interval [360 : 650]x 0 . In view of this variation of the local slope a comparison with the experimental data of Ref. [27] is not unambiguous. In Fig. 12 we compare the experimental data for the evolution of L loc as a function of t with our numerical prediction for scenarios B and C in two different extraction intervals, close to the peak [35 : 100]x 0 [ Fig. 12 (a)] and at larger distance [100 : 360]x 0 [ Fig. 12 (b)]. For the latter case, we find reasonable agreement with experimental data up to the largest times observed in the experiment (t = 68t 0 or 2 seconds). The significance of this agreement should, however, be viewed with caution. Apart from the uncertainty in the interval over which the experimental decay of the density has been determined the experimental data does not cover the following evolution where scenarios B and C begin to diverge (t > 100t 0 ) and where scenario B yields, indeed, a stationary value of L loc while scenario C does not [ Fig. 12 (a) and (b)]. The particle-particle interactions included in scenario C render the slope of the central peak time-dependent yielding a monotonically growing L loc . Extrapolation on a t −1 plot (Fig. 13) suggests that L loc (t) may diverge. To test this assumption we fit the data points for scenario C to a function f (1/t) = a (1/t) b for the extraction interval x ∈ [100 : 360] x 0 (purple dashed lines in Fig. 13). We find a finite value of b ≈ 0.35 suggesting that L loc does, indeed, diverge. As for ∆x(t), also for L loc (t) a definite conclusion on the asymptotic behavior cannot be drawn from the present simulation. The point to be noted is that the apparent divergences for ∆x(t) and L loc (t) occur on different length scales and are due to fundamentally different processes. While the continued growth of ∆x(t) is mostly unrelated to particle-particle interactions in the late phase of the expansion but originates from excitation of energetically high-lying Anderson modes with very large (or even diverging) localization length in the early phase of the expansion, the continued growth of L loc (inverse slope of the central peak) is a true interaction effect at late times. This observation strongly suggests that particle-particle interactions, indeed, destroy Anderson localization. We emphasize that the difference between the non-interacting and interacting case (scenario B and scenario C) becomes visible only after relatively long times (t = 150t 0 which corresponds to 4.5s) beyond the observation time in [27]. Observation of this effect in future experiments may pose a challenge. The present results provide quantitative predictions for the time scale on which the role of interactions on Anderson localization in the expansion of a BEC become experimentally observable.

V. SUMMARY AND CONCLUSIONS
We have investigated the expansion of a onedimensional (1D) Bose-Einstein condensate (BEC) in Gaussian correlated as well as speckle disorder closely following the experimental setup in Ref. [27]. For the theoretical and numerical description of the system we have used the Gross-Pitaevskii equation (GPE). We addressed the role of interactions (i.e. of nonlinearities) on Anderson localization by studying three different scenarios that allow to disentangle effects due to disorder and interactions during the expansion process. The conversion of interaction energy into kinetic energy in the early phase of the expansion results in a broad energy distribution of the wavepacket and fast moving tails. In the presence of a disorder landscape this broad excitation spectrum corresponds to the excitation of energetically high-lying Anderson modes with large (or even diverging) localization lengths. The presence of these ballistic components renders one of the frequently involved indicators for Anderson localization, the saturation of the width ∆x of a wavepacket as a function of time, largely inapplicable, even in the absence of interactions at large times. When both disorder and interactions are simultaneously present, ∆x grows even faster with time. In this case also the momentum distribution continuously grows indicating that new Anderson modes with larger energies and localization lengths are excited due to particle-particle interactions. We have investigated whether the width obeys a power law ∆x ∝ t α with α < 1/2 predicted analytically and verified numerically for the discrete system [38,39]. Our numerical results show that even in the absence of interactions, α(t) still grows at late times rendering convergence to a definite exponent α < 1/2 for subdiffusive spreading unlikely due to the rapidly expanding tails originating from high-lying Anderson modes. These have not been observed in previous works and may be absent in discrete systems. When both disorder and interactions are present, the time dependence of α(t) becomes markedly weaker. The saturation value, if it exists, is however much higher than for discrete systems and close to the diffusive limit α ≈ 0.5. On a shorter distance scale the inverse slope of the approximately exponential decay of the central peak in the real-space particle density can provide a useful local measure for the localization length L loc . In the absence of interactions at late times the numerically extracted values for L loc become (approximately) timeindependent and are of the order of the analytic prediction for the localization length for the Anderson modes with energies close to the maximum energy realized in the early phase of the expansion. However, with interactions present, the inverse slope and, thus, L loc continues to increase up to long times. We show that the extracted values of L loc is strongly dependent on the interval over which the slope is measured. Differences between the interacting and the non-interacting case become significant only after t = 150t 0 which corresponds to 4.5s with the parameters from the experiment [27]. This time is substantially longer than the previously accessed time and sets a benchmark for future experiments that can address the role of interactions on Anderson localization in ultra-cold atoms. The present numerical results suggest that, indeed, interactions destroy localization also on this length and time scale. SFB042-VICOM, FWF-SFB049-NEXTlite, FWF Doctoral College Solids4Fun (W1243) as well as the IMPRS-APS.

Appendix A: Numerical method for time propagation
For the numerical time integration of the GPE which belongs to the class of nonlinear Schrödinger equations of the type we employ high-order adaptive multiplicative operator splitting methods. These s-stage exponential splitting methods for the integration of Eq. (A1) use multiplicative combinations of the partial flows φ A (t, ψ) and φ B (t, ψ). For a single step (0, ψ 0 ) → (dt, ψ 1 ) with timestep t = dt, this reads where the coefficients a j , b j , j = 1 . . . s are determined according to the requirement that a preset order of consistency is obtained [58]. The dot in the brackets symbolizes that during each application of a new partial flow propagator φ the most current realization of the state vector is used. φ B • φ A represents the subsequent application of the partial flow φ B after the partial flow φ A .
The system in our case contains two vector fields, the kinetic energy part and the potential energy consisting of both the nonlinearity and the disorder potential, of different stiffness. The stiff flow is associated with the kinetic energy. If these are treated separately the resulting subproblems can typically be integrated with more efficient schemes. For the GPE the kinetic part can be integrated efficiently after (pseudo)spectral space discretization by exponentiation of a diagonal matrix while the nonlinear part allows an exact integration of the resulting ordinary differential equation in real space. Thus, the numerical effort effectively reduces to transformations between real and Fourier space which can be implemented with low cost also in a parallel environment. In the present case we use an equidistant grid and employ fast Fourier transforms within the package FFTW [59]. The efficiency of the time discretization can be improved if high-order time propagators are employed. These provide a more efficient approximation if the exact solution is sufficiently regular. In the present case we have used a fourth-order propagator as introduced in Ref. [60,61]. Adaptive choice of the time steps yields a further reduction in the computational effort. As a basis for time-step adaption asymptotically correct a posteriori estimators of the local time-stepping error are required. We can choose from four classes of local error estimators which have different advantages depending on the context in which they are applied (for a review see e.g. Ref. [60]  and are based on reusing a number of evaluations from the basic integrator. In our simulations, we have used an embedded pair of orders 4(5) referred to as Emb 5/4 AK (ii) in the collection [61]. Based on a local error estimator the step size is adapted such that the prescribed tolerance is expected to be satisfied in the following step. If dt old denotes the current step-size the next step-size dt new in an order p method is predicted as (see [63]) dt new = dt old ·min α max , max α min , α tol P(dt old ) 1 p+1 , (A3) where α = 0.8, α min = 0.25, α max = 4.0, and P denotes an estimate of the local error. This strategy incorporates safety factors to avoid an oscillating and unstable behavior. In Fig. 14 (b) we present the adaptive time step as a function of time for a disorder realization for different tolerances. For the calculations in the present work we use a tolerance of tol = 10 −10 corresponding to an average time step of 5.5 × 10 −4 t 0 . For this time step the energy is conserved up to ten significant digits which is comparable to the results using a tolerance of tol = 10 −11 . For speckle potentials the time step is about a factor of two larger for the same error tolerance than for Gaussian disorder potentials. The error in energy conservation is of the same order. Note that the present propagator allows for a time step which is ≈ 10 2 larger than the one required by the leap frog propagator used previously [ 49,50]. Only through this speed up long time propagations as presented in this paper become possible. For the three different tolerance levels shown in Fig. 14 (a) the width of the wavepacket differs from each other by much less than one standard deviation of the ensemble average (see also the discussion in App. B).

Appendix B: Convergence tests
In this section we discuss the convergence of our results with respect to spatial discretization and box size. The results presented here correspond to individual realizations of the disorder and do not contain any ensemble averaging. We show representative results only for scenario C and the Gaussian disorder potential as it is the numerically most challenging one. The results for speckle potentials are similar. Perfect agreement of each individual wavefunction for different error tolerances, spatial discretizations, and box sizes cannot be expected because of the underlying chaoticity of the GPE [49]. However, we regard the results as converged if the differences are smaller than the variance within the ensemble of disorder realizations. Comparing different spatial discretizations we observe that they agree very well within the error bars (one standard deviation) of the ensemble average (Fig. 15). For the results presented in the main text we use a spatial discretization of dx ≈ 0.0293 x 0 . Comparing different box sizes [−x max : x max ] we see a rather large deviation after long times for the smallest box (Fig. 16). All larger boxes give results within the standard deviation of the ensemble average. The influence of reflections off the walls at −x max , x max are thus negligible. For the calculations in this work we use x max = 2.4 × 10 4 . With the spatial discretization given above this corresponds to a grid of 1638400 points. In this appendix we briefly review the key ingredients required to calculate the asymptotic densities (Eq. 21). The density of an Anderson localized state with fixed energy E is given by the integral [64,65] n(λ(E), x) = π 2 γ(E) 8 gives the probability amplitude to scatter into a state with energy E at momentum k. The self-energy Σ can be calculated in first order Born approximation yielding for the imaginary part The real part Σ describing the deviation from the dispersion relation of a free particle is assumed to be small [53] and is neglected.