Subdiffusion of nonlinear waves in quasiperiodic potentials

We study the spatio-temporal evolution of wave packets in one-dimensional quasiperiodic lattices which localize linear waves. Nonlinearity (related to two-body interactions) has destructive effect on localization, as recently observed for interacting atomic condensates [Phys. Rev. Lett. 106, 230403 (2011)]. We extend the analysis of the characteristics of the subdiffusive dynamics to large temporal and spatial scales. Our results for the second moment $m_2$ consistently reveal an asymptotic $m_2 \sim t^{1/3}$ and intermediate $m_2 \sim t^{1/2}$ laws. At variance to purely random systems [Europhys. Lett. 91, 30001 (2010)] the fractal gap structure of the linear wave spectrum strongly favors intermediate self-trapping events. Our findings give a new dimension to the theory of wave packet spreading in localizing environments.


I. INTRODUCTION
In one dimension, a wave packet of noninteracting particles subject to a random potential does not diffuse because of Anderson localization, due to the exponential localization of the eigenstates of the underlying Hamiltonian [1,2]. Instead the presence of interactions is expected to act against localization, though the actual mechanism may be highly nontrivial and may depend on the type of disorder and interaction. This is a problem of fundamental importance for many systems in different contexts. From the theoretical side, it has been largely studied by using discrete lattices. The interaction is often included by means of mean-field theories, where it enters in the form of a nonlinear local term in a nonlinear Schrödinger-like equation [3][4][5][6][7][8][9][10][11][12][13].
Numerical simulations of wave packets propagating in a random potential -with the interaction included via a nonlinear mean-field term -showed that the presence of interaction indeed destroys localization and leads to a subdiffusive growth of the second moment of the wave packet in time as t γ [4][5][6][7][8][10][11][12]. In particular it was predicted that at large t, the coefficient γ should converge to 1/3 in a regime of so-called "weak chaos", as opposed to normal diffusion where γ = 1 and the wave packet width grows as √ t. A transient regime of "strong chaos" was also identified, where γ = 1/2 [9,11,12]. The occurrence of these different regimes can be predicted by comparing the nonlinear frequency shift introduced by the expanding wave packet to the typical energy scales of linear spectra of random models.
Exponential localization for noninteracting quantum particles (or linear waves) can be found also in systems which are not truly disordered. An example is provided by quasiperiodic potentials, which are of great interest by themselves [14,15]. These systems can be considered to lie between the two extreme cases of a perfectly periodic system and a pure random potential. By tuning the parameters of a quasiperiodic system, the localization properties can change dramatically -from having all extended states to all localized states. In recent years, exponential localization has been observed with light propagating in quasiperiodic photonic lattices [16], as well as with ultracold atoms propagating in a bichromatic optical lattice [17]. Notably in both cases the inclusion of interaction is experimentally feasible, by using a Kerr medium for light and tuning the scattering length by means of a suitable magnetic field for atoms.
Numerical simulations studying nonlinear dynamics of wave packets have been performed also in the case of quasiperiodic systems [18][19][20][21]. In particular for exponentially localized linear waves, nonlinearity yields subdiffusive spreading of wave packets as well [20]. However, there are clear indications that the coefficient γ, at least at finite spreading times, is significantly larger than the one observed in random systems. Nonlinear effects have been also studied in experiments using ultracold atoms and light propagating in photonic lattices. In both cases it has been shown that nonlinearity acts against localization [16,22].
The purpose of the present work is to clarify the details of the spreading mechanism leading to the destruction of the localization in quasiperiodic systems, and to address differences and similarities between quasiperiodic and purely random potentials. We extend and refine previous numerical investigations by pushing the simulations to much longer times, thus allowing for the identification of the strong and weak chaos regimes in quasiperiodic systems and compare the situation with known properties of purely random systems. For this purpose, we use two different models; namely, a discrete nonlinear Schrödinger equation (DNLS) and a quasiperiodic version of the quar-tic Klein-Gordon (KG) lattice model.
A key result of the present work is that a regime of weak chaos is indeed observed in the long time spreading of nonlinear wave packets propagating in quasiperiodic systems; in particular we find that the asymptotic value of the spreading coefficient γ is 1/3 as in purely random systems, thus showing that this behaviour is rather general and model independent. Another similarity with purely random systems is the occurrence of selftrapping: when the nonlinear interaction is large enough to shift the mode frequencies so strongly that they are tuned out of resonance with all nonexcited neighbouring modes, a part of the wave packet remains spatially localized [3,6,20]. However as opposed to the random system, in the quasiperiodic case partial self-trapping is also possible for weaker nonlinearities. This is due to the complexity of the linear wave spectrum which exhibits a fractal gap structure of sub-bands. Self-trapping gives rise to transient spreading regimes characterized by an intermediate large exponent γ; we call this effect "overshooting". Finally, we have also observed signatures of strong chaos, but detection of this regime is difficult in quasiperiodic systems, since it is often masked by overshooting and partial self-trapping, which occur on the same temporal scales.

II. MODELS
We consider two different models: the first is the one-dimensional (1D) quasiperiodic discrete nonlinear Schrödinger equation (DNLS), defined by the Hamiltonian where j labels the lattice sites and V j = λ cos(2παj + ϕ). The quantity |ψ j | 2 gives the probability to find a particle at site j. The first term in Eq. (1) describes the hopping between nearest neighbouring sites, the second term describes the quasiperiodic on-site energy, while the last term represents the mean-field interaction energy and introduces the nonlinearity. The key parameters of this Hamiltonian are the strength of the quasiperiodic potential λ (for λ = 0 the lattice is periodic with period 1), the strength of the nonlinearity β (for β = 0 the particles are non-interacting), and the irrational number α which causes the underlying potential to be quasiperiodic. In fact, when α is irrational the cosine adds a second periodicity which is incommensurate with respect to the underlying periodicity given by the discreteness of the system. Let us note that, without any loss of generality, one can always choose −0.5 < α ≤ 0.5 since Eq. (1) and the subsequent Eq. (3) are invariant under a shift of α by an integer number. Choosing the value of α in this interval has the additional advantage that the wavenumber associated to V j , k α = 2πα, rests in the Brillouin zone of the underlying discrete lattice. As a convenient choice for this work, we use α = ( √ 5 − 3)/2 [23]. The phase ϕ is a phase shift between the two lattices. The equations of motion associated with Eq. (1) are i∂ψ j /∂t = ∂H/∂ψ * j , or The above set of equations conserve the energy of Eq. (1), as well as the total norm S = j |ψ j | 2 of the initial wave packet (we always assume S = 1). The set can be used for a wide range of applications, including ultracold atoms expanding in optical lattices [17,20,24,25] and light propagating in photonic lattices [16,26]. The second model is a quasiperiodic version of the quartic Klein-Gordon (KG) lattice, given by where u j and p j are the generalized coordinates and momenta on the site j andṼ j = 1 + (1/2) cos(2παj + ϕ). The energy associated with lattice site j is The equations of motion are generated by ∂ 2 u j /∂t 2 = −∂H/∂u j , yielding This set of equations conserve the total energy H = j E j ≥ 0 only. The KG model has also been extensively studied, since it can give a simple description of the nondissipative dynamics of anharmonic optical lattice vibrations in molecular crystals [27].
The total energy of the system H serves as a control parameter of nonlinearity, analogous to β for the DNLS model. In fact, for small amplitudes the equation of the KG chain can be approximately mapped onto a DNLS model [28] using βS ≈ 6λH. Further analytics will be discussed only in terms of the DNLS chain, since it is then straightforward to project to the KG model. For the DNLS model we measure the spreading of the wave packet by tracking the quantity n j ≡ |ψ j | 2 /S, hereafter named norm density consistently with the notation of Refs. [6,7,9,11]. The key quantities that we use to describe the time evolution of the expanding wave packet are the second moment m 2 = j n j (X − j) 2 , (X = j n j j) which quantifies the spatial extent of the wave packet, and the participation number P = 1/ j n 2 j which measures the number of significantly populated sites. A combination of these two quantities ξ = P 2 /m 2 , called compactness index, gives a measure of the sparsity of a wave packet [7]. For the KG we do exactly the same, but replacing the norm density n j with its counterpart E j /H, which is the normalized energy density.

III. RELEVANT ENERGY SCALES
Neglecting the nonlinear term in Eq. (2) reduces to an eigenvalue problem Here the index ν labels the different normal modes A ν,j and eigenvalues E ν . The coefficient 1/(2λ) in Eq. (5) was chosen so that the linear parts of H DNLS and H KG would correspond to the same eigenvalue problem: the linear KG model can then likewise be identically reduced to Eq. (6), under the substitution E ν = 2λ ω 2 ν − 1/λ − 1 , where ω ν are the corresponding eigenfrequencies.
Eq. (6) is also known as the Aubry-Andrè model [15]. The localization properties of this model are well known and extensively studied both analytically and numerically [14,15,20,25,[29][30][31]. A transition occurs from an extended regime to a localized regime at λ = 2. For λ < 2 all normal modes A ν,j are extended over the entire lattice, at λ = 2 they are critical, while for λ > 2 they are exponentially localized in the form A ν,j ∼ e −|j−jν |/ℓ , where j ν is the localization center and ℓ = 1/ ln(λ/2) is the localization length (notice that it is the same for all the modes) [15]. Since we are interested in the interplay between localization and nonlinearity, we will focus exclusively on the regime λ > 2.
In order to quantify the spatial extent of a given eigenstate, it is convenient to define a localization volume is the second moment of |A ν,j | 2 and X ν = j j|A ν,j | 2 is its center of norm [32]. The localization volume is used to estimate the number of modes which interact with a given mode ν. We show its meaning schematically in Fig. 1a. The modes that interact with a given reference mode ν are those whose center of norm lies in an area V ν around it. The average localization volume V is then found by numerically diagonalizing the linear system, calculating V ν for each eigenmode, and then averaging over all eigenmodes. A plot of this quantity as a function of the potential strength λ is shown in Fig. 1b.
The spectrum for λ > 2 is purely dense-point, characterized by the presence of an infinite number of gaps and bands. A plot of the Aubry-Andrè model's spectrum as a function of λ is shown in Fig. 1c. In this figure, one clearly sees the presence of two major gaps dividing the spectrum in three parts, each of them divided in turn in three smaller parts, and so on [33]. We call these portions of spectrum separated by the largest gaps "mini-bands". For our purposes, it is enough to consider a division of the spectrum in M = 3 or at most in M = 9 mini-bands. Smaller mini-bands have vanishingly small effects on the time evolution of wave packets.
Let us introduce two energy scales associated with the linear system [6,32]. The first one, ∆, is the full width of the spectrum, defined as the difference between the largest and the smallest eigenvalues: ∆ = max{E ν } − min{E ν }. The second one, d, is the mean spacing of eigenvalues within a single mini-band and within the range of a localization volume. Let us explain how we calculate this quantity. We consider a given mini-band and all the eigenstates that lie in it. For each eigenstate ν, we calculate its localization volume V ν and then we form the subset of the other eigenstates, {µ}, belonging to the same mini-band and interacting with it, namely, those fulfilling the condition |X ν − X µ | < V ν /2. The average number of states in the subset can be estimated as V /M . Then we calculate the energy spacings within this subset. This procedure is repeated for each eigenstate in the band and the average gives the mean spacing d.
The number of mini-bands M to be used in the calculations of d depends on λ. For a given λ we choose M in such a way that the localization volume V satisfies the condition V /M > 2. This implies that, on average, there are at least two eigenstates within the subset {µ} that we can use to calculate the average energy spacings. We always consider λ > 2.1; therefore it is enough to divide the spectrum at most in nine mini-bands. As λ is increased the average localization volume of the eigenstates V decreases -therefore at some point we have to consider the spectral separation into smaller mini-bands. In practice we consider M = 9 mini-bands for 2.1 λ 2.2, M = 3 mini-bands for 2.2 λ 2.75 and just one band (i.e., the full spectrum) for λ 2.75. A plot of the energy scales ∆ and d as a function of λ is shown in Fig. 2. The dashed vertical lines represent the values of λ where the number of mini-bands changes in the calculation of d.
Similarly to the case of disordered systems [6,7,9], the scales ∆ and d of the linear spectrum (which are frequencies in the present setting of nonlinear wave equations) must be compared to the frequency shift caused by the nonlinearity. Indeed a single oscillator which satisfies the equation of motion iψ = V ψ + β|ψ| 2 ψ experiences a nonlinear frequency shift δ = β|ψ| 2 away from its linear frequency V . For many oscillators, we can conveniently use the eigenstates of the linear Aubry-Andrè model as a decomposition basis of the wave function ψ j : (2) can then be rewritten for the evolution of the normal mode amplitudes: where I ν,ν1,ν2,ν3 is an overlap integral involving four normal modes: As discussed in Ref. [11], one can introduce a norm density also in the normal mode space, n ν = |φ ν | 2 ; as the packet spreads and after averaging over many realization, this quantity becomes almost identical to the norm density n j = |ψ j | 2 and the frequency shift can be expressed as δ ∼ βn, where n is a characteristic norm density. In the KG model δ is proportional to the energy density E and, within our formalism, can be obtained by the small amplitude mapping. When δ < d, the mode frequencies in a wave packet are only weakly shifted, and a small fraction of these modes will resonantly and strongly interact with each other. Following the terminology of Refs. [6,7,9], we say that this is a regime of weak chaos. Conversely, when d < δ < ∆, the mode frequencies in a packet are strongly shifted and almost all of them will resonantly and strongly interact with each other. This is labeled a regime of strong chaos. When finally δ > ∆, the mode frequencies are shifted so strongly that they are tuned out of resonance with all nonexcited neighbouring modes. An excited mode in this condition may stay localized, i.e., self-trapped, for long or even infinite times. The meaning of these regimes will be further clarified in the next section.

IV. EXPECTED SPREADING REGIMES
As one can see in Eq. (7), the presence of nonlinearity in the DNLS model introduces a coupling between eigenstates of the underlying linear spectrum. It has already been observed numerically and experimentally that this Comparing the nonlinear frequency shift δ with the energy scales ∆ and d one can predict the different spreading regimes of weak chaos (δ < d), strong chaos (d < δ < ∆) and self-trapping (δ > ∆). The separation between the three regimes should not be interpreted as a sharp boundary, but as a smooth crossover.
leads to a subdiffusive spreading of wave packets, i.e. its second moment grows asymptotically as m 2 ∼ t γ with γ < 1 [20,22]. However, a systematic investigation of the behaviour of the exponent γ in different regimes of strong and weak chaos, and self-trapping, have not been done so far. In this section, we approach this issue by first comparing the nonlinear frequency shift δ = βn with the energy scales ∆ and d, in such a way as to introduce the different spreading regimes expected to be observed in the subsequent numerical simulations.
Let us consider an initial wave packet with norm density n and localization volume L larger than the average localization volume of the eigenstates of the linear spectrum, L ≥ V . If δ > ∆, nonlinearity is so strong that all the participating normal modes within the wave packet are shifted out of resonance with respect to the non-excited neighbourhood; therefore spreading is largely suppressed and a significant part of the wave packet remains self-trapped [34]. If instead δ < ∆, we are no longer in the self-trapping regime and can distinguish two sub-cases: on one hand, when δ > d, strong chaos is realized; that is, all the modes in the packet are resonantly interacting with each other, thus producing an efficient spreading. On the other hand, when δ < d, weak chaos is obtained: only a fraction of modes interact resonantlythe localization is still destroyed, but spreading is slower.
If L < V the estimate of the self-trapping transition is done as before, that is by comparing δ = βn with the spectrum width ∆. If self-trapping is avoided, however, the wave packet initially spreads also in absence of nonlinearity, eventually filling the localization volume V . Consequently the initial norm density n is reduced toñ ≈ nL/V , due to linear time evolution -the relevant nonlinear frequency shift must now be calculated by using this reduced densityñ. Apart from this detail, which originates from the initial dynamics at short times, the asymptotic spreading regimes are the same as before.
Studies performed on random systems have shown that the basic mechanism that destroys localization is the presence of resonances in mode-mode interaction [6,7,9]. This leads to chaotic dynamics within a part of the wave packet, and to a subsequent subdiffusive spreading. Here we apply the same theory to the case of quasiperiodic systems. For the reader interested in the formal details of the theory, we refer to the previous articles [7,9].
Let us estimate the number of resonant modes in the packet, which is a key quantity that determines the type of spreading behaviour. According to Eq. (7), due to nonlinearity, the evolution of a given normal mode is affected by any three (triplet) modes. The coupling is the largest if the triplet modes have large amplitudes and if the overlap integrals are large, i.e., if the triplet modes are close enough in space to the given normal mode. Some of these triplet modes may affect the dynamics of the chosen mode ν strongly, some weakly. To distinguish these triplet groups, we follow [9] and apply perturbation theory to first order in β. It follows that the amplitude of a normal mode ν inside the wave packet is changed by a given triplet of other wave packet modes µ = {µ 1 , µ 2 , µ 3 } as where From now on, we assume that all the modes that belong to the packet (i.e., that are located between the two exponential tails of the wave packet) have the same norm density equal to n. The perturbation approach breaks down and resonance sets in when √ n < |φ (1) ν |. Substituting Eq. (9) for φ (1) ν one can rewrite the last inequality as This expression tells us that the resonance condition, for a given normal mode ν, is fulfilled if there is at least one triplet of modes µ that satisfies inequality (11). The probability for the onset of a resonance can therefore be calculated with the following statistical numerical analysis [6,32]. For a given normal mode ν, we define R ν, µ0 = min µ {R ν, µ }. Collecting R ν, µ0 for many modes and many values of the phase ϕ, we find the probability density distribution W(R ν, µ0 ). From this quantity we can calculate the probability P for a mode, with norm density n, to be resonant with at least one triplet of other modes at a given value of the interaction parameter β. This is obtained by integrating W(R ν, µ0 ) from zero to βn An example of probability density W(R ν, µ0 ) for λ = 2.5 is shown in Fig. 3 (red line). For comparison we also show the same quantity for the random DNLS model (black line), as discussed in [6,32]. Except for fine structures, like small sharp peaks appearing in the quasiperiodic case, the overall behaviour is qualitatively very similar in the two cases. In particular, the probability density W tends to a finite constant value C when R ν, µ0 → 0. As a consequence, for small values of βn, a non-zero fraction of modes in the packet is resonant. The probability to be resonant is given by P ∼ Cβn, thus we are in the weak chaos regime. For large values of βn, instead all the modes interact resonantly and P = 1; we are then in the strong chaos regime. Following the reasoning presented in [9], this implies that also in the quasiperiodic case, as in disordered systems, we may expect to find m 2 ∼ t 1/3 in the weak chaos regime and m 2 ∼ t 1/2 in the strong chaos regime. Note the strong chaos regime can only exist as a transient regime: as the wave packet spreads, its norm density n decreases, and eventually will reach a situation where βn < d. At this point, a crossover from strong to weak chaos is expected to occur during the time evolution [11].
Let us finally stress that the "transition lines" that we have introduced by comparing the nonlinear frequency shift with the typical energy scales of the linear spectrum do not define sharp phase transitions between different spreading regimes. Instead, we may expect to see a relatively smooth crossover, such that the regimes of selftrapping, strong chaos and weak chaos should be clearly identified only far from the transition lines.

V. TIME EVOLUTION
We perform extensive numerical simulations solving Eqs. (2) and (5) for different sets of parameters {λ, β} and {λ, E}, respectively. For each choice of parameters we average over N different realizations of the quasiperiodic potential obtained by randomly changing the phase shift ϕ. For initial conditions, we use compact wave packets that lie in the center of our computational box, taking care that during the time evolution the wave packet never reaches the box boundaries. The number of realizations considered varies between 100 and 500 and the number of lattice sites between 200 and 2000. To solve the equations of motion, we use symplectic integration schemes of the SABA family [7,35] that allow us to reach large integration times with good accuracy [36].
In order to quantify the type of subdiffusive behaviour, we calculate the exponent γ by considering the logarithm of the second moment log 10 m 2 for different realizations of the potential. We compute the average value log 10 m 2 and its statistical error, given by the standard deviation divided by the square root of the number of realizations N . Then the value of γ at a given time t is calculated by applying a linear fitting procedure to the curve log 10 m 2 within a fixed time interval around log 10 t. By repeating this procedure at different t, we extract the behaviour of γ as a function of time and its relative statistical error.

A. Results of the DNLS model
Let us first show our results for the DNLS model. For the initial wave packet, we choose a square shaped distribution which equally populates L lattice sites with norm density n j = n = 1/L. In Fig. 4 we present a representative set of simulations for λ = 2.5. We choose L = 13, which gives an initial localization volume larger than V .
The different panels show the time evolution of the second moment log 10 m 2 , the spreading exponent γ, the participation ratio log 10 P , and the compactness index ξ . The width of the curves for log 10 m 2 , log 10 P and γ corresponds to the statistical error. The values of the nonlinear frequency shift δ induced by the initial wave packets used in these simulations are shown in Fig. 2 (empty downward triangles) in order to compare them to the relevant energy scales ∆ and d.
In all simulations we observe that nonlinearity causes the wave packet to spread. The spreading starts earlier when β is larger. We find that the spreading is always subdiffusive (γ < 1), confirming the result of previous works [20,22]. Subdiffusion is seen both in the second moment m 2 and in the participation ratio P , except for the largest value of β (yellow curves in Fig. 4). In the latter case, P saturates to a constant value after a transient time -a clear signature of self-trapping. This observation of self-trapping only for β = 100 is consistent with the energy scale arguments schematically represented in Fig. 2. In the absence of self-trapping, the compactness index ξ saturates to a constant value, indicating that the wave packet spreads but does not become more sparse. Conversely, in the presence of self-trapping the central part of the wave packet remains spatially trapped while its tails keep expanding, thus resulting in a wave packet that becomes more sparse during the evolution -nicely quantified by the compactness index which decreases to zero. We notice that the portion of packet that is expanding is characterized by a value of γ larger than 1/2. After an initial increase, γ reaches a maximum and then decreases to smaller values. In this regime, the evolution is rather complex. A similar behaviour was previously obtained also in random systems [11,12]. The transient large values of γ may be due to a nontrivial interacting mechanism that takes place between the expanding part and the self-trapped portion, resulting in faster spreading -an effect labeled as "overshooting".
For the lowest values of β the energy scale arguments suggest the occurrence of weak chaos. Indeed for β = 0.5 and 1 the exponent saturates asymptotically around the theoretical value γ = 1/3 (red and green curves in Fig. 4), as expected. It is worth mentioning that this asymptotic exponent is the same as in random systems [9,11]; meaning that the mechanism leading to destruction of exponential localization is rather universal.
In difference to the random case, here during the time evolution, the value of γ temporarily increases above 1/3, eventually reaching its asymptote only at longer times. This is an overshooting similar to the one that we have discussed above for the self-trapping regime, but occurring also for weaker nonlinearity. This effect is unique to the quasiperiodic system and is likely due to the presence of an infinite number of mini-bands and gaps in the linear spectrum of the Hamiltonian, which causes a temporary self-trapping of portions of the expanding wave packet in one or more energy gaps between mini-bands. This partial self-trapping is different from the self-trapping that occurs when δ > ∆, where all the packet modes are simultaneously shifted out of resonance. For this reason partial self-trapping is not detectable as a saturation of the participation number P and can only be seen indirectly as an overshooting in the exponent γ. The two simulations for β = 5 and 10 lie in a range of energy were we expect to see strong chaos (blue and magenta curves in Fig. 4). As already said in the previous section, the strong chaos regime is transient: one should find a value of γ around 1/2 for a few decades of time, eventually decreasing towards the asymptotic value 1/3. The two corresponding curves in Fig. 4 indeed exhibit a behaviour which qualitatively agrees with this expectation. The value of γ first rises up to 1/2, oscillates around this value and then starts to decrease as predicted. However, especially for large β, we also observe values of γ larger than 1/2. As in the weak chaos regime, this overshooting again is evidence of partial self-trapping. Its mechanism is also transient and occurs in the same time intervals where strong chaos is expected. For this reason, while weak chaos is clearly observed in our simulations, strong chaos and partial self-trapping tend to overlap, thus producing a more complex evolution of the wave packet in quasiperiodic systems than in random systems.
In Fig. 5 we show the results of simulations for λ = 2. the predicted behaviour is either strong chaos or a regime in between strong and weak chaos. What we observe numerically is a growth of the spreading exponent γ up to 1/2 and even to larger values, followed by a decrease towards 1/3. In most cases, our simulations show a significant overshooting due to partial self-trapping. It is worth mentioning that this effect is larger for weaker disorder strength λ, consistent with the fact the linear spectrum exhibits larger mini-gaps in this regime (see Fig. 1). Finally for {λ, β} = {3.5, 50}, we observe self-trapping, as expected.
In conclusion, from the analysis of the results of the DNLS model for different values of λ we find that the energy scale arguments and the model discussed in Section IV correctly explain the overall trend of the numerical simulations and the separation between different spreading regimes in the parameter space.

B. Results of the KG model
Due to the existence of a mapping between KG and DNLS, we expect to observe the same spreading regimes in the two models. This has been already proven in purely random systems where the two models reveal similar qualitative results in a wide range of parameters [6,7,9,11,12]. Despite this similarity, the study of the KG model remains interesting for at least two reasons. On one hand, it allows for testing the generality of the result in a case where there is just one conserved quantity. This is highly nontrivial -especially for selftrapping -for which rigorous results have been recently derived only in the case of Hamiltonians conserving both energy and norm [3]. On the other hand, the KG model is advantageous from a numerical point of view. The fact that there is just one conserved quantity results in two orders of magnitude faster integration speed within the same integration error.
Similarly to what was done for the DNLS model, we initially set the compact wave packets to span a width L = 13 (unless otherwise stated) centered in the lattice, such that each site has equal energy E j = E = H/L. This is implemented by setting initial momenta of p = ± √ 2E with randomly assigned signs and zero coordinates. The values of initial energy densities E are chosen to give expected spreading regimes of asymptotic weak chaos, intermediate strong chaos, and dynamical crossover from strong chaos to the slower weak chaos subdiffusive spreading [9].
The results of the time simulations are shown in Fig. 6, while the expected spreading regimes are given in Fig. 2 (full upward triangles) [37]. As one can see by comparing Fig. 6 with Fig. 4, the qualitative behaviour of the two models is rather similar. After initial transients, which increase with decreasing nonlinearity, all KG simulations reveal subdiffusive growth of the second moment m 2 according to power law m 2 ∼ t γ with γ < 1. If selftrapping is avoided, all simulations show a similar subdiffusive behaviour for the participation numbers; moreover, the wave packets remain compact as they spread, since compactness indices at the largest computational times saturate around a constant ξ ≈ 3.5 ± 0.25. For the two smallest values of initial energy density E = 0.05 and E = 0.01, the characteristics of the weak chaos regime are observed, namely, the exponent γ saturates around 1/3 (red and green curves in Fig. 6) after a transient time. We stress that the only difference from the purely random systems is the overshooting phenomenon at transient times. This effect is an inherent property of quasiperiodic systems which inevitably manifests itself in all spreading regimes, while in the disordered case it was shown to occur only in the regime of self-trapping [11,12].
For the two energy densities E = 0.055 and 0.075 we suggest strong chaos, with characteristics similar to the DNLS case. The simulation with E = 0.055 (blue curves in Fig. 6) indeed exhibits the typical behaviour of the strong chaos scenario: the characteristic exponent γ increases up to predicted value 1/2 and remains so for about two time decades, followed by a crossover with γ decreasing to the weak chaos dynamics. There is also another possibility for larger E = 0.075, when intermediate strong chaos is masked due to partial self-trapping (magenta curves in Fig. 6). Thus, γ shows values larger then 1/2 but still with subsequent decay to slower subdiffusion. Here, we would like to strongly emphasize that none of the simulations exhibit pronounced deviations from strong or weak chaos regimes of spreading, i.e. longlasting overshooting with γ > 1/2, or significant slowing down to values γ < 1/3.
Finally, for E = 1.0 the dynamics enter the selftrapping regime, as our theory predicts. There a major part of the initial excitation stays localized, while the remainder spreads (yellow curves in Fig. 6). The participation numbers, therefore, do not grow significantly and log 10 P starts to level off at large time (Fig. 6, left panel, bottom, yellow curve). In contrast, the small spreading portion yields a continuous increase of the second moment m 2 (Fig. 6, left panel, top, yellow curve), which initially is characterized by large values of γ > 1/2 (howbeit, for larger time γ decreases). Consequently, the compactness index ξ (Fig. 6, right panel, bottom, yellow curve) drops down to small values indicating deep self-trapping regime. Note that a similar behaviour has been observed before in purely random systems [11,12].
Unusually large values of m 2 can be explained by local trapping-detrapping processes in the evolving wave packet. The corresponding dynamics is in strong nonequilibrium -its theoretical description has yet to be developed.
C. Role of the shape of the initial wave packet In this subsection we show that the results discussed so far do not depend on the shape of the initial wave packet. Besides its theoretical interest, this issue is also relevant from the point of view of experiments, where it is not always possible to design the wave packets at will.
In the previous sections, we have used a square distribution as the initial wavepacket. Now, inspired by the experiments with ultracold atoms, we consider initial wave packets with the shape of a Gaussian distribution or a Thomas-Fermi (TF) distribution. The Gaussian wave packet centered around the site j = 0 has the form where σ is a parameter controlling the width of the packet while C 1 is a constant factor that can be determined by using the normalization condition j |ψ j | 2 = 1. A Thomas-Fermi wave packet is instead defined by in the region where |j| < R and ψ j = 0 otherwise. The parameter R is the Thomas-Fermi radius characterizing the width of the distribution, while the constant C 2 is a normalization factor. These two distributions are of interest when considering ultracold bosons initially released from an harmonic trap in the Gross-Pitaevskii regime [38].
In Fig. 7 we show the time evolution in the DNLS model of the second moment of the expanding wave packet, log 10 m 2 (top row) and of the spreading exponent, γ (bottom row), using initially a Gaussian (left column) and a TF (right column) wave packet distribution. In the insets we also show the compactness index ξ , in order to identify the self-trapping regime. We choose the width of the initial distributions (σ and R) so that the nonlinear frequency shift is similar to the one already used for the simulations in Fig. 4 [39]. In particular we use σ = 5 and R = 7.5, yielding a nonlinear frequency shift δ ≈ β/13. The values of β used in Fig. 7 are the same as those previously considered. From the comparison between the results of Fig. 7 and Fig. 4, we can conclude that the shape of the initial wave packet does not affect the overall behaviour of the time evolution, nor its interpretation in terms of regimes of weak and strong chaos, self-trapping, and overshooting. This suggests the results that we have obtained are rather general and that the nonlinear frequency shift δ is the only key parameter controlling the dynamics of the wave packet.

D. Application to cold atoms
The DNLS model can be used to simulate the dynamics of bosons in optical lattices at zero temperature [24] and in the tight-binding regime, where the DNLS equation corresponds to a discretized version of the Gross-Pitaevskii (GP) equation for the dynamics of a Bose-Einstein condensate in the single-band approximation. The validity of this mean-field theory is not ensured for those dynamical regimes where GP predicts chaos [40], which can be viewed as a signature of a large depletion of the condensate. For this reason, in the presence of disorder the theory fails to predict the long time evolution of observables directly related to small scale fluctuations and long-range coherence. However, for coarse-grained observables, like the width of the wave packet in real and momentum space, or the participation number, the predictions of the theory remain very good even in regimes where the depletion is expected to be large, long after the random fluctuations prevent the prediction of fine scale structures. This has been recently shown in Ref. [41] by comparing the predictions of the GP equation with one beyond mean-field theory in numerical simulations within timescales of the order of typical experiments with cold atoms and long enough to observe the effects of depletion and chaotic dynamics. Indeed our analysis is essentially based on coarse-grained observables. In addition, for each set of parameters we also average over many realizations and this extends the validity of the present approach even for longer times, as any residual dependence on small scale fluctuations is further suppressed by the averaging procedure.
When applied to bosons expanding in bichromatic optical lattices, our results provide a consistent interpretation of the experimental data of Ref. [22], where a Bose-Einstein condensate, initially confined in an harmonic trap, is let free to expand in a bichromatic potential. In our dimensionless units, the expansion lasts for times of the order of 10 4 (see [20] for more details) and the width of the atomic cloud increases up to 50 − 100 lattice sites. In this experiment a subdiffusive spreading is observed with exponents γ significantly larger than 1/3 already for weak nonlinearities and even larger than 1/2 for larger nonlinearities [42]. Our work suggests that such large values of γ can be explained in terms of a transient overshooting caused by partial self-trapping in mini-bands.

VI. SUMMARY AND CONCLUSIONS
In this work we have considered the problem of the interplay between localization and interaction in onedimensional quasiperiodic systems. We have investigated the expansion of initially localized wave packets in two different quasiperiodic models, DNLS and KG. We have confirmed that interaction destroys localization, giving rise to a subdiffusive growth of the second moment of the wave packet (m 2 ∼ t γ with γ < 1).
We have interpreted the spreading process in terms of resonances in the mode-mode coupling. In particular, we have identified the different spreading regimes of self-trapping, strong chaos and weak chaos by comparing the frequency shift induced by the nonlinearity with the energy scales extracted from the spectrum of the underlying linear system. For weak and strong chaos regimes we have also predicted the expected spreading exponents γ = 1/3 and γ = 1/2 respectively.
We have performed numerical simulations, which last for much longer times than existing simulations, and we have averaged our results over many realizations. This gave us the possibility to accurately calculate the spread-ing exponent γ and observe the weak chaos regime. A key difference with respect to random systems [6,9] is the occurrence of transient overshooting regimes that we interpret as due to the peculiar structure of the linear spectrum of the quasiperiodic system, which is separated into mini-bands. These mini-bands are responsible for mechanisms of partial self-trapping. Signatures of strong chaos have also been observed, but the temporal overlap of strong chaos and partial self-trapping makes the analysis of the spreading more complex than for random systems. We have also verified that our main results do not depend on the details of the shape of the initial wave packet. This suggests that the nonlinear frequency shift, δ, is the only parameter that controls the dynamics.
Finally, our results provide a consistent interpretation of the subdiffusive spreading observed in experiments with ultracold atoms propagating in bichromatic optical lattices [22].