{'fraction_non_alphanumeric': 0.037886784067551804, 'fraction_numerical': 0.014249208377312372, 'mean_word_length': 4.277484608619173, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 8, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Wave chaos as signature for depletion of a Bose-Einstein condensate

We study the expansion of repulsively interacting Bose-Einstein condensates (BECs) in shallow one-dimensional potentials. We show for these systems that the onset of wave chaos in the Gross-Pitaevskii equation (GPE), i.e. the onset of exponential separation in Hilbert space of two nearby condensate wave functions, can be used as indication for the onset of depletion of the BEC and the occupation of excited modes within a many-body description. Comparison between the multiconfigurational time-dependent Hartree for bosons (MCTDHB) method and the GPE reveals a close correspondence between the many-body effect of depletion and the mean-field effect of wave chaos for a wide range of single-particle external potentials. In the regime of wave chaos the GPE fails to account for the fine-scale quantum fluctuations because many-body effects beyond the validity of the GPE are non-negligible. Surprisingly, despite the failure of the GPE to account for the depletion, coarse grained expectation values of the single-particle density such as the overall width of the atomic cloud agree very well with the many-body simulations. The time dependent depletion of the condensate could be investigated experimentally, e.g., via decay of coherence of the expanding atom cloud.


I. INTRODUCTION
The workhorse for describing the non-equilibrium dynamics of Bose-Einstein condensates (BECs) of ultracold gases is the Gross-Pitaevskii equation (GPE) (for a review see e.g. Ref. 1,2). Replacing the true many-body wave function by a single-particle orbital for the macroscopically occupied condensate (particle number N 1) results in an equation of motion that belongs to the class of nonlinear Schrödinger equations (NLSE). The GPE provides an appropriate starting point to investigate the underlying many-body system on the mean-field level. Effects beyond the GPE have been observed in BECs, for example in optical lattices with deep wells and small occupation numbers per site. 3 Other finite-number condensate effects include the demonstration of atom-number squeezing 4-6 and of Josephson junctions in a double well. [7][8][9] Meanwhile, progress has been made in exploring the time-dependent many-bosons Schrödinger equation. One approach is the multiconfigurational time-dependent Hartree for bosons (MCTDHB) method which is a numerically efficient and, in principle, exact method for the time-dependent many-body problem. [10][11][12] In practice, limitations are imposed by the finite yet large number of configurations (millions) and orbitals (tens) that can be handled. We investigate repulsively interacting BECs after release into shallow one-dimensional (1D) potentials. The Bose gas is dilute and, initially, practically all particles are in one single-particle state. The external potential is weak compared to the single-particle energy. Comparison be-tween the MCTDHB method and the GPE for the expansion of the BEC provides detailed insights to what extend the GPE is capable of describing the condensate dynamics and may be capable of mimicking excitations out of the condensate state. One case in point is our recent observation of true (physical) wave chaos in the GPE, 13 as opposed to numerical chaos 14 due to discretizations which has been exploited to study e.g. thermalization in the Bose-Hubbard system at the mean-field level. 15 Two wave functions nearby in Hilbert space are exponentially separating from each other, as measured by the L 2 norm. Chaotic wave dynamics within the GPE is a mathematical consequence of the non-integrability resulting from the interplay between the external (one-body) potential and the nonlinearity which replaces the inter-particle interactions. Its physical implications are, however, less clear as the original many-body Schrödinger equation is strictly linear and, thus, regular and non-chaotic. Previously, a connection between chaotic dynamics within the GPE and growth in the number of non-condensed particles has been made for time-dependent external driving which can be seen as a source of energy. 16 In the present study the external potentials are time independent such that the total energy is conserved. While wave chaos is likely associated with instabilities known for dynamics in periodic potentials (see e.g. Ref. 17) it is a much more general effect since it occurs in a broad variety of potentials ranging from harmonic oscillators with defects to periodic and disordered potentials. The aim of the present paper is to shed light on the physical meaning of wave chaos in the GPE for time evolution of BECs. By comparing the mean-field GPE and many-body MCT-DHB dynamics we show that the build-up of random fluctuations within the GPE is directly related to the depletion of the condensate, i.e., the population of excited orbitals out of the condensate state. In this respect, the mean-field description surprisingly contains some information on many-body effects. Chaos on the mean-field level may open the door to identify depletion and thermalization and the comparison with MCTDHB may shed light on how such a thermalization process takes place in a quantum many-body system. The outline of the paper is as follows. After first introducing the system under investigation in Sec. II we briefly review the mean-field GPE and the many-body MCT-DHB method and identify relevant observables (Sec. III). The initial state whose dynamics we study upon release from the initial trapping is discussed in Sec. IV. We present numerical results for the dynamics in Sec. V followed by conclusions and remarks.

II. SYSTEM UNDER INVESTIGATION
We consider in the following a system of N bosons interacting via a pseudo-potential which captures the scattering dynamics of the real interaction potential in the limit of small wave numbers k → 0. In a 1D system with tight transverse harmonic confinement with oscillator frequency ω r the pseudo-potential is given by the contact where a s is the 3D scattering length provided that a s /mω r such that the scattering can still be regarded as a 3D process. 18 The dynamics of the bosonic system is then determined by the many-body Hamiltonian (in second quantization) The field operators fulfill the usual commutation rules for bosons. We study in the following the expansion of a Bose gas that is initially trapped also longitudinally (i.e. in the direction of expansion) by a harmonic potential with frequency ω 0 . These initial conditions serve to define characteristic scales for length, time, and energy. We use the units l 0 = /mω 0 for length, t 0 = 1/ω 0 for time, and e 0 = ω 0 for energy. For a trap with ω r = 2π × 70Hz and ω 0 = 2π × 5.4Hz used in a recent experiment on Anderson localization 19 our units take on the numerical values l 0 = 4.6µm and t 0 = 29.47ms. Upon release from the trap, the particles move in an external potential V (x) which we specify in the following -10 The initial (t < 0) harmonic trapping potential of radial frequency ωr and longitudinal frequency ω0 (black line). At t = 0 the longitudinal harmonic trap is switched off (the radial trap is preserved). Simultaneously, either a periodic (dashed red) or a disorder (blue solid) potential is switched on. The system expands for t ≥ 0 in the 1D external potential. The length scale l0 corresponds to 4.6µm.
to be a periodic potential of the form with l = 0.54811l 0 (corresponding to l ≈ 5.8ξ with ξ = / √ 4mµ being the healing length and µ the chemical potential after release from the trap for N = 1.2 × 10 4 87 Rb atoms) and varying potential amplitude V A . The periodic potential is realized in experiments by crossed laser beams in linear polarization along the same axis. For a realistic laser wave length tuned out of resonance with the 87 Rb 5S → 5P transition, λ L = 810nm, the above potential period of l corresponds to two linearly polarized crossed beams enclosing an angle of θ ≈ 0.1π. Alternatively, we also consider Gaussian correlated dis- , and the correlation length σ. The interplay between the inter-particle interaction and the external potential plays a key role for chaotic dynamics resulting from non-integrability. We investigate in the following the dynamics of the expanding Bose gas in the mean-field approximation within the GPE and compare to the corresponding many-body dynamics within the MCTDHB method.

A. Gross-Pitaevskii equation
In the mean-field approximation the existence of a macroscopic occupation of one state is assumed such that the expectation value ψ (x, t) = ψ(x, t) takes on finite values and can be treated as the classical field describing the dynamics of the BEC. Further, requiring that the expectation value of the product of four field operators factorizes i.e., assuming second-order coherence, one arrives together with andĤ from Eq. 2 at the GPE Renormalization of the particle density to dx|ψ(x, t)| 2 = 1 leads to the explicit dependence of the nonlinearity on the particle number N : Since N enters as a parameter into the effective nonlinearity, the GPE predicts the same dynamics for different N as long as the product g 1D N is kept constant.
In the limit N → ∞ with g 1D N = const. the (timeindependent) GPE is expected to give exact results for the many-body system (at least for the ground state in three dimensions 20 ). The numerical simulations presented in the following are performed for N = 1.2 × 10 4 ultracold 87 Rb atoms in a cigar-shaped trap with frequency ω r (see Sec. II). The associated nonlinearity is Note that the rather high numerical value of g 0 is due to the explicit inclusion of the number of particles N and does not contradict the assumption of weak interactions. Nevertheless, the interaction strength is sufficiently strong such that in the presence of an external potential depletion and fragmentation of the condensate may occur.
To propagate the GPE, we use a finite element discrete variable representation (DVR) to treat the spatial discretization (see e.g. Ref. 21,22). The propagation in time is performed by a second-order difference propagator (for details see Ref. 13 and references therein).

B. Multi-configuration time-dependent Hartree for bosons (MCTDHB) method
The MCTHDB method 10,11 allows one to describe many-body effects beyond the mean-field description for the condensate. Briefly, the many-body wave function is taken as a linear combination of time-dependent permanents where | n; t corresponds to states with occupation numbers n = (n 1 , ..., n M ) and M is the number of singleparticle orbitals. The sum runs over all sets of occupation numbers { n} which fulfill N = M i=1 n i . In the limit M → ∞ the ansatz Eq. 9 gives the exact manybody wave function. MCTDHB efficiently exploits the fact that ultracold atoms may occupy only few orbitals above the condensate state. By dynamically changing the expansion amplitudes C n (t) and the orbitals {Φ k ( r, t)}, even large many-body systems can be treated accurately. MCTDHB involves the solution of coupled linear differential equations in C n (t) and coupled nonlinear differential equations in {Φ k ( r, t)}. The MCTDHB equations of motion reduce in the case of M = 1 to the GPE Eq. 7 [with nonlinearity g 1D (N − 1), the difference between N and N − 1 can be neglected in the limit of large N ]. Within the MCTDHB method kinetic operators are treated via a fast Fourier transform which is equivalent to an exponential DVR. 23 The nonlinear differential equations in the orbitals are propagated via a 5 th order Runge-Kutta algorithm. The linear differential equations for the amplitudes C n (t) are propagated via a shortiterative Lanczos algorithm. The propagations are parallelized using openMP and MPI. The MCTDHB orbital equations of motion are stiff differential equations for the interaction strength considered in this work. The applicability of the Runge-Kutta algorithm has been checked by comparing to an especially adapted integrator 24,25 ZVODE which relies on the Gear-type backwards differentiation formula for stiff ordinary differential equations and gives the same results as the faster Runge-Kutta algorithm. Numerical checks give a very good agreement between the initial and the backwards propagated density per particle. The difference is of the order of 10 −4 and less (in units of l −1 0 ).

C. Observables
The simplest and most important benchmark observable for a comparison between the mean-field and the many-body dynamics is the single-particle density. Within the MCTDHB method the density is given by where the elements ρ m,n (t) are readily accessible as a combination of the amplitudes C n (t) and the correspond-ing occupation numbers contained in n. Upon diagonalization of Eq. 10 the natural orbitals Φ NO i (x, t) and their occupation numbers n NO i (t) are obtained as: In the presence of a BEC the occupation of one state has to be "macroscopic" 26 (of order N ). In the following we denote this condensate state as Φ NO i=1 (x, t) and its occupation as n NO i=1 (t). All other states Φ NO i (x, t) with i > 1 are referred to as excited states. The Fourier spectrum is obtained by Fourier transforming the orbitals Φ n (x, t) to giveΦ n (k, t). In the limit of a long-time expansion of the BEC in free space when the initial interaction energy is converted into kinetic energy,ρ(k, t) corresponds to the experimentally observed momentum distribution. Within the GPE,ρ(k, t) is given by the absolute square of the Fourier transform of the condensate wave function |ψ(k, t)| 2 . We utilize coherence as measured by the normalized twoparticle correlation function 27,28 to analyze the breakdown of the GPE on the length scales of the random fluctuations which develop in the wave function in the regime of wave chaos. In g (2) the reduced two-body density matrix enters. For a fully second-order coherent system g (2) fulfills |g (2) Within the GPE the reduced two-body density matrix is a product of onebody wave functions. Thus, |g (2) | = 1 for all times, i.e., full second-order coherence is a generic feature of the GPE. In the many-body case the departure of |g (2) | from unity gives a measure for how well the system is described by a single-orbital product state and how correlated (|g (2) | > 1) or anticorrelated (|g (2) | < 1) the measurement of two coordinates is. (Anti-)Correlation indicates the degree of fragmentation in the system. As a measure for wave chaos, i.e. the build-up of random local fluctuations on the length scale comparable to that of the external potential, we have introduced 13 the Lyapunov exponent characterizing the exponential increase of the distance in Hilbert space of two initially nearby GPE wave functions ψ 1,2 (x, t). The distance is measured by the L 2 norm The distance function takes on values d (2) ∈ [0, 2] and is 1 for orthogonal wave functions. In terms of d (2) , the Lyapunov exponent indicating chaos is given by d (2) is invariant for unitary time propagation of linear systems: if ψ 1 (x, t) and ψ 2 (x, t) would be solutions of the linear Schrödinger equation, d (2) would be constant. Similarly, inserting instead of ψ 1 (x, t) and ψ 2 (x, t) two many-body wave functions Ψ 1 (x 1 , ..., x N , t) and Ψ 2 (x 1 , ..., x N , t), and integrating over all spatial coordinates gives a constant of motion. By contrast, construction of a reduced one-particle wave function from an initial N -body state of system 1 by (see e.g. Ref. 29) where |Ψ 1 (N ) denotes a many-body state with N particles leads to a many-body measure analog to the d (2) (t) function that is not conserved as a function of time. This can be seen by inserting Eq. 17 into Eq. 15, taking the time derivative of d (2) , and using the many-body Hamiltonian. In the time derivative of d (2) contributions originating from the kinetic energy and the external potential V (x) cancel, while contributions from the interaction term lead to d[d (2) (t)]/dt = 0. The tracing out of unobserved degrees of freedom leads to the violation of the distance conserving evolution. In the particular case of the GPE the nonlinearity present can cause exponential divergence of d (2) . It is now of key importance to relate the behavior of the d (2) function within the GPE to properties of the unitary time evolution of the underlying many-body system, as described by the MCTDHB method. The working hypothesis is that the random fluctuations developing within the GPE are the signature for its failure to properly account for the depletion of the condensate state, i.e. excitation of the BEC during expansion in an external potential. In turn, within the MCTDHB method the population of all natural orbitals beyond that describing the condensate should grow. While at t = 0 the MCTDHB method and the GPE closely agree which each other with only one natural orbital occupied, n NO 1 (0) = n NO 1 (0)/N ≈ 1 (see next Sec. IV), with increasing time all other occupation numbersn NO i (i > 1) should increase. In the following we study the dynamics of N = 10 3 to 10 5 particles for which the ground state densities closely agree with each other (see Fig. 2). For N = 10 4 and N = 10 5 only M = 2 orbitals allow a numerically feasible number of configurations (N c = 10 4 +1 to N c = 10 5 + 1, respectively). Already adding one more orbital (M = 3) leads to a configuration size of N c = 50, 015, 001 for N = 10 4 which may be at the border of feasibility and requires a massive parallelization over a large number of processors. The system with N = 10 5 and M = 3 resulting in N c ≈ 5 × 10 9 is out of reach for the current implementation of the MCTDHB method. For N = 10 3 a number of orbitals up to M = 3 is numerically feasible and allows to quantify the effect of adding one more orbital to the case M = 2. Due to the numerical limitations we focus on the early stages of the depletion process when the depletion is still relatively weakn NO 1 (t) ≥ 0.95. As a quantitative measure for the depletion process we introduce the state entropy for a general many-body state wheren NO i (t) = n NO i (t)/N and n NO i (t) are the occupation numbers of natural orbitals obtained by diagonalizing the reduced one-body density at each time (see Eq. 11). For the initial conditions used in the present study we have S N (t) ≥ S N (0) ≈ 0. S N (t) measures the entropy of the occupation numbers of the many-body system along distinct quantum states. Within the GPE, where n NO 1 (t) = 1 and n NO i>1 (t) = 0, S N (t) remains strictly zero. Deviations of S N (t) from zero within the manybody theory thus mark the breakdown of the mean-field theory within the GPE. In the following we will focus on the time evolution of Eq. 18 and investigate the time scale of depletion, t d , defined by the occurrence of an abrupt change of S N from S N ≈ 0 to S N > 0. We associate this observable with the onset of exponential growth, t e , of d (2) (t) within the GPE. t e is determined from the crossing-over point between the free-space expansion behavior of d (2) (t) to the exponential increase. t d is implicitly dependent on N through the degree of coherence of the condensate. The larger N , the smaller the depletion (n NO i , i ≥ 2) at the same time. Consequently, the depletion time is size dependent t d (N ) (see below).

IV. THE INITIAL STATE
The initial state of the bosonic gas corresponds to the ground state of the harmonic trap. For this ground state the GPE with nonlinearity g 0 predicts a BEC in the Thomas-Fermi regime. Applying both the MCTDHB method and the GPE to the same system requires a careful choice of system parameters, in particular the particle number N . While the validity of the GPE calls for the limit of large N → ∞, such a case is numerically prohibitive for the MCTDHB expansion Eq. 9. Since the GPE results are invariant for varying N but fixed g 0 = g 1D N we adjust the particle number such as to remain in the Thomas-Fermi limit of the longitudinally trapped BEC (see Fig. 2). In such a way it is assured that discrepancies between the GPE and the MCTDHB method during the time evolution are not caused by incompatible initial conditions. The ground state of an interacting system of bosons trapped by a harmonic potential is governed by three length scales: the characteristic length of the harmonic trap l 0 , the mean inter-particle distance r s = n −1 with n the particle number per unit length, and l δ = 2 mg1D a measure for the zero-point fluctuations (or anti-correlation length) of the repulsive two-body deltafunction interactions of strength g 1D . The regimes obtained range from a non-interacting "Gaussian" shaped BEC, over a Thomas-Fermi BEC, to a strongly interacting fermionized Tonks-Girardeau gas. 30 The presence or absence of a BEC is determined by the ratio referred to as the Lieb-Lininger parameter. 31 If l δ is much larger than the inter-particle spacing r s the particles are allowed to occupy the same state and form a BEC. The condition for the presence of a BEC thus is: In order to distinguish between a Gaussian and a Thomas-Fermi BEC the harmonic oscillator length l 0 must be included. The regimes are controlled by the parameter 30 In our case α 1 (α = g 0 /N with g 0 from Eq. 8), i.e. l 0 l δ . In this case the system is in the Thomas-Fermi regime if in addition l 0 r s which is equivalent to the statement that the Thomas-Fermi length L TF = 2µ/mω 2 0 corresponding to the width of the BEC in the Thomas-Fermi regime is much larger than l 0 . The condition l 0 r s implies N α −1 for the Thomas-Fermi limit to hold. 30 For all systems with N ≥ 1000 in Tab. I the criteria N α −1 and γ 1 are well fulfilled and, indeed, the many-body ground state density takes on the Thomas-Fermi shape (see Fig. 2). The density is practically indistinguishable from the GPE prediction. For comparison, we also show in Fig. 2 a system with N = 100 for which the criterion of a Thomas-Fermi BEC is only marginally fulfilled because γ ≈ 1 and deviations become apparent. The requirement of large N places a severe limit on the number of orbitals that allow for a numerically feasible configuration space. Convergence in the orbital number is controlled by the occupancy n NO M of the least occupied state. While for the ground state calculations n NO M is sufficiently low, we expect this number to rapidly increase during expansion since strong depletion may occur which will spread over a large number of orbitals. We, therefore, expect only the onset of depletion to be quantitatively reliable while the occupation numbers of excited orbitals can be considered to be an indication of the excitation process as the orbital expansion ceases to converge (M > 3 time-dependent orbitals would be needed) with increasing propagation time.

V. NUMERICAL RESULTS
We first consider the expansion of the BEC initially formed within the harmonic trap (Fig. 2) which is, after the release, subject to a periodic potential (Eq. 3) with l = 0.54811l 0 and V A = 0.2e with e the total energy per particle. After the release an explosion-like process takes place: the interaction energy is rapidly transformed into kinetic energy. In free space the cloud expands keeping its Thomas-Fermi shape with the characteristic length increasing in time. 32 This process is modified by the presence of the periodic potential. Practically immediately the density is modulated by standing waves with the same spatial periodicity as the potential. The local maxima of the density coincide with the local minima of the potential and lead to an increase of kinetic and interaction energy at cost of potential energy. As soon as the Fourier spectrum is sufficiently broad, inelastic processes set in. At that point a fraction of particles have acquired a momentum k > k L = mv L with v L the Landau velocity which leads to excitations of phonons, i.e. friction of superfluid flow. For a homogeneous system v L is given by v L = ng1D m with n the particle density and g 1D the interaction strength (prefactor of contact interaction). An estimate for the Landau velocity in case of an inhomogeneous system can be derived from the chemical potential µ which contains an average over the small scale fluctuations of the density. For a homogeneous system µ = ng 1D . We thus use v L = µ m to determine the numerical value of k L . At t ≈ 3t 0 the width ∆k of the Fourier spectrum is approximately as large as k L and we observe the development of strong density modulations [ Fig. 3 (a)]. These spatial density modulations go hand in hand with reduced density in the Fourier spectrum near ±k L since those particles lose their momentum by phonon excitations [Fig. 3 (b)]. Friction leads to the separation of a strongly fluctuating central part of the density from its fast tails [ Fig. 3 (a)]. The tails expand nearly freely and are modulated by the potential. We point out that this process is fully accounted for within the GPE (i.e., the system remains condensed) since it gives practically the same density and spectrum for t = 3t 0 as MCTDHB in Fig. 3. For longer times we have previously observed for this system signatures of wave chaos: 13 two nearby effective onebody wave functions ψ 1 (x, t) and ψ 2 (x, t) (with initially large overlap) propagated by the GPE become orthogonal to each other after an exponential increase in distance in Hilbert space [see Fig. 4 (c)]. This orthogonality results from the build-up of random local fluctuations in the wave functions on length scales comparable to the period of the potential. It is now instructive to compare the growth in d (2) within the mean-field description with the growth of S N (t) (or depletion) within the MCTDHB method which the GPE cannot follow. For vanishing potential V (x) = 0 we find that the explosion-like expansion with a rapid transformation of interaction energy to kinetic energy does not lead to depletion of the condensate [ Fig. 4 (b) dashed line]. The GPE accounts for the expansion dynamics since S N (t) remains approximately zero as a function of time. For vanishing potential the GPE is integrable 33 such that d (2) saturates after a short universal increase [ Fig. 4 (c) dashed line]. For periodic potentials (Fig. 4) we find over a wide range of potential strengths (0.04e ≤ V A ≤ 0.2e) that the exponential separation on the mean-field level and the depletion on the many-body level correlate well with each other. To extract the depletion time t d and the rate of depletion η we fit the curves to functions of the form with fit parameters c, a, and t d (in Fig. 4 t d as well as t e are marked for V A = 0.2e). Θ in the above equation is the Heaviside step function. We introduce the depletion rate η as The depletion rate η is equal to the slope of S f N (t) at t = t d + t 0 , i.e. after the abrupt increase of S f N (t) at t d . η has the same dimension as the Lyapunov exponent λ. An increasing Lyapunov exponent λ with increasing V A goes hand in hand with an increasing η (Fig. 4). Up to a constant numerical factor (≈ 10) η closely follows λ function of two close wave functions within the GPE expanding in the periodic potentials of (a). The parameter VA is indicated next to the curves in units of the energy per particle e. The manybody system corresponds to a particle number of N = 10 4 and an orbital number of M = 2. The onset of exponential growth te as well as the onset of depletion t d are marked by arrows. t d corresponds to a time of ≈ 60ms. Note that SN is strictly zero within the GPE.
as a function of V A and even a threshold-like behavior is present (Fig. 5, compare also to Ref. 13 Fig.7). We note that both t d and η are sensitive to the specifics of the fit function Eq. 22 which results in an uncertainty of the fit.
Since we are interested in the qualitative behavior here, the error bars in t d and η are of no concern. The variation of t d from the onset of exponential divergence of d (2) in the GPE is comparable to the variation of t d within the MCTHB method for different particle numbers N [compare Fig. 4 (b), (c) with Fig. 6]. t d increases with increasing N and approaches t d = 6.48t 0 which we obtain from the GPE by associating the onset of chaos t e Another important example is propagation in a disorder potential. We compare the results for periodic potentials to propagation in a disorder potential with parameters for which previously Anderson localization has been observed. 19 Within the GPE for V A = 0.2e and correlation length σ = 0.7ξ we obtain on a time interval of t ∈ [6,15]t 0 an exponential increase with practically the same Lyapunov exponent as for the periodic potential with V A = 0.2e and l = 0.54811l 0 (not shown). Only initially the slope is reduced [ Fig. 6 (c)]. Accordingly, we find within the many-body system that the curves for S N (t) follow each other with a reduced rate for the disordered system [ Fig. 6 (b)]. We point out that our results indicate a destruction of the BEC as signified by the occupation of excited states during expansion in disorder potentials in the regime of Anderson localization. The onset of depletion of the condensate is mirrored in the fine scale oscillations of the density (Fig. 7). Substantial deviations of the density within the GPE from the density obtained within the MCTDHB method emerge at different instants of time for different particle numbers. For N = 10 4 deviations in the local fluctuations of the density emerge at t ≈ 4t 0 (see Fig. 7) monitored by S N (t) > 0. The occupation numbers are n NO 1 ≈ 0.96 and n NO 2 ≈ 0.04. While the condensed part [given by n NO 1 (t)|Φ NO n1 (x, t)| 2 ] still closely follows the GPE prediction |ψ(x, t)| 2 , the total density ρ(x) shows smoothing of the local fluctuations near the center. This smoothing is due to excited atoms whose density partially fills in the local minima. For the system with N = 10 5 the picture is very similar except that the occupation of the excited state is lower at t = 4t 0 , n NO 2 ≈ 0.003 instead of n NO 2 ≈ 0.04. Spatially strongly localized structures emerge which the GPE overestimates as the transition into excited orbitals is forbidden. The initially spatially localized excitations spread over the entire system with increasing time. One can expect the fine scale structure of the density of the full many-body system to strongly differ from the prediction of the GPE. The discrepancies in the particle density go hand in hand with the breakdown of coherence as measured by the normalized two-particle correlation function (Fig. 8): g (2) directly quantifies the degree of second-order coherence in the system 27 Color online) Two-particle normalized correlation function g (2) (x1, x2, x1, x2, t = 4t0) for N = 10 4 bosons as well as the density of the second natural orbital |φ NO 2 (x, t = 4t0)| 2 (inset) obtained from MCTDHB. White corresponds to no correlations and full coherence, red to correlation and blue to anti-correlation, i.e. reduced coherence. Note that g (2) (x1, x2, x1, x2; t) is a real function (see Eq. 14) and is equal to unity within the GPE. second natural orbital) the two-particle coherence is lost, i.e. g (2) strongly differs from 1. The deviation of g (2) from unity indicates that the many-body state is no longer representable by a product of a single complex-valued function, i.e. the GPE ceases to be a valid description. This is a fingerprint of the emerging fragmentation of the many-body system. Roughly speaking, the pattern of g (2) (x 1 , x 2 , x 1 , x 2 ; t) is similar to the functional form of |φ NO 2 (x 1 , t) × φ NO 2 (x 2 , t)|. For longer time intervals our MCTDHB calculations indicate a destruction of the condensate or at least a strong fragmentation. At t = 20t 0 , e.g, the occupation of both orbitals is approximately 50% indicating that many more orbitals would be required for convergence. Nevertheless, current experiments indicate remarkable agreement with the prediction of the GPE for large scale observables such as the width of the atom cloud or the average position (see e.g. Ref. 19,[35][36][37]. For the width where x n = dxρ(x)x n , we previously observed that within the GPE it is independent of wave chaos: two close wave functions ψ 1 (x, t) and ψ 2 (x, t) develop random local fluctuations and become orthogonal. The width, however, for both ψ 1 (x, t) and ψ 2 (x, t) agrees. If we now compare the prediction for the width within the GPE and within the MCTDHB method, we observe the same trend. Despite its failure to account for the state entropy ( Fig. 4 and Fig. 6) and the coherence properties (Fig. 8), the GPE remains predictive in describing the expansion of a BEC in external potentials on longer time scales for coarse-grained observables such as the width and the average position, long after the random fluctuations prevent the prediction of fine scale structures in ρ(x). Up to now local small-scale fluctuations have not been investigated experimentally because of the difficulty of (sub) µm resolution. While the fine scale structures of the wave function within MCTDHB have not fully converged for the small number of orbitals (M ≤ 3 included in the simulation), the coarse-grained distribution remains essentially unchanged compared to the GPE [ Fig. 9 (a)] and for varying N and M . We thus expect that the time dependence of the width of the full many-body system is well accounted for by the curves in Fig. 9 (a). The same excellent agreement we observe for the average over momenta k 2 as accessible in time-of-flight experiments [see Fig. 9 (b)]. The average over k 2 is determined via and is proportional to the kinetic energy per particle. The GPE thus reproduces the mean kinetic energy of a highly excited system despite its failure to account for breakdown of coherence, fragmentation, and small scale fluctuations. The latter observation indicates that thermalization may be within the realm of the GPE despite its failure to account for two-body scattering which is key to any thermalization process.

VI. CONCLUSIONS
By comparative simulations invoking the Gross-Pitaevskii equation (GPE) and the multiconfigurational time dependent Hartree for bosons (MCTDHB) method we have uncovered a direct connection between wave chaos in the GPE and depletion of the occupation of a BEC during expansion in the presence of weak ex-ternal 1D potentials. We have checked that this connection holds for a large class of external potentials including a harmonic potential with short-ranged perturbation, an aperiodic potential with incommensurate frequencies, and disordered and periodic potentials explicitly discussed in this paper. This connection has farreaching consequences: while the depletion and fragmentation process is an intrinsic many-body effect outside the realm of the GPE, the mean-field theory allows one to monitor its onset through the development of random local fluctuations. The measure for the random local fluctuations, d (2) , can be used to delimit the applicability of the GPE to approximate the many-body dynamics. On the many-body level the depletion process manifests itself through the loss of coherence as measured by deviations of g (2) from unity. The connection between wave chaos and depletion indicates that a mean-field description may allow to identify and describe relaxation and thermalization processes (see e.g. Ref. [38][39][40][41], when the system reaches equilibrium for coarse-grained quantities while non-universal fluctuations about the equilibrium persist. The depletion process, the onset of which we have investigated, can be experimentally studied provided a sufficient spatial resolution is achieved. Observables include higher-order coherence, i.e. deviations of g (2) from unity as measured e.g. in Ref. 34. It would be of considerable interest to verify experimentally our predictions by exploring the fine-scale fluctuations and coherence properties of expanding BECs in external potentials and thus gain deeper insight into the involved many-body effects.