Density Resolved Wave Packet Spreading in Disordered Gross-Pitaevskii Lattices

We perform novel energy and norm density resolved wave packet spreading studies in the disordered Gross-Pitaevskii (GP) lattice to confine energy density fluctuations. We map the locations of GP regimes of weak and strong chaos subdiffusive spreading in the 2D density control parameter space and observe strong chaos spreading over several decades. We obtain a renormalization of the ground state due to disorder, which allows for a new disorder-induced phase of disconnected insulating puddles of matter due to Lifshits tails. Inside this Lifshits phase, the wave packet spreading is substantially slowed down.


Introduction
Disorder is inevitable naturally in all materials due to the impurities or defects caused by external fields. Interacting wave dynamics in disordered and incommensurate lattice structures are actively studied due to its complex properties [1]. Interacting waves can be modeled with nonlinear terms approximating true quantum many-body physics of interacting bosons. In the absence of interactions, Anderson localization (AL) leads to exponential localization of eigenstates and a coherent halt of wave packet spreading [2]. The complete suppression of wave propagation has been manifested by a bevy of experimental observation; including localization of light waves [3], photonic crystals [4], sound waves [5], microwaves [6], and atomic matter waves [7,8]. In the presence of nonlinear wave interaction terms, delocalization can arise and ultimately lead to chaotic dynamics, which destroys Anderson localization through incoherent spreading [9][10][11]. Notably, this phenomenon was also studied experimentally with ultracold atomic gases [12]. As interesting is the fact that experimental efforts are limited by the time of atomic condensate control, which -assuming a natural dimensionless time scale of order t = 1 -yields the largest realizable times of the order of t ∼ 10 4 (the typical natural time scale is of the order of 10ms, and condensate control usually extends up to 10s). At these experimental time horizons, the observation of the onset of incoherent spreading was possible, but a quantitative assessment of this process was not. Theoretical simulations can achieve t ∼ 10 8 − 10 9 [1], and in special settings of unitary maps with discrete-time quantum walks are reaching unprecedented times t ∼ 10 12 [13]. These numbers demonstrate the rare opportunity for computational studies being the superior testbed of the first choice.
Let us recap nonlinear wave packet spreading in simple terms. At the initial time t = 0, a wave packet is assumed to have a compact distribution of finite norm A and energy H, which extends over L lattice sites. Deprived of its nonlinear terms, the system manifests AL through exponentially localized eigenstates, of which roughly L are excited by the wave packet. Transforming into normal mode space yields a harmonic oscillator equation for each AL eigenstate, with the nonlinear terms inducing a short-range coupling between them. Assume that this dynamical system will be nonintegrable, characterized by nonzero Lyapunov coefficients, and evolving chaotically in time. The consequent phase decoherence of the normal modes removes the basis of existence for AL, and normal modes in the boundary layer at the edges of the wave packet will be incoherently excited. The wave packet will spread, and L 2 (t) ∼ t α increase in time. The assumption of complete dephasing of all AL normal modes yields the strong chaos regime α s = 1/2 for nonlinearity originating from two-body interactions [11]. However, the computationally tested asymptotic weak chaos regime was observed to yield α w = 1/3. This result can be derived assuming that the probability of a normal mode being resonant and chaotic is proportional to energy and/or norm density in the wave packet [11] (notably this assumption results in dependence of both α w and α s on the lattice dimension and different choices of N -body interactions [1]). Whether the observed weak chaos spreading is asymptotic or will slow down, has been a topic of debates and discussions [14], with still no answer in sight. What was confirmed in computations, is the potentially long-lasting inter-mediate strong chaos regime [15] -but notably only for systems with one integral of motion (energy) [15].
The spreading wave packet is assumed to thermalize on time scales shorter than the one on which it spreads. For systems conserving energy only, the energy density h(t) = H/L(t) corresponds to some inverse temperature β(t), and while spreading h(t → ∞) → 0 and thus β(t → ∞) → ∞, so the packet spreads and cools down. A thermal Gross-Pitaevskii (GP) wave packet however conserves, in addition, the total norm A and must be characterized by two energy h(t) and norm a(t) = A/L(t) densities which are related to some inverse temperature β(t) and chemical potential µ(t) [16]. The spreading dynamics then correspond to moving along a line in the density parameter space {a, h} which connects an initial point {a 0 , h 0 } with the origin. At variance to systems with only one integral of motion, a GP lattice supports non-Gibbs phases [17], and hence the outcome will depend on the chosen line, including heating upon spreading and possibly reaching infinite temperatures at finite values of L. Computational studies of GP wave packet dynamics [15] involved disorder averaging which was controlling a 0 but not h 0 , therefore averaging over fans of lines in the density parameter space.
In this work, we will unfold the density-resolved dynamics, which allows us to finally identify a clean strong chaos regime, and the potential slowing down Lifshits phase regimes. We observe strong chaos and map strong and weak chaos in the density parameter space, including a localization regime coined Lifshits phase (LP) due to a disorder-induced ground state renormalization.

Model definition
We consider the disordered Gross-Pitaevskii chain on N sites with Hamiltonian where ψ = √ n e iφ are complex scalars, and the integer enumerates the lattice sites. is a quenched uncorrelated on-site random potential taken to be uniformly distributed within a box of size W : W is a measure of the disorder strength, J the tunneling amplitude between neighboring sites, and g > 0 is the nonlinearity parameter resulting from e.g. the repulsive two-body interaction between atoms. Energy (and the inverse of time) can be measured in units of J which leaves us with two parameters: W, g. Uniform rescaling of the norm |ψ | 2 is used to tune the nonlinearity parameter into g = 1.
For g = 0, the system (1) is reduced to an eigenvalue problem using ψ (t) = Ψ e −iλt : λΨ = Ψ − (Ψ +1 + Ψ −1 ). It follows |λ ν | ≤ 2 + W/2, with all eigenvectors Ψ ,ν being exponentially localized in space for any W = 0 [2]. The localization length ξ(λ ν ) takes its largest value ξ 0 ≈ 96/W 2 for λ ν = 0 [18]. The typical core size of a corresponding eigenvector at energy λ = 0 fluctuates with an average of the order of V loc ≈ 3ξ 0 [11]. This follows from the numerical observation that the participation number of such a normalized eigenstate P ν = |Ψ | 4 ≈ 1.5ξ ν , whose random amplitude fluctuations result in another factor of 2 for the size of the state (which accounts for a site amplitude being either large or small with equal probability). The localization volume hosts V loc eigenstates with amplitudes of order V That group of significantly overlapping eigenstates has V loc eigenvalues which are separated from each other by the average level spacing d(W ) = (4 + W )/V loc .
With J ≡ g ≡ 1 the final equations of motionψ = ∂H/∂iψ * read The above dynamics conserves the total norm A ≡ n and the total energy H. The partition function where β and µ are Lagrange multipliers associated with the total energy and the total norm, respectively. Let us recap the equilibrium properties of the GP model assuming nonzero overall densities h = H/N and a = A/N . For the ordered case W = 0 the microcanonical ground state line h = −2a + a 2 /2 equals the grand-canonical zero temperature β = ∞ line, and the line h = a 2 corresponds to the infinite temperature line β = 0 [17] (see Fig.1). Pairs of densities in the Gibbs range −2a + a 2 /2 ≤ h ≤ a 2 are addressable by pairs of positive inverse temperature and chemical potential {β, µ}. The entire non-Gibbs density range h > a 2 is not captured by a positive temperature, while negative temperature assumptions lead to a divergence of partition functions, and microcanonical dynamics show strong deviations from expected ergodic behavior including self-trapping [17,19]. The Gibbs-nonGibbs separation line in the density space was recently shown to persist for entire classes of generalized GP lattice equations as well as their Bose-Hubbard quantum counterparts, for any lattice dimension, and in the presence of disorder [20]. Remarkably the addition of disorder for (1) leaves the infinite temperature line h = a 2 invariant [16].

Ground state renormalization
The zero-temperature ground state line of the ordered case h = −2a + a 2 /2 is renormalized in the presence of disorder. This happens in the regime of small norm density (i.e. weak nonlinearity) a < 1 due to the presence of Lifshits states which are sparsely distributed AL eigenstates with eigenvalues close to the bottom of the AL spectrum, i.e. their distance from the bottom ∆ λ = λ + 2 + W/2 1. Such Lifshits states exist due to rare disorder fluctuations with l + W/2 < ∆ λ over a simply connected chain segment of length L = π/ √ 2∆ λ . The average distance between such regions d L ≈ (W/∆ λ ) L . As a result, one can expect a set of disjoint puddles of norm distribution in real space for small norm density a. Note also that for any finite system the ground state is bounded by h = −(2 + W/2)a + a 2 /2 which is generated by the disorder realization = −W/2.
Contrary, in the large norm density limit (i.e. for strong nonlinearity) the ground state correction becomes weak since the nonlinear terms a 2 /2 are of leading order and disorder has a minor impact.
In order to numerically compute the ground state, we note that ψ can be gauged into real variables as all the phases φ = φ to minimize the Hamiltonian (1). The remaining task is to minimize a real function H for real variables ψ for a given disorder realization. We choose an initial set of ψ under the constraint N a = ψ 2 . We define a window of l w = 5 adjacent sites and minimize the energy varying the amplitudes on these adjacent sites using the Nelder-Mead simplex algorithm [21]. As the algorithm changes the total norm in general, we perform a homogeneous renormalization of all amplitudes ψ to restore the required norm density a. We then shift the window by one lattice site and repeat the procedure, until the whole lattice with N sites has been covered by minimization windows. The procedure is repeated 10 times, after which full convergence is obtained. The chemical potential is defined through local relations and yields a ratio of the standard deviation to mean which is less than 10 −3 , indicating the quality of our ground state computation. Finally, we repeat the procedure for 100 different disorder realizations and compute the average ground state energy density h and its standard deviation. The result is shown as red solid circles in Fig.1 with their standard deviation for N = 1000, and W = 4.

Initial conditions
We prepare a wave packet on L ≈ V loc consecutive sites in the center of a disordered lattice (see table 1). For a 2 /2 − 2a ≤ h ≤ a 2 /2 + 2a, we choose an initial state with a homogeneous norm distribution ψ = √ ae iφ . We fix the phase differences We then adjust the phase on one of the sites to tune the total energy such that H = Lh (and disregard disorder realizations for which the adjustment can not be realized).
For h > a 2 /2 + 2a or h < a 2 /2 − 2a, we use the ground state renormalization method, replacing the original energy with |H−Lh| and strictly varying only the amplitudes of the wave function on the L sites of the wave packet. With that, we can prepare density resolved wave packets which are characterized by a pair of initial density values {a 0 , h 0 }, and a corresponding point in the phase diagram Fig.1. A spreading wave packet is characterized by a pair of timedependent densities a(t) and h(t) moving along a straight line connecting the initial phase diagram point with the origin.
Parametrizing the line as h(t) = ca(t) and assuming a(t) 1/g, we approximate the partition function by its interaction-free limit g = 0 and derive h = − µ β (µ 2 − 4) −1/2 + 1 β and a = 1 β (µ 2 − 4) −1/2 for the ordered case. The line parametrization then finally yields β = − 2c a(4−c 2 ) and µ = − 4+c 2 2c . We can expect a number of qualitatively different spreading outcomes for the disordered case depending on the choice of the initial density values. If h 0 > 0 and the initial point is in the Gibbs phase, the origin-connecting line will cross the infinite temperature line. Therefore the wave packet will heat up to infinite temperatures, enter the non-Gibbs phase, and is expected to show features of self-trapping, fragmentation, and condensation. For h 0 = 0 and c = 0 it follows h(t) = 0, the temperature will gradually increase and will reach infinite values at infinite times. If −2a 0 + a 2 0 /2 < h 0 < 0 (i.e. c < 0), the wave packet may heat up to some finite temperature, but will asymptotically gradually cool down and reach potentially zero temperature upon approaching the origin at infinite times. Finally, if h 0 < −2a 0 + a 2 0 /2, the wave packet can evolve in the Lifshits phase.

Computational details
We integrate Eq.2 using a symplectic method SBAB 2 [22,23] for several different disorder realizations with time step size ∆t = 0.1 unless mentioned otherwise.
The wave spreading is characterized by two main ingredients. The second moment m 2 ≡ (l −l) 2 |ψ (t)| 2 /A 2 measures the width of the wave packet. Herel = l|ψ (t)| 2 /A 2 is the center of the wave packet. Participation number P ≡ A 2 / |ψ | 4 characterizes the bulk of the wave-packet and allows us to compute the compactness index C ≡ P 2 /m 2 . If C ∼ 1, then the wave packet is close to a thermal distribution, while C 1 indicates either fragmentation of the wave packet into a set of essentially disjoint pieces or a part of the wave packet staying localized with another part spreading [1]. In order to observe whether asymptotic subdiffusion m 2 ∼ t α is recorded, we compute The obtained function α(t) is additionally smoothened using a Hodrick-Prescott filter with a standard deviation of the output of less than 2% [24].

Observation of weak and strong chaos
In the absence of nonlinear interactions g = 0, a wave packet will evolve in time without appreciable spreading. We plot the evolution of its norm density versus space and time in Fig.2(a). After some short initial dynamics during which the field established exponentially localized tails in space, the wave packet evolution essentially halts, signaling Anderson localization.
To avoid non-Gibbs dynamics or a subsequent cooling of the wave packet, we choose h 0 = 0. A localized Anderson mode of the linear system will be interacting with d(W ) neighboring other modes in the presence of nonlinear interactions. For a > d, that interaction is expected to exhibit strong resonances and lead to efficient norm mixing among the participating modes [1]. Upon further wave packet spreading, the norm a(t) will decrease, and cross over into the asymptotic regime of weak chaos a < d. The evolution of the norm density of a wave packet in the weak chaos regime is plotted versus space and time in Fig.2(b). At variance to the case of Anderson localization Fig.2(a), the wave packet grows in size with increasing time.
The weak chaos regime is characterized by m 2 ∼ t 1/3 , as reported in Ref. [11]. We confirm these findings in our computations as shown in Fig.3. We launch wave packets for various values of W and a(t = 0) < d(W ). Both second moment m 2 (t) and participation number P (t) indicate subdiffusive spreading ( Fig.3 (a,b)). The compactness index C(t ∼ 10 8 ) ≈ 3, as expected for a thermalized system (Fig.3(c)). A quantitative evaluation of the time dependence of the exponent α(t) indicates its asymptotic convergence to α(t → ∞) ≈ 1/3.
In order to observe intermediate strong chaos m 2 ∼ t 1/2 , we simply need to increase the initial norm density a > d. The evolution of the norm density of a wave packet in the strong chaos regime is qualitatively similar to the weak chaos evolution Fig.2(b). The wave packet will spread in the strong chaos regime until some time at which a(t) ≈ d, which will mark the crossover to weak chaos. That crossover was observed for Klein-Gordon lattices in Ref. [15]. However the attempt to observe strong chaos and the crossover to weak chaos for the GP lattice failed [15]. The main reason is that the initial norm density a of the wave packet was controlled, but the initial energy density h was fluctuating when choosing and averaging over different disorder realizations [25,26]. As a result, qualitatively different regimes of strong chaos, self-trapping, and Lifshits phases were mixed into one curve. Interestingly the results in Ref. [11] were obtained for single-site excitations which did control both norm and energy, but the strong chaos regime was not observed because single-site excitations get self-trapped. Fig.4 shows our results when controlling both initial densities and performing density-resolved spreading studies. Both second moment m 2 (t) and participation number P (t) indicate subdiffusive spreading ( Fig.4 (a,b)). The compactness index C(t ∼ 10 8 ) ≈ 3, as expected for a thermalized system (Fig.4(c)). We observe strong chaos with α(t) reaching the value 1/2 and keeping this value for several decades in time, before slowly decaying, presumably to its asymptotic value 1/3.

Lifshits phase and self-trapping
In the Lifshits phase, the spreading of the wave packet is slowing down dramatically. The evolution of the norm density of a wave packet in the Lifshits phase is plotted versus space and time in Fig.2(c). It shows very little difference to the linear case of Anderson localization in Fig.2(a). The second moment grows slightly (blue line Fig.5(a)), while the participation number is essentially frozen (blue line Fig.5(b)). The compactness index drops with increasing time (blue line Fig.5(a)). The derivative α(t) (blue line Fig.5(d)) shows a slow increase.
In the non-Gibbs regime, the GP dynamics turn non-ergodic, with the field ψ l (t) forming strongly localized self-trapped large-amplitude excitations which are persisting on a background of delocalized waves [17,[27][28][29]. The self-trapped field part appears to condense in a way such that the remaining background field part can evolve at an infinite temperature in a Gibbs regime. Since the exact amount of the norm density in the background is not known, its wave packet spreading dynamics can not be precisely determined and classified as weak or strong chaos. The evolution of the norm density of a wave packet in the self-trapping regime is plotted versus space and time in Fig.2(d). We clearly observe a remaining strongly localized self-trapped part of the wave packet at the origin, which is much more localized than its Anderson localization counterpart in Fig.2(a). The background part of the wave packet instead spreads qualitatively similar to the weak (and strong) chaos case Fig.2(b). Consequently, the second moment of the wave packet grows in time due to the spreading of the background part (red line Fig.5(a)). At the same time, the participation number drops as time increases which indicates a slow but steady formation of self-trapped excitations on very few lattice sites (red line Fig.5(b)). The compactness index evolution (red line Fig.5(c)) confirms these findings very well. Finally, the derivative α(t) (red line Fig.5(d)) indicates that the background wave packet part may reach asymptotic weak chaos characteristics at times which are not accessible by our computational resources.
The wave packet spreading in both the Lifshits phase and the self-trapping regime is characterized by a substantial slowing down from the subdiffusive spreading as observed for weak and strong chaos. At the same time, the self-trapping enforces highly localized almost single site excitations to be formed, at variance to the Lifshits phase dynamics where the wave packet structure resembles the Anderson localization case.

Conclusion
We have analyzed energy and norm density resolved wave packet spreading studies in the disordered Gross-Pitaevskii (GP) lattice. We managed to confine energy density fluctuations, which were not controlled in previous studies. We mapped the locations of the GP regimes of weak and strong chaos sub-diffusive spreading in the two-dimensional density control parameter space and observed strong chaos spreading over several decades. A number of qualitatively different wave packet spreading outcomes were identified. Positive energy densities are doomed to bring the wave packet into a non-Gibbs regime with potential fragmentation of the packets into a self-trapped condensate part and an infinite temperature background capable of spreading infinitely. One of the intriguing quantities is the ratio of the norm in the two field components and its asymptotic time dependence. Will the self-trapped component take over the entire wave packet norm at large enough times, or will some finite remain in the infinite temperature background? Zero energy densities keep the wave packet in the Gibbs regime and may lead to the entire packet heating up to infinite temperatures upon infinite spreading. Negative energy densities in the Gibbs regime can bring the wave packet closer to the ground state, and therefore zero temperatures upon spreading. Finally, initializing the wave packet in the Lifshits regime shows strong suppression of subdiffusive spreading. But even in this case, we notice a speedup of the spreading process with some potential fragmentation of the wave packet. It appears that there are a number of interesting and hard open problems to be addressed in future work.