Spectral description of the dynamics of ultracold interacting bosons in disordered lattices

We study the dynamics of a nonlinear one-dimensional disordered system from a spectral point of view. The spectral entropy and the Lyapunov exponent are extracted from the short time dynamics, and shown to give a pertinent characterization of the different dynamical regimes. The chaotic and self-trapped regimes are governed by log-normal laws whose origin is traced to the exponential shape of the eigenstates of the linear problem. These quantities satisfy scaling laws depending on the initial state and explain the system behaviour at longer times.


Introduction
The motion of non-interacting particles in disordered lattices has been intensively studied over the last few decades. In one and two dimensions, and for a sufficient amount of disorder in three dimensions, it has been shown that the spreading of quantum wavepackets is suppressed, a phenomenon known as Anderson localization (AL) [1,2]. However, the celebrated Anderson model that leads to this prediction is a simplified model that, in particular, neglects particle-particle interactions, so a crucial question is: do interactions destroy (or, in contrast, enhance) AL?
A burst of interest in the subject was recently seen, driven by the demonstration that these questions could be studied experimentally and theoretically with an unprecedented degree of cleanness and precision using ultracold atoms. The study of the problem in the absence of interactions led to impressive results, concerning both AL [3][4][5][6][7][8] and the Anderson transition [9][10][11][12], observed in three-dimensional (3D) systems. The effects of interactions on ultracold bosons turned out to be very well modelled by mean-field approaches [13,14], in contrast to fermionic systems where such a simplification of the corresponding many-body problem is not possible. Moreover, the level of interactions can be easily controlled in an experiment using Feshbach resonances [15].
From the theoretical and numerical points of view, ultracold bosons in a one-dimensional (1D) optical disordered lattice can be described rather realistically by a simple generalization of the original Anderson model including a nonlinear term taking into account interactions in a mean-field description. Numerical studies suggest that the long-term motion of the particles is subdiffusive [16][17][18][19][20][21]. The same conclusions hold for the Aubry-André model, where disorder is simulated by a bichromatic lattice [22,23]. These studies pointed out, in particular, the central role of chaotic dynamics in the destruction of AL. A subdiffusive dynamics has recently been observed experimentally in the case of a bichromatic lattice [24].
From the numerical point of view, there are two important difficulties when one tries to characterize the dynamics resulting from the interplay of interactions and disorder. The first is that the timescale of subdiffusion is larger by several orders of magnitude than the timescale of the single-particle dynamics, implying very long computer calculations. The second is that the nonlinearity leads to a strong dependence on the initial conditions, which also makes it difficult to give a 'global' characterization of the different dynamical regimes. In previous works, we have shown the existence of scaling laws with respect to the width of the initial state, allowing, to some extent, such a global characterization [25], and demonstrated that these scaling laws are robust with respect to decoherence effects [26] (which can also destroy AL).
In this work, we tackle the problem of interacting ultracold bosons in a 1D disordered lattice using a spectral analysis that does not require such long computational times, as the information on the chaotic behaviour is 'inscribed' even in early times in the spectrum of the dynamics. This allows us to perform a more complete and precise study of the problem over a large range of parameters. The use of the spectral entropy has proved very useful in characterizing the dynamic behaviours, which is confirmed by comparing the information extracted from the spectral entropy to the Lyapunov exponent, another well-known measure of chaotic behaviour. Moreover, we show that such behaviours are described by 'log-normal' laws which can be scaled with respect to the initial conditions, and we propose a simple physical interpretation of these findings.

The model
We use the discrete nonlinear Schrödinger equation with diagonal disorder. It is essentially the same model used in previous works [25,26], so we only give an outline of it here.
The mean-field theory applied to ultracold bosons in an optical (ordered) lattice leads to the so-called Gross-Pitaevskii equation where φ is the macroscopic wavefunction (or order parameter) describing a Bose-Einstein condensate well below the critical temperature, V (x) = −V 0 cos (2k L x) is the optical potential and g 1D is the 1D coupling constant [27]. Tight-binding equations are obtained by projecting the wavefunction φ(x) = n c n (t)w n (x) onto the set of localized Wannier functions of the first band w n (x) associated with the nth lattice site where we kept only (symmetric) nearest-neighbours couplings, which is justified if the Wannier functions are strongly localized, that is, for large enough V 0 . According to the usual conventions, we measure distances in steps of the lattice and energies in units of the coupling constant of neighbour sites T = − dxw n (x)[ p 2 /2m + V (x)]w n+1 (x), 2 and as a consequence, times are measured in units ofh/T . The coefficient v n = dxw n (x)[ p 2 /2m + V (x)]w n (x)/T is called diagonal on-site energy and the effect of interactions is taken into account, in a mean-field approach, by adding the nonlinear term g |c n | 2 where the dimensionless interacting strength is In the absence of disorder, the on-site energy v n does not depend on the site index n and can thus be set to 0. According to Anderson [1], we introduce diagonal disorder by picking random on-site energies v n uniformly in an interval [−W/2, W/2]. For g = 0, equation (1) describes the standard Anderson model, and we shall call the corresponding eigenstates (eigenvalues) 'Anderson' eigenstates (eigenvalues). A few facts about it will be useful in what follows. The eigenstates (labelled by the index ν in the lowest band and associated with the eigenvalue ν ) are exponentially localized in the average c ν n ∼ exp (−|n|/l ν ) (the overbar indicates the averaging over realizations of the disorder) with a localization length l ν ( ν , W ) ∼ 96(1 − 2 ν /4)/W 2 [28]. For g = 0 the equation becomes nonlinear, and it is useful, if not strictly rigorous, to interpret the nonlinear term as a 'dynamical correction' to the onsite energy v NL n = g|c n | 2 . Previous works [16-19, 21, 25] provided evidence for the existence of three main dynamical regimes: for g W , AL is expected to survive for very long times, a regime that we shall call 'quasi-localized'. For g ∼ W , the nonlinear correction v NL n induces chaotic dynamics leading to subdiffusion and the destruction of AL. For g W , the very large nonlinear term v NL n decouples all sites whose populations are not nearly equal (even in the absence of the disorder), suppressing the diffusive behaviour and leading to another (nonlinear) type of localization, called self-trapping. This dynamics does not rely on quantum interference and is therefore very robust against external perturbations, including decoherence [26].
Our goal is to characterize the global dynamics based on a numerical integration restricted to a relatively short timescale. We thus study the evolution according to (1) of bosons in a 1D box containing L sites (typically, L = 101) and put an exponential absorber at each end of the box in order to prevent wavepacket reflection. The norm of the wavepacket is thus no longer conserved as soon as it 'touches' the borders, and we characterize the diffusive behaviour by calculating the survival probability p(t) = n |c n | 2 ; a value p(t) < 1 indicates that the packet has diffused outside the box. The use of absorbing conditions implies neglecting the possibility for a nonlinear excitation to spread over a large distance and then come back towards the centre of the system, but one can reasonably consider that the probability of such a phenomenon is relatively small, because, in most cases we are considering here, the dominant dynamics is known to be either (sub)diffusive or localized. Moreover, we have verified that doubling the size of our 'box' does not affect our results. As in [25], we restrict the analysis to initial wavepackets of the form 3 The typical behaviour of the survival probability is illustrated in figure 1, which represents p(g, t = 10 5 ) as a function of the interaction strength g for W = 4 and for various values of L 0 . The survival probability was averaged over typically 500 realizations of the disorder and of the initial phases θ n . Such long times are necessary in order to clearly provide evidence for the three different regimes mentioned above, which are then easily identifiable in figure 1: at low g, the survival probability is close to 1, corresponding to the quasi-localized regime in which AL survives. For intermediate values of g, the decrease of p indicates that the wavepacket has spread along the box and part of it has been absorbed at the borders, corresponding to the onset of a chaotic regime induced by the nonlinearity, which leads to the destruction of AL. For large values of g, the survival probability increases again and gets close to 1, indicating that the packet is self-trapped [16,17,25] and never touches the borders if it is thin enough. In section 3, we show that a spectral description of the dynamics allows a very precise characterization of the dynamics with much shorter computation times.

Characterizing chaotic dynamics with the spectral entropy
Spectral analysis is a very useful way of analysing a chaotic dynamics, either in classical, 'quantum' 4 or quantum nonlinear chaotic systems [29]. The spectral entropy is a measure of the 'richness' of a spectrum. Chaotic behaviours are associated with continuous spectra and thus to a high spectral entropy. Given a quantity M(t), its power spectrum is defined as The choice of the value of t max determines the resolution of the spectrum. In our case, we chose t max = 200, which is a very small value compared to the typical time of the dynamics we are interested in here (t ∼ 10 3 -10 4 ) but which will be shown to be sufficient to characterize the dynamics from a spectral point of view. As we do not expect excitations whose timescale is inferior to the tunnelling time (1 in our rescaled units), we set f max = 1. 5 From the power spectrum, one defines the spectral entropy as For a perfectly monochromatic signal The spectral entropy 'counts' the number of frequencies which are present in the signal, and is a good indicator of the chaoticity of the system [30]. In addition to a significant gain of computation time, the use of spectral analysis makes the system rather insensitive to the finite size of the box. We checked that for times t ∼ 200 the survival probability remains close to 1, which means that the influence of the absorbers is negligible. The spectral entropy obviously depends on the choice of the observable. A good observable in the present problem is the so-called 'participation number' P with respect to the Anderson eigenstates, which is defined, in the present case, as follows. For a given realization of disorder {v n }, we calculate the Anderson eigenstates in the Wannier basis, φ ν (x) = n d (ν) n w n (x) corresponding to an energy ν , which are solutions of (1) with g = 0 (the d (ν) n replacing the c n ). In the g = 0 case, equation (1) allows us to calculate the evolution of the wavepacket ψ(t) = c n (t)w n (x) under the action of both disorder and nonlinearity 6 . At any time, we can express the wavepacket in the Anderson eigenstates basis from which, trivially, q ν (t) = n d (ν) n c n (t). The participation number is then defined as If g = 0, the φ ν are the exact eigenstates of the problem, so that the populations |q ν | 2 are constant. In the non-interacting case g = 0, the |q ν | 2 evolve under the action of nonlinearity. The participation number roughly 'counts' the number of Anderson eigenstates participating significantly in the dynamics 7 and its time evolution thus reflects the excitation of Anderson eigenstates that were not initially populated. Figure 2 shows the spectral power S( f ) of the participation number P in the three interacting regimes: (a) g = 1 in the quasi-localized regime, (b) g = 100 in the chaotic regime and (c) g = 1000 in the self-trapping regime. For g = 1, the dynamics is very similar to the linear case: the populations of Anderson eigenstates practically do not evolve in time, as well as the participation number P and the spectrum is dominated by low frequencies. As a consequence, the spectral entropy, calculated according to (2), is relatively small: H = 10 −2 . For g = 100, Anderson eigenstates are strongly coupled by the nonlinear term: each pair of coupled states generates a Bohr frequency (shifted by the nonlinear correction), that is, In this regime, most Anderson eigenstates are coupled to each other, so that the power spectrum is almost flat (with important local fluctuations) at high frequencies, and the spectral entropy increases by almost an order of magnitude, H = 10 −1 , 6 We use a standard Crank-Nicholson scheme with time step 0.01 < dt < 0.1. 7 To see this intuitively, consider two limit cases: if only one ν = ν 0 Anderson eigenstate is populated, |q ν | 2 = δ ν,ν 0 and thus, from (4) with respect to the preceding case. For g = 1000 the wavepacket is self-trapped and only a relatively small number of Anderson eigenstates whose populations happen to be close enough can interact. The spectrum is therefore dominated by a finite number of frequencies and the spectral entropy is reduced, H = 2 × 10 −2 . However, although populations are stable in this case, quantum phases may evolve chaotically under the action of the nonlinearity [31]. Our definition of the spectral entropy from the participation number-which does not directly depend on the phases-excludes this 'phase dynamics' from the corresponding spectrum. We display in figure 3(a) the averaged spectral entropy as a function of the nonlinearity parameter g for different widths L 0 of the initial state. One clearly sees the crossover from the quasi-localized to the chaotic regime, signalled by a marked increase of H . The smaller the value of L 0 , the smaller the value of the crossover. This is easily understandable, as a more concentrated wavepacket leads to a stronger nonlinear term v NL . On the right side of the plot, one also sees, especially for low values of L 0 , the beginning of a decrease of H due to self-trapping.
It is interesting to compare the information obtained from the spectral entropy to another relevant quantity characterizing chaos, the Lyapunov exponent, which indicates how exponentially fast neighbour trajectories diverge. This quantity is usually defined for classical systems, but can be extended, with a little care, to quantum nonlinear systems [29]. We present in the appendix a method for calculating the Lyapunov exponent of a quantum trajectory defined by amplitudes c n (equation (1)). Figure 3(b) displays the Lyapunov exponent λ obtained with the same parameters as in figure 3(a). It displays a monotonic increase with the nonlinear parameter, even in the region where the spectral entropy decreases due to self-trapping, evidencing the presence of a regime of 'phase chaos' mentioned above. This clearly shows that H and λ provide different information on the dynamics of the system, and that H is clearly more adapted to distinguishing the three dynamical regimes.
In section 4, we will discuss the shape of these curves and their scaling properties and suggest a physical mechanism for explaining such properties.

The log-normal law and scaling
As we shall see below, log-normal functions are ubiquitous in the dynamics described in this work and so in various contexts which are not obviously related to each other. A first example can be found in figure 1, where each curve representing p(logg) can be fitted around its minimum with a rather good accuracy by an inverted Gaussian function. More generally, a log-normal function is defined as In physics, such a function appears more often as a 'log-normal distribution', related to the statistics of quantities which are a product of randomly distributed terms [32]. Log-normal statistics is very different from normal (Gaussian) statistics; for example, the most probable value of a log-normally distributed quantity is different from its average value. In the following, we shall not only consider log-normal curves that correspond to statistical distributions over the realizations of the disorder. Figure 1 displays the average of the survival probability as a function of the interacting strength, so does figure 3(a) for the spectral entropy. However, and more interestingly, we shall show that the statistical distributions of the spectral entropy and of the Lyapunov exponent are log-normal.
Let us thus study the distribution of these two quantities over the realizations of the disorder v n and of the initial phases θ n . Figure 4 A simple, heuristic model can give a hint on the origin of such shapes. The destruction of AL is due to the nonlinear coupling between Anderson eigenstates. Initially unpopulated eigenstates which, in the absence of nonlinearity, would never be populated can be excited thanks to a transfer of population induced by the nonlinearity. The goal of the simple model developed below is to characterize the statistical distribution of such excitations. Projecting the wavepacket in the Anderson eigenbasis (cf (3)), equation (1) then reads with n . Population exchanges are controlled by the overlap I of four coefficients. The population transfer from eigenstate µ to eigenstate ν depends on J 1 = I (µ, µ, ν, ν), J 2 = I (µ, µ, µ, ν) and J 3 = I (µ, ν, ν, ν), the term due to J 1 being, for instance, In order to evaluate the probability distribution of J 1 , we make the assumption that the corresponding Anderson eigenstates are exponentially localized with the same localization length ξ . The exponential localization is valid on average and the localization length is the same if the eigenstates have close enough eigenenergies. We thus write The overlap sum J 1 can then be written as where, for simplicity, we have supposed that the spatial distance between the eigenstates, l(µ, ν), is an integer. The most important term in (7) is e −2l(µ,ν)/ξ . In the limit ξ → 0: It is known that the inverse localization length = 1/ξ follows a normal distribution [28,33] we obtain for the distribution of values of the overlap J 1 where G = log [l(µ, ν) + 1] − 2l(µ, ν) 0 andσ 2 = 4l(µ, ν) 2 σ 2 . The overlap sum J 1 thus obeys a log-normal distribution.
We conjecture that overlap sums such as J 1 , controlling the coupling between Anderson eigenstates, in fact control the destruction of the AL [18]. Given that nonlinear excitations are directly linked to the values of overlap sums such as J 1 , as shown in (6), it is reasonable to think that the spectral entropy and the Lyapunov exponent, which evaluate the number of nonlinear excitations, have the same statistical properties as J 1 . The above heuristic argument undoubtedly presents various assumptions that are not rigorously justified, but it has the merit of providing evidence for the intimate relation between the log-normal distribution and the exponential localization of the Anderson eigenstates. The link between the exponential shape and the emergence of log-normal statistics has also been studied in the case of the conductance of disordered systems [28,34,35].
In previous works [25,26], we showed that suitable scaling with respect to the initial state width L 0 allowed a classification of dynamic regimes independently of the shape of the initial state. In particular, the interacting strength g was scaled asg = gL −s 0 with s ≈ 3/4. This scaling is meaningless for low values of L 0 because in this case, the initial participation number is not of the order of L 0 but is of the order of the maximum AL length 0 (W ) ∼ 96/W 2 . Figure 5(a) shows that using the scaled nonlinearityg and the scaled spectral entropyH = H L s 0 makes the curves H (g) corresponding to L 0 20 collapse to a single curve, except in the self-trapped region. Figure 5(b) shows the Lyapunov exponent as a function ofg; the curves also collapse for L 0 12. The Lyapunov exponent itself is independent of L 0 (no scaling of the variable λ), which is not surprising as it does not measure the absolute distance between two quantum trajectories but the timescale of their exponential divergence.
As shown by the black solid line in figure 5(a), the scaled spectral entropyH is very well fitted (except in the self-trapped regiong > 100) by the log-normal law with three free parameters, the amplitude h 0 , the centre G of the distribution and the 'standard deviation' σ . The study of these fitting parameters as a function of the disorder W provides a full characterization of the dynamic regime. Instead of representing the fit parameters h 0 and G, we prefer to use more physical quantities, namely the maximum value ofH , 6(a)) and the rescaled interaction strengthg c = exp(G − σ 2 ) (figure 6(b)) corresponding to this maximum. The dependence of the standard deviation σ on W is displayed in figure 6(c). The quantityH max is a decreasing function of W : as the localization length decreases with W , so does the overlap between two neighbour Anderson eigenstates. Conversely,g c is an increasing function of W : the number of resonances is maximum when v NL n is comparable to the typical energy between two neighbour states, itself of the order of the bandwidth ∼ 4 + W ;g c is thus independent of W at low disorders and increases with W for larger disorders. Finally, in figure 6(c), one can notice that the log-normal  curve becomes sharper as disorder increases, which can be attributed to the fact that AL is more robust against interactions at high disorders. The fact that the spectral entropy can be determined from a relatively short time interval also allows one to follow the evolution of the dynamics. Figure 7 shows the evolution of the spectral entropy H (t), defined as the spectral entropy calculated in an interval [t, t + 200]; we shall call this, for short, dynamic spectral entropy. For the low value of g = 10 (figure 7(a)), after the destruction of AL, the wavepacket diffusion produces a dilution and a consequent diminution of the nonlinearity that reduces the chaoticity and thus the spectral entropy of the system. For the high value of g = 1000 ( figure 7(b)), one sees a more complex interplay of different regimes: the initial state is initially frozen in a self-trapping regime, but this regime is unstable (for this particular set of parameters): if a fraction of the packet escapes the self-trapped region, the nonlinear contribution v NL n decreases, which leads to a weakening of the trapping. The destruction of the self-trapping regime takes a much longer time than the destruction of the AL. One first observes a small decrease of the spectral entropy, as some eigenstates have been absorbed and no longer interact. Then, when v NL n has decreased sufficiently in the centre of the packet, the system enters the chaotic regime.
Plots (a) and (b) of figure 8 represent the dynamic spectral entropy H (t) for three values of L 0 , corresponding to the same value ofg (0.62 in (a) and 6.17 in (b) and (c)). The curves converge for long times, providing evidence for the existence of a universal shows that additional ad hoc scaling gives an even more universal character to this behaviour. asymptotic regime, independent of the initial state, once the proper scaling ong is applied. Additional scaling can even be used to describe the transition from the self-trapping regime to the chaotic regime, as shown in figure 8(c). Although we have at present no physical explanation for the exponents appearing in these scaling laws, these results strongly support the idea that the features of the nonlinear dynamics should be scaled with respect to the initial state [25].

Conclusion
We have shown that spectral entropy is a pertinent indicator of the chaoticity of the dynamics and allows for a very good characterization of the dynamics regimes. Moreover, it can be calculated dynamically and thus also gives information on the evolution of the dynamics. The Lyapunov exponent gives, in the present context, a less complete characterization of the dynamics. It describes very well the progressive destruction of the AL by the onset of chaotic behaviour, as the spectral entropy does, but it does not give a good description of the self-trapping regime, where 'phase chaos' is still present.
The argument presented in section 4 suggests that the log-normal shape is linked to the exponential localization of the Anderson eigenstates. It is a little surprising to find this link between the behaviour of a strongly nonlinear system and the eigenstates of the corresponding linear system. From a mathematical viewpoint, spectral analysis (in the sense of the determination of eigenvalues and eigenvectors) cannot be directly applied to nonlinear systems and, at present, there is no alternative analytical method for analysing nonlinear systems. The results discussed in this work and in [25,26] indicate that scaling over the initial state is a promising tool for the study of the emerging field of nonlinear quantum mechanics. In this context, it is worth noting that we have at present no theory for the values of the exponents appearing in these scaling laws and that we have considered only a particularly simple family of initial states. Moreover, there is no experimental verification of these scaling laws. It thus appears that a large field of investigation is opened up for both theoreticians and experimentalists in the near future. to finally deduce β, γ and δ. This method allows an accurate calculation of the Lyapunov exponent in the discrete system considered in this work, but it can, in principle, be applied also to the continuous Gross-Pitaevskii equation.