Subdiffusion of nonlinear waves in quasiperiodic potentials

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, although 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 mostly 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–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–8, 10–12]. In particular, it was predicted that at large t, the coefficient γ should converge to 1/3 in a regime of weak chaos, as opposed to normal diffusion where γ = 1 and the wave packet width grows as $\sqrt {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 also be found in systems that 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 also been performed in the case of quasiperiodic systems [18–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 also been 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 this work is to clarify the details of the spreading mechanism leading to the destruction of 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 quartic Klein–Gordon (KG) lattice model.

A key result of this 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 self-trapping: 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 that 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.

We consider two different models: the first is the one-dimensional (1D) quasiperiodic DNLS, defined by the Hamiltonian

$$\begin{equation} H_{\rm DNLS}=\sum_j \left[-\psi_{j+1}\psi_j^*-\psi_{j-1}^*\psi_j +V_j|\psi_j|^2 + \frac{\beta}{2}|\psi_j|^4 \right], \end{equation} \tag{ 1 }$$

where j labels the lattice sites and V_j = λ cos(2παj + φ). The quantity |ψ_j|² gives the probability to find a particle at site j. The first two terms in equation (1) describe the hopping between nearest-neighbouring sites and the quasiperiodic on-site energy, respectively, while the third 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 not 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 equation (1) and the subsequent equation (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 with V_j, k_α = 2πα, rests in the Brillouin zone of the underlying discrete lattice. As a convenient choice for this work, we use $\alpha =(\sqrt {5}-3)/{2}$.⁷ The phase φ is a phase shift between the two lattices. The equations of motion associated with equation (1) are i∂ψ_j/∂t = ∂H/∂ψ*_j, or

$$\begin{equation} {\rm i}\frac{\partial \psi_j}{\partial t} = -(\psi_{j+1}+\psi_{j-1})+V_j\psi_{j}+\beta|\psi_j|^2\psi_j. \end{equation} \tag{ 2 }$$

The above set of equations conserve the energy of equation (1), as well as the total norm $S=\sum_j |\psi_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, 23, 24] and light propagating in photonic lattices [16, 25].

The second model is a quasiperiodic version of the quartic KG lattice, given by

$$\begin{equation} H_{\rm KG}=\frac{1}{2}\sum_j \left[ p_j^2+\tilde{V_j}u_j^2+\frac{1}{2}u_j^4+\frac{1}{2\lambda}(u_{j+1}-u_j)^2\right], \end{equation} \tag{ 3 }$$

where u_j and p_j are the generalized coordinates and momenta on the site j and $\tilde {V_j}=1+(1/2)\cos (2\pi \alpha j+\varphi )$. The energy associated with lattice site j is

$$\begin{equation} \mathcal{E}_j = \frac{p_j^2}{2}+\frac{\tilde{V_j}u_j^2}{2}+\frac{u_j^4}{4}+\frac{(u_{j+1}-u_j)^2}{8 \lambda} +\frac{(u_{j-1}-u_j)^2}{8 \lambda}. \end{equation} \tag{ 4 }$$

The equations of motion are generated by ∂²u_j/∂t² = −∂H/∂u_j, yielding

$$\begin{equation} \frac{\partial^2 u_j}{\partial t^2}=-\tilde{V_j}u_j-u_j^3+\frac{1}{2\lambda}(u_{j+1}+u_{j-1}-2 u_j). \end{equation} \tag{ 5 }$$

This set of equations conserves the total energy $\mathcal {H}=\sum_j \mathcal {E}_j \geqslant 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 [26].

The total energy of the system $\mathcal {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 [27] using $\beta S \approx 6 \lambda \mathcal {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|²/S, hereafter named norm density consistently with the notation of [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=\sum_j n_j(X-j)^2$ ($X=\sum_j n_j j$), which quantifies the spatial extent of the wave packet, and the participation number $P=1/\sum_j n_j^2$, which measures the number of significantly populated sites. A combination of these two quantities ξ = P²/m₂, called the 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 $\mathcal {E}_j/\mathcal {H}$, which is the normalized energy density.

For β = 0 equation (2) reduces to an eigenvalue problem

$$\begin{equation} -A_{\nu,j+1}-A_{\nu,j-1}+\lambda\cos(2\pi\alpha j+\varphi)A_{\nu,j}=E_\nu A_{\nu,j}. \end{equation} \tag{ 6 }$$

Here, the index ν labels the different normal modes A_ν,j and eigenvalues E_ν. The coefficient 1/(2λ) in equation (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 equation (6), under the substitution E_ν = 2λ(ω²_ν − 1/λ − 1), where ω_ν are the corresponding eigenfrequencies.

Equation (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, 24, 28–30]. 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, whereas for λ > 2 they are exponentially localized in the form A_ν,j ∼ e^{−|j−j_ν|/ℓ}, where j_ν is the localization centre and ℓ = 1/ln(λ/2) is the localization length (note 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 $V_\nu =1+\sqrt {12 m_{2}^{(\nu )}}$, where $m_{2}^{(\nu )}=\sum_j (X_\nu -j)^2|A_{\nu ,j}|^2$ is the second moment of |A_ν,j|² and $X_\nu =\sum_j j |A_{\nu ,j}|^2$ is its centre of norm [31]. The localization volume is used to estimate the number of modes which interact with a given mode ν. We show its meaning schematically in figure 1(a). The modes that interact with a given reference mode ν are those whose centre 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 figure 1(b).

Zoom In Zoom Out Reset image size

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 figure 1(c). In this figure, one clearly sees the presence of two major gaps dividing the spectrum into three parts, each of them divided into turn in three smaller parts, and so on⁸. We call these portions of the spectrum separated by the largest gaps 'mini-bands'. For our purposes, it is enough to consider a division of the spectrum into M = 3 or at most into 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, 31]. The first one, Δ, is the full width of the spectrum, defined as the difference between the largest and the smallest eigenvalues: $\Delta =\max \{E_{\nu }\}-\min \{E_{\nu }\}$. 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 into at most nine mini-bands. When λ 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 figure 2. The dashed vertical lines represent the values of λ where the number of mini-bands changes in the calculation of d.

Zoom In Zoom Out Reset image size

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 that satisfies the equation of motion ${{\rm i}}\dot {\psi }=V \psi + \beta |\psi |^2\psi$ experiences a nonlinear frequency shift δ = β|ψ|² 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: $\psi_j=\sum_\nu \phi_\nu A_{\nu ,j}$. Equation (2) can then be rewritten for the evolution of the normal mode amplitudes:

$$\begin{equation} {\rm i} \frac{\partial \phi_\nu}{\partial t}= E_\nu \phi_\nu +\beta \sum_{\nu_1,\nu_2, \nu_3} I_{\nu,\nu_1,\nu_2,\nu_3}\phi_{\nu_1}^*\phi_{\nu_2}\phi_{\nu_3}, \end{equation}\noindent \tag{ 7 }$$

where I_{ν,ν1,ν2,ν3} is an overlap integral involving four normal modes:

$$\begin{equation} I_{\nu,\nu_1,\nu_2,\nu_3}=\sum_j A_{\nu,j}A_{\nu_1,j}A_{\nu_2,j}A_{\nu_3,j}. \end{equation}\noindent \tag{ 8 }$$

As discussed in [11], one can introduce a norm density also in the normal mode space, n_ν = |ϕ_ν|²; as the packet spreads and after averaging over many realizations, this quantity becomes almost identical to the norm density n_j = |ψ_j|² 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 $\mathcal {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 [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 labelled 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.

As one can see in equation (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 leads to a subdiffusive spreading of wave packets, i.e. its second moment grows asymptotically as m₂ ∼ 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, has 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⁹. If instead δ < Δ, we are no longer in the self-trapping regime and can distinguish two sub-cases: on the 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 resonantly—localization is still destroyed, but spreading is slower.

If L < V , estimation of the self-trapping transition is performed as before, that is, by comparing δ = βn with the spectrum width Δ. If self-trapping is avoided, however, the wave packet initially spreads also in the absence of nonlinearity, eventually filling the localization volume V . Consequently, the initial norm density n is reduced to $\tilde {n}\approx nL/V$, due to linear time evolution—the relevant nonlinear frequency shift must now be calculated by using this reduced density $\tilde {n}$. Apart from this detail, which originates from the initial dynamics at short times, the asymptotic spreading regimes are the same as before.

Studies carried out 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 formal details of the theory, see the previous papers [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 equation (7), due to nonlinearity, the evolution of a given normal mode is affected by any three (triplet) modes. The coupling is 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 $\vec {\mu }=\{\mu_1,\mu_2,\mu_3\}$ as

$$\begin{equation} |\phi_\nu^{(1)}|=\beta\frac{\sqrt{n_{\mu_1} n_{\mu_2} n_{\mu_3}}}{R_{\nu,\vec{\mu}}}, \end{equation} \tag{ 9 }$$

where

$$\begin{equation} R_{\nu,\vec{\mu}}\sim \left|\frac{E_\nu+E_{\mu_1}-E_{\mu_2}-E_{\mu_3}}{I_{\nu,\mu_1,\mu_2,\mu_3}}\right|. \end{equation} \tag{ 10 }$$

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 $\sqrt {n}<|\phi_\nu ^{(1)}|$. Substituting equation (9) for ϕ⁽¹⁾_ν, one can rewrite the last inequality as

$$\begin{equation} R_{\nu,\vec{\mu}} < \beta n. \end{equation} \tag{ 11 }$$

This expression tells us that the resonance condition, for a given normal mode ν, is fulfilled if there is at least one triplet of modes $\vec {\mu }$ that satisfies inequality (11).

The probability for the onset of a resonance can therefore be calculated with the following statistical numerical analysis [6, 31]. For a given normal mode ν, we define $R_{\nu ,\vec {\mu_0}}= \min_{\vec {\mu }} \{R_{\nu ,\vec {\mu }}\}$. Collecting $R_{\nu ,\vec {\mu_0}}$ for many modes and many values of the phase φ, we find the probability density distribution $\mathcal {W}(R_{\nu ,\vec {\mu_0}})$. From this quantity, we can calculate the probability $\mathcal {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 $\mathcal {W}(R_{\nu ,\vec {\mu_0}})$ from zero to βn

$$\begin{equation} \mathcal{P}=\int_0^{\beta n} \mathcal{W}(R)\, {\rm d}R. \end{equation} \tag{ 12 }$$

An example of probability density $\mathcal {W}(R_{\nu ,\vec {\mu_0}})$ for λ = 2.5 is shown in figure 3 (red line). For comparison we also show the same quantity for the random DNLS model (black line), as discussed in [6, 31]. Except for fine structures, such as small sharp peaks appearing in the quasiperiodic case, the overall behaviour is qualitatively very similar in the two cases. In particular, the probability density $\mathcal {W}$ tends to a finite constant value C when $R_{\nu ,\vec {\mu_0}}\rightarrow 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 $\mathcal {P}\sim C\beta n$; thus we are in the weak chaos regime. For large values of βn, instead all the modes interact resonantly and $\mathcal {P}=1$; we are then in the strong chaos regime.

Zoom In Zoom Out Reset image size

Following the reasoning presented in [9], this implies that also in the quasiperiodic case, as in disordered systems, we may expect to find m₂ ∼ t^1/3 in the weak chaos regime and m₂ ∼ t^1/2 in the strong chaos regime. Note that 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 self-trapping, strong chaos and weak chaos should be clearly identified only far from the transition lines.

We perform extensive numerical simulations solving equations (2) and (5) for different sets of parameters {λ,β} and $\{\lambda ,\mathcal {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 centre 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, 32] that allow us to reach large integration times with good accuracy¹⁰.

In order to quantify the type of subdiffusive behaviour, we calculate the exponent γ by considering the logarithm of the second moment log₁₀ m₂ for different realizations of the potential. We compute the average value 〈log₁₀ m₂〉 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₁₀ m₂〉 within a fixed time interval around log₁₀ t. By repeating this procedure at different t, we extract the behaviour of γ as a function of time and its relative statistical error.

5.1. Results of the discrete nonlinear Schrödinger equation model

Let us first show our results for the DNLS model. For the initial wave packet, we choose a square-shaped distribution that equally populates L lattice sites with norm density n_j = n = 1/L. In figure 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₁₀ m₂〉, the spreading exponent γ, the participation ratio 〈log₁₀ P〉 and the compactness index 〈ξ〉. The width of the curves for 〈log₁₀ m₂〉, 〈log₁₀ 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 figure 2 (empty downward triangles) in order to compare them to the relevant energy scales Δ and d.

Zoom In Zoom Out Reset image size

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₂ and in the participation ratio P, except for the largest value of β (yellow curves in figure 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 figure 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 note that the portion of the 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. 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 labelled 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 figure 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 contrast 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 probably 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 where we expect to see strong chaos (blue and magenta curves in figure 4). As already mentioned 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 figure 4 indeed exhibit behaviour which qualitatively agrees with this expectation. The value of γ first rises 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 figure 5, we show the results of simulations for λ = 2.2 and λ = 3.5; the corresponding values of nonlinear frequency shift are reported as triangles in figure 2. The values of L are L = 31 for λ = 2.2 and L = 5 for λ = 3.5, both larger than V . For {λ,β} = {2.2,0.18} and {λ,β} = {3.5,5.5} energy scale arguments predict weak chaos. We indeed find a spreading exponent which approaches asymptotically the value 1/3. For {λ,β} = {2.2,1}, {λ,β} = {2.2,6.5} and {λ,β} = {3.5,15} 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 figure 1). Finally, for {λ,β} = {3.5,50}, we observe self-trapping, as expected.

Zoom In Zoom Out Reset image size

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 4 correctly explain the overall trend of the numerical simulations and the separation between different spreading regimes in the parameter space.

5.2. Results of the Klein–Gordon 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 already been 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 the one hand, it allows for testing the generality of the result in the case where there is just one conserved quantity. This is highly nontrivial—especially for self-trapping—for which rigorous results have recently been 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) centred in the lattice, such that each site has equal energy $\mathcal {E}_j = \mathcal {E}=\mathcal {H}/L$. This is implemented by setting initial momenta of $p=\pm \sqrt {2 \mathcal {E}}$ with randomly assigned signs and zero coordinates. The values of initial energy densities $\mathcal {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 figure 6, while the expected spreading regimes are given in figure 2 (full upward triangles)¹¹. As one can see by comparing figure 6 with figure 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₂ according to power law m₂ ∼ t^γ with γ < 1. If self-trapping is avoided, all simulations show 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 $\mathcal {E}=0.05$ and $\mathcal {E}=0.01$, the characteristics of the weak chaos regime are observed, namely the exponent γ saturates around 1/3 (red and green curves in figure 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].

Zoom In Zoom Out Reset image size

For the two energy densities $\mathcal {E} = 0.055$ and 0.075, we expect strong chaos, with characteristics similar to the DNLS case. The simulation with $\mathcal {E} = 0.055$ (blue curves in figure 6) indeed exhibits the typical behaviour of the strong chaos scenario: the characteristic exponent γ increases up to the 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 $\mathcal {E}=0.075$, when intermediate strong chaos is masked due to partial self-trapping (magenta curves in figure 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. long-lasting overshooting with γ > 1/2, or significant slowing down to values γ < 1/3.

Finally, for $\mathcal {E}=1.0$ the dynamics enter the self-trapping regime, as our theory predicts. There a major part of the initial excitation stays localized, while the remainder spreads (yellow curves in figure 6). The participation numbers, therefore, do not grow significantly and 〈log₁₀ P〉 starts to level off at large time (figure 6, left panel, bottom, yellow curve). In contrast, the small spreading portion yields a continuous increase of the second moment m₂ (figure 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 〈ξ〉 (figure 6, right panel, bottom, yellow curve) drops to small values indicating a deep self-trapping regime. Note that similar behaviour has been observed before in purely random systems [11, 12]. Unusually large values of m₂ can be explained by local trapping–detrapping processes in the evolving wave packet. The corresponding dynamics is in strong non-equilibrium—its theoretical description is yet to be developed.

5.3. 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 centred around the site j = 0 has the form

$$\begin{equation} \psi_j(0)=C_1 \,\mathrm{e}^{-\frac{j^2}{2\sigma^2}}, \end{equation} \tag{ 13 }$$

where σ is a parameter controlling the width of the packet, while C₁ is a constant factor that can be determined by using the normalization condition $\sum_j |\psi_j|^2=1$. A TF wave packet is instead defined by

$$\begin{equation} \psi_j(0)=C_2 \sqrt{1-\frac{j^2}{R^2}} \end{equation} \tag{ 14 }$$

in the region where |j| < R and ψ_j = 0 otherwise. The parameter R is the TF radius characterizing the width of the distribution, while the constant C₂ is a normalization factor. These two distributions are of interest when considering ultracold bosons initially released from a harmonic trap in the Gross–Pitaevskii (GP) regime [33].

In figure 7, we show the time evolution in the DNLS model of the second moment of the expanding wave packet, 〈log₁₀ m₂〉 (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 figure 4¹². In particular, we use σ = 5 and R = 7.5, yielding a nonlinear frequency shift δ ≈ β/13. The values of β used in figure 7 are the same as those previously considered.

Zoom In Zoom Out Reset image size

From the comparison between the results of figures 7 and 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.

5.4. Application to cold atoms

In this work, we have considered the problem of the interplay between localization and interaction in 1D 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₂ ∼ 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 spreading 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].

ML acknowledges the Max Planck Institute for the Physics of Complex Systems, Dresden, for their hospitality. We are indebted to G Modugno and E Lucioni for fruitful discussions. This work was supported by the European Research Council through the QGBE grant and by the UPV/EHU under the UFI 11/55.