Localization in one-dimensional chains with Lévy-type disorder
LLocalization in one-dimensional chains with L´evy-type disorder
Sepideh S. Zakeri, ∗ Stefano Lepri,
2, 3, † and Diederik S. Wiersma
1, 4, 5 European Laboratory for Non-linear Spectroscopy (LENS),University of Florence, Via Nello Carrara 1, I-50019 Sesto Fiorentino, Italy Consiglio Nazionale delle Ricerche, Istituto dei Sistemi Complessi,via Madonna del Piano 10, I-50019 Sesto Fiorentino, Italy Istituto Nazionale di Fisica Nucleare, Sezione di Firenze,via G. Sansone 1, I-50019 Sesto Fiorentino, Italy Consiglio Nazionale delle Ricerche, Istituto Nazionale di Ottica, Largo Fermi 6, I-50125 Firenze, Italy Universit´a di Firenze, Dipartimento di Fisica e Astronomia,via G. Sansone 1, I-50019 Sesto Fiorentino , Italy
We study Anderson localization of the classical lattice waves in a chain with mass impuritiesdistributed randomly through a power-law relation s − (1+ α ) with s as the distance between two suc-cessive impurities and α >
0. This model of disorder is long-range correlated and is inspired by thepeculiar structure of the complex optical systems known as L´evy glasses. Using theoretical argu-ments and numerics, we show that in the regime in which the average distance between impurities isfinite with infinite variance, the small-frequency behaviour of the localization length is ξ α ( ω ) ∼ ω − α .The physical interpretation of this result is that, for small frequencies and long wavelengths, thewaves feel an effective disorder whose fluctuations are scale-dependent. Numerical simulations showthat an initially localized wavepacket attains, at large times, a characteristic inverse power-law frontwith an α -dependent exponent which can be estimated analytically. PACS numbers: 05.45.-a,05.60.-k,42.25.Dd,63.20.PwKeywords: Lyapunov exponent, L´evy distribution, localization, power-law correlation
I. INTRODUCTION
The spatial distribution of disorder plays a key roleon transport properties in disordered media, allowinganomalous laws like superdiffusion to arise. This inter-esting topic has been basis of numerous theoretical andexperimental studies in a wide variety of complex sys-tems. Recently, one of the studies reports on realiza-tion of some engineered materials named L´evy glasses inwhich light rays propagate through an assembly of trans-parent microspheres embedded in a scattering medium[1]. If the diameter of the microspheres φ is designedto have a power-law distribution p ( φ ) ∼ φ − ( α +1) , where α is the so-called L´evy exponent defining degree of theheterogeneity of the system, light can indeed perform su-perdiffusion. Although a random walk model is rich andprecise to describe light transport through L´evy glassesand successfully explains the experimental observations[2], it cannot address wave properties such as polariza-tion and interference. Therefore, an open question tounderstand is how interference of waves can affect thepropagation and what might be the possible role of An-derson localization [3, 4] in such materials.Although our main motivation stems from the abovedescribed setup, it is worth mentioning some relatedstudies of systems with L´evy-type disorder that includetransport in quantum wires [5], photonic heterostructures[6], and disordered electronic systems [7, 8]. The purpose ∗ [email protected]fi.it † [email protected] of this work is to provide a framework to investigate local-ization in power-law correlated disordered systems and toillustrate how the localization length is affected by dif-ferent characteristic features of the system such as fre-quency, degree of heterogeneity, disorder strength, etc.To address the above issues we consider a very simplemodel, a harmonic chain of coupled oscillators with ran-dom impurities separated by random distances s havinga power-law distribution p ( s ) ∝ s − ( α +1) . This model ofdisorder is clearly inspired by the peculiar structure ofL´evy glasses. As usual, the one-dimensional case allowsfor a detailed study. In particular, the Lyapunov expo-nent γ ( ω ), the inverse of the frequency-dependent local-ization length ξ ( ω ) = γ ( ω ) − can be computed straight-forwardly. Our main result is that, in the small-frequencyregime, the scaling relation γ ( ω ) ∝ ω α is numerically esti-mated for the Lyapunov exponent in the range 1 ≤ α ≤ (cid:104) s (cid:105) is infinite. In-stead, for α >
2, the usual scaling γ ( ω ) ∝ ω , typical ofrandom uncorrelated disorder [9] is recovered.The model and some theoretical arguments, based onthe Hamiltonian map formulation of the transfer method,will be presented in section II. We then present and dis-cuss the complete numerical steady-state analysis in sec-tion III of this paper. For a more comprehensive under-standing of the transport in such disordered media, weinvestigate and report in section IV the time evolutionof an initially localized wave packet. In particular, weconsider the time and disorder averaged energy profile (cid:104) e n ( t ) (cid:105) . Here, we show that in the range 1 ≤ α < α -dependent power-law exponent. In contrast, for therange α ≥
2, depending on the type of the initial excita- a r X i v : . [ c ond - m a t . s t a t - m ec h ] D ec tion, the power-law exponent is independent of α . More-over, we pay close attention to the time evolution of thedifferent moments m ν ( t ) of the (cid:104) e n ( t ) (cid:105) and discuss the α -dependent properties in the range 1 ≤ α <
2. Finally,we summarize our results in the concluding section V.
II. THEORETICAL DESCRIPTIONA. Model: 1D L´evy-type disordered lattice
We considered a one-dimensional harmonic disorderedlattice of n sites with masses m n . The governing equa-tions of motion are m n ¨ u n = k ( u n +1 − u n ) + k ( u n − − u n ) , (1)where k is spring constant ( k (cid:54) = 0), u n is displacement ofthe n th site of the lattice from its equilibrium positionand n = 1 , , ..., N .In the present work, we consider dichotomic disorderwhereby m n assumes two possible values, either m as thebackground mass of the lattice or M to represent massof the defects. Illustrated in Fig. 1, disordered lattices oflength N were constructed by arrangement of N d defectsthrough a power-law distribution p ( s ) ∝ s − ( α +1) , where s was the random number of lattice sites with mass m between two successive defects. As a direct consequenceof such choice of p ( s ), the mass of the entire lattice wasthe parameter with α -dependent statistical properties.Obviously, its average density can be written as (cid:104) m (cid:105) = ρM + (1 − ρ ) m, (2)where ρ = N d N is the fraction of defects on the lattice.In the range 1 ≤ α < (cid:104) m (cid:105) is finite. On the otherhand, for α < ρ vanishes in the thermodynamic limit N → ∞ since the average distance between consecutivedefects diverges in this regime. In this paper, we mainlyfocus on the case α > (cid:104) s (cid:105) is finite) and willcomment only briefly on the case α < α , fundamentally differentregimes of transport are achieved: indeed walkers su-perdiffuse for α <
1. Several distinguished features of thequenched nature of disorder like the importance of initialconditions and the consequences on higher-order statis-tics of the diffusive process are discussed in the abovementioned papers. It can be thus envisaged that alsolocalization properties may display unusual features.
B. Statistical properties of the disorder
Using the algorithm described in [16], positive integerpower-law random numbers were generated to construct disordered lattices. For a preliminary statistical char-acterization, we computed the ensemble-averaged powerspectrum: S ( k ) = 1 N (cid:42)(cid:12)(cid:12) N (cid:88) n =1 m n exp( − ikn ) (cid:12)(cid:12) (cid:43) , (3)where the average is over different realization of the dis-order. The low-wavenumber behavior of S provides in-formation on the large-scale decay law of the disordercorrelation function which will be necessary for the sub-sequent analysis. Fig. 2 indicates that, for α < S has apower-law singularity | k | − Θ for small k . This implies thatthe disorder correlation is indeed decaying as an inversepower of the relative distance with an exponent 1 − Θ.The exponent will be important in the subsequent anal-ysis: fitting of the numerical data suggests the followingrelation with the L´evy parameter α (see the Inset of Fig.2) Θ( α ) = (cid:40) α if 0 < α < − α if 1 < α < . (4)A theoretical justification of the above relation can begiven based on the similarity of this process with the cor-relation of a L´evy walk as discussed by Geisel in Ref. [17].Note that in all the examined cases Θ <
1, so the powerspectrum is integrable and the associated process can beconsidered as stationary for any α .For what concerns the dependence on the latticesize (data not shown), by comparing data for differentlengths, we found that for α > N in-dependent, while for α < N roughly like N α − . This is because thedensity of defects vanishes in the thermodynamic limitand almost-ordered realizations dominate the statisticalaverages.To conclude this subsection, let us mention that a sim-ilar model has been studied in Ref. [18] (see also therelated work [19]). The main differences with our workis that the sequence of masses is generated as a trace ofa fractional Brownian motion that, by construction, hasalso a a power-law singularity at small wavenumber. Itwas there shown that in the non-stationary case (Θ > C. Transfer matrix approach
As it is well known, the localization properties can bestudied by the transfer matrix method [20]. Assumingthat u n oscillates harmonically in time at an angular fre-quency ω , so that u n = v n e iωt , Eq. (1) can be written asan eigenvalue problem − m n ω v n = k ( v n +1 − v n ) + k ( v n − − v n ) , (5) FIG. 1. (color online) Schematic illustration of the model used to study localization in a one-dimensional lattice with L´evy-typedisorder. Small blue balls represent the background mass of the lattice m and large red balls are the defects with mass M distributed on the lattice through a power-law relation p ( s ) ∝ s − ( α +1) , where s is the number of lattice sites between twosuccessive defects and the index j = 1 , ..., N d . An initial displacement w is given to the first site of the lattice to initiateiteration of the transfer matrix over the entire lattice. α Θ -5 -4 -3 -2 -1 k/2 π -3 -2 -1 P o w e r s pe c t r u m S ( k ) ( a r b . un i t s ) α=0.75α=1.25α=1.5α=1.75 FIG. 2. (color online) Power spectra, Eq. (3), for differentvalues of the exponent α ; each spectrum is obtained for se-quences of length 2 over an ensemble of 10 realizations.The inset shows the dependence of the exponent of the small-frequency singularity | k | − Θ on α . that, upon defining w n ≡ (cid:18) v n v n − (cid:19) ; T n ≡ (cid:18) − m n ω /k −
11 0 (cid:19) , for n = 2 , , ..., N − w n +1 = T n w n , (6)where T n is a 2 × n th site onthe lattice. By iteration of the transfer matrix over theentire lattice and applying the appropriate boundary con-ditions of the system, solutions of Eq. (5) as a functionof frequency ω can be obtained.The central quantity to be computed is the Lyapunovexponent γ ( ω ) which gives the inverse of the localizationlength ξ ( ω ). It is convenient to define a new parameter R n = v n /v n − . Inserting R n into Eq. (5), one obtainsthe following recursive relation [21] R n +1 = 2 − m n ω k − R n . (7) Eq. (7) can be interpreted as a ”discrete time” stochasticequation. The mass m n plays the role of a noise source(with bias) whose strength is gauged by the frequency ω . In the present model, Lyapunov exponent γ ( ω ) as theinverse of localization length can be computed by ξ ( ω ) − = γ ( ω ) = (cid:104) ln R n (cid:105) . (8)Also the integrated density of states I ( ω ) follows fromnode counting arguments, i.e. I ( ω ) = f , where f is thefraction of negative R n values.For deeper analysis of the systems we are interestedin, let us rewrite Eq. (5) in a different notation, intro-ducing the new variables q n = v n and p n = q n − q n − ,the transfer map can be reformulated as follows. Letting ω m n = Ω (1+ µ n ), Ω = ω (cid:104) m (cid:105) and µ n = m n / (cid:104) m (cid:105)− (cid:104) µ (cid:105) is finite at leastfor α >
1. Plugging q n and p n into Eq. (5) and setting k = 1, one finds the following two-dimensional map [22] p n +1 = p n − Ω (1 + µ n ) q n (9) q n +1 = q n + p n +1 . (10)The following transformation relations can be intro-duced by using the canonical variables ( r n , θ n ) where r n and θ n are amplitude and phase of the eigenvector at site n , respectively: p n = √ r n sin θ n ; q n = (cid:114) r n cos θ n . (11)Substituting Eqs. (11) back into Eqs. (9) and (10) andneglecting terms of the order Ω and higher, one thenobtains the following map (cid:18) r n +1 r n (cid:19) = 1 − µ n Ω sin 2 θ n + O (Ω ) (12)tan θ n +1 = tan( θ n − Ω) − Ω µ n . (13)Eq. (13) describes evolution of the phase as it is per-turbed by disorder and can be approximated as [23] θ n +1 = θ n − Ω + Ω µ n sin θ n . (14)Long-range correlations of µ n lead to sub-diffusive phasefluctuations. Within the same approximation, the Lya-punov exponent is obtained by expanding the logarithm γ = 12 (cid:104) log (cid:18) r n +1 r n (cid:19) (cid:105) ≈ −
12 Ω (cid:104) µ n sin 2 θ n (cid:105) + 14 Ω (cid:104) µ n sin θ n (cid:105) + . . . (15)The main issue is to estimate the correlators in Eq.(15). Following the same path as in ref. [23] (see inparticular section (5.2)), it can be shown that (cid:104) µ n sin 2 θ n (cid:105) = 2Ω ∞ (cid:88) m =1 (cid:104) µ n µ n − m (cid:105) cos(2Ω m ) (16) (cid:104) µ n sin θ n (cid:105) ≈ (cid:104) µ n (cid:105) . (17)We should get something like γ = 18 Ω (cid:104) µ n (cid:105) (cid:32) ∞ (cid:88) m =1 (cid:104) µ n µ n − m (cid:105)(cid:104) µ n (cid:105) cos(2Ω m ) (cid:33) . (18)Note that this formula agrees with the long-wavelengthlimit of Eq.(22) in [24] which was obtained for weak dis-order (where Ω ≈ k ). The sum is basically the Fouriertransform of the correlation function, i.e. the power spec-trum of the sequence of masses studied above. We thusexpect three cases: • For short-range correlations as in the case α > γ ( ω ) ∝ ω (cid:104) m (cid:105) − (cid:104) m (cid:105) (cid:104) m (cid:105) . (19) • For long-range correlations and in the case 1 < α <
2, the fluctuation spectrum diverges at small fre-quencies so the second term in Eq. (18) dominatesso that γ ∼ Ω S (2Ω) . (20)Using the result Eq. (4) that the spectrum goes as k α − and recalling the definition of Ω this meansthat the Lyapunov exponent should be propor-tional to ω α . • Furthermore, if we assume that the same argumentholds also for α <
1, we may conclude that γ ∼ ω − α /N − α , which means that γ vanishes for N →∞ . This argument may not be entirely correct sincein this case, (cid:104) µ n (cid:105) also vanishes. To conclude this section, we mention that John andStephen [25] studied a related model of a classical wave inthe continuum limit where the mass fluctuations µ werequenched random variables but power-law correlated inspace (cid:104) µ ( x ) µ ( x (cid:48) ) (cid:105) ∼ ( x − x (cid:48) ) − m . According to their Eq.(7.5), in 1D the Lyapunov exponent should be γ ( ω ) ∼ (cid:40) ω if m > / ω − m if m < / . (21)Although the models are pretty similar at large scales,our result is different from this last estimate. The originof the discrepancy is not clear at this stage, however it istrue that the approach in ref. [25] relies on some approx-imate self-consistent calculations while the calculationsreported here are somehow exact. III. NUMERICAL RESULTS
In this section, numerical results of the model de-scribed in the previous section are presented. We havecarefully studied the dependence of the Lyapunov ex-ponent on different characteristic parameters of the sys-tems. Data analysis was performed mainly for systemswith 1 ≤ α <
2, however, some computations were donefor systems with α ≥ k = 1 and, unless otherwise stated, the mass ratiowas fixed to M/m = 3.
A. Lyapunov exponent
According to the definition described in Eq. (8), theLyapunov exponent is the exponential growth or decayrate of a vector in the limit N → ∞ .Due to the long-range spatial correlations (especially inthe regime 1 ≤ α < γ could vary from a single realization to anotherand usually requires a long series of iterations. In orderto ensure that our obtained results were independent ofthe disorder realization (i.e. that self-averaging holds), aseries of lattices with different number of defects N d werestudied and compared. Results led us to the conclusionthat for N d ≥ × , fluctuations of γ ( ω ) versus N d are relatively small (better than (cid:39) ± α parameters. Clearly, as α approaches to2, even smaller values of N d can yield a fast convergence γ ( ω ). However, to improve stability of our results evenfurther, N d = 5 × defects were used to construct ourdesired disordered systems. We pursue this number con-sistently in the rest of this paper for the analysis of verylarge lattices (asymptotic limit). In all the studied cases,same initial displacement w must be assigned to startiteration of the transfer matrix. Note that Lyapunov ex-ponent has this intrinsic property to be independent of −2 −1 −7 −6 −5 −4 −3 Frequency ω L y apuno v e x ponen t γ ( ω ) α =1.3 α =1.5 α =1.7 α =2.5 α =3 α σ AnalyticalFit
FIG. 3. (color online) Variations of the Lyapunov exponent γ ( ω ) versus ω for systems with different α parameter, setting N d = 5 × , M/m = 3 and k = 1. Lyapunov exponentswere normalized after N it = 5 × iterations of the transfermatrix. Frequency ω was varied from 0 .
01 to 0 .
3. The in-set shows the scaling exponent σ (of the frequency-dependentLyapunov exponent γ ( ω ) ≈ ω σ ) versus α . Numerical resultsfor the scaling are obtained by power-law fits on the dataset for Lyapunov exponent and an excellent agreement is ob-tained with the theoretical predictions γ ( ω ) ≈ ω α in the low-frequency regime for the range 1 < α < w and our numerical results were in excellent agreementwith that.Since we mostly focus on the behaviour of the Lya-punov exponent in the small frequency regime, we firstlychecked that in all the presented cases the integrateddensity of states is linearly changing with the frequency, I ( ω ) ∼ ω . The physical interpretation of this fact issimply that the spectrum is effectively equivalent to anordered chain with a renormalized wave speed.In Fig. 3, variations of the Lyapunov exponent γ ( ω )versus frequency ω for systems with different α parameterare shown in doubly logarithmic scale. In each system,as the frequency was increased from 0 .
01 to 0 .
3, greatervalues were obtained for Lyapunov exponent and accord-ingly, localization length was decreased. As illustrated inthe inset of Fig. 3, the scaling relation predicted in theprevious section is in excellent agreement with the dataobtained by power-law fitting on the small-frequency re-gion of the curves. In some cases, we also checked thatthe fit that includes the leading order correction ω fol-lows the data in the whole range, giving further supportto the theoretical arguments above.Another interesting physical effect that should be as-sessed is the dependence of the Lyapunov exponent ondisorder strength or equivalently mass ratio M/m in ourmodel. In an intuitive picture, one expects that highermass ratio leads to stronger disorder, shorter localiza-tion length and thereby greater Lyapunov exponents. Inaddition to that, the scaling of γ ( ω ) ∝ ω σ intrinsically −2 −1 −7 −6 −5 −4 −3 −2 Frequency ω L y apuno v e x ponen t γ ( ω ) M/m=0.5M/m=2M/m=3M/m=5
FIG. 4. (color online) Comparison of the Lyapunov exponent γ ( ω ) versus ω for systems with α = 1 . N = 14648365but different mass ratios M/m . Other parameters such as N d , N it and the frequency range are the same as those usedto generate Fig. 3. −4 Levy parameter α L y apuno v e x ponen t γ ω =0.05 ω =0.24 α γ / γ m a x FIG. 5. (color online) Dependence of the Lyapunov exponenton the L´evy exponent α at two fixed frequencies ω = 0 . ω = 0 .
24 (orange triangles). The insetshows that by increase of the frequency, peak of the Lyapunovexponent is shifted to right. For better illustration, two curvesare normalized by their maxima. depends on transport properties of the system and henceshould be independent of the mass ratio. In other words,the mass ratio should at most change the prefactor inthe above scaling relation. To confirm this prediction,four different mass ratios
M/m = 0 . , , , α = 1 . N = 14648365composed of N d = 5 × defects. As shown in Fig.4, scaling of the Lyapunov exponent with frequency ω (slope of the curves) is independent of the mass ratios.However, larger values of M/m yield an increase in theLyapunov exponent by orders of magnitude, indicatingthat the prefactor changes rapidly with
M/m .Fig. 5 shows and compares dependence of the Lya-punov exponent on the L´evy exponent α at two fixedfrequencies ω = 0 .
05 and ω = 0 .
24. Upon decreasing α , γ was initially increased. This is in agreement with theintuition that the localization is enhanced by increasingfluctuations of the distances between defects. However,by approaching α = 1, γ started decreasing and a peakwas observed in both curves. We can therefore concludethat in the range 1 . ≤ α ≤ .
3, a transition occurs inthe behaviour of the Lyapunov exponent. The value of α , at which the transition happens, depends slightly onthe frequency ω . According to our numerical data, thetransition appeared at α = 1 . ω and at α = 1 . ω (see the inset in Fig. 5, showing the Lyapunov ex-ponents normalized by their maximums). Moreover, noappreciable dependence on the lattice length is observedin this regime. Although it seems plausible that γ van-ishes below α = 1, the origin of the peak is not clear. Apossibility is that sample-to-sample fluctuations becomeso important that averaging over a much larger ensembleof lattices would be required. B. Eigenfunctions
Moreover, it is interesting to look at eigenfunctions ofa L´evy-type disordered system. In general, spatial partof the frequency dependent solutions of Eq. (1) can beobtained by solving( K + ω Λ ) u = 0 , (22)where Λ is the diagonalized mass matrix Λ =diag( m n )and K is the matrix of spring constants expressed by K = − k k . . .k − k k . . . k − k k . . . ... . . .0 . . . k − k N × N Stationary eigenfunctions of disordered systems can beobtained by numerically solving Eq. (22). For a systemwith α = 1 . N = 1187 consisted of N d = 200defects, eigenfunctions u n at three different frequenciesare presented in Fig. 6. Top panel shows distribution ofmass m n on the entire system which is a step-functionof m n = 1 or 3. Eigenfunctions of the systems at threedifferent frequencies are depicted in the other panels. Ac-cording to our results, at ω = 0 . ω = 1 . ω = 0 . m n | u n | | u n | | u n | n ω =0.9083 ω =1.6853 ω =0.47925 FIG. 6. (color online) Top panel shows distribution of themass on a system with α = 1 . N = 1187 com-posed of N d = 200 defects (red curve). Other panels (from topto bottom) show eigenfunctions of the system at ω = 0 . ω = 0 . ω = 1 . C. Phase dynamics
Next we present our numerical results for the spatialevolution of the phase that, using Eq. (11) and the defini-tion p n = q n − q n − , can be computed during the iterationof the transfer matrix by the following relation θ n = tg − (cid:18) q n − q n − Ω q n (cid:19) . (23)The computed phase has to be unwrapped to the realaxis. Top panel in Fig. 7 illustrates variations of θ n insystems of equal length N = 4 × and different L´evyexponents α = 1 . , . , . ω = 0 . a ofthe phase ( a ≈ Ω) was numerically computed and re-moved by defining φ n = θ n − an . In order to studyevolution of φ n over the entire system, the lattice was di-vided into 50 sub-lattices of length n s = 8 × . For eachsub-lattice, the quantity ( φ n − φ ) , where n = 1 , .., n s was computed. Averaged over all the sub-lattices, evo-lution of the rms-phase (cid:104) φ n (cid:105) was consequently obtainedas shown in the bottom panel of Fig. 7. According to −10−50 x 10 n θ n −2 n s < φ n2 > α =1.5 α =1.7 α =2.5 α =1.5 α =1.7 α =2.5 FIG. 7. (color online) Top panel shows variations of the un-wrapped phase θ n versus n for three different systems with α = 1 . , . , . N = 4 × at a fixed frequency ω = 0 .
1. In the bottom panel, variations of (cid:104) φ n (cid:105) over a sub-lattice of length n s = 8 × are presented and compared forthe studied systems. Note that Ω . = 0 . . = 0 . . = 0 . the data, growth of (cid:104) φ n (cid:105) is faster in systems with smallervalues of α . This can be explained by the fact that disor-der fluctuations become stronger which is in agreementwith the superdiffusive behaviour of the correspondingrandom-walk problem [13]. IV. EVOLUTION OF WAVEPACKETS
In this section, we investigate the consequences of theabove results on the spreading of an initially localizedperturbation on an infinite lattice. Although every sin-gle eigenstate is exponentially localized, the wave propa-gation can be non trivial as low-frequency states form acontinuum with arbitrarily large localization lengths.The starting point of this analysis is a Hamiltonianformulation of the system: H = (cid:88) n [ P n m n + 12 k ( u n +1 − u n ) ] , (24)where u n and P n are displacement and momentum ofthe mass m n , respectively. From the Hamiltonian, timeevolution of any excitation of interest can be derived.Evidently, local energy at each site of the lattice at time t can be written as: −8 −6 −4 −2 |n−n | A v e r age ene r g y den s i t y α =2.5 α =1.5 −5 −4 −3 |n−n | −5/3 (a) −8 −6 −4 −2 |n−n | A v e r age ene r g y den s i t y α =2.5 α =1.5 −6 −4 |n−n | −3 (b) FIG. 8. (color online) Spreading of the energy density after t = 3 × averaged over 2 × realizations. Results arecomparison for two systems with equal lengths N = 8 × but different L´evy exponents α = 1 . , . n = N/
2. (a) Momentum excitation with B = 2. Inthe inset, obtained energy density data of the system with α = 1 . | n − n | − / in the asymptotic regime. (b)Displacement excitation with A = 2. In the inset, obtainedenergy density data of the system with α = 1 . | n − n | − in the asymptotic regime. e n ( t ) = e kinn ( t ) + e potn ( t ) (25) e kinn ( t ) = 12 m n ˙ u n ( t ) (26) e potn ( t ) = 12 k [ u n ( t ) − u n − ( t )] u n ( t ) − k [ u n +1 ( t ) − u n ( t )] u n ( t ) , (27)where e kinn ( t ) and e potn ( t ) represent kinetic and potentialenergies, respectively. Together with the equation of mo-tion Eq. (5), these relations constitute the set of govern-ing equations of one-dimensional linear discrete systems.In our analysis, disordered lattices of length N =8 × with an excitation at the center n = N/ u n (0) = Aδ n,n time < m ν ( t ) > (a) time < m ν ( t ) > (b) time < m ν ( t ) > (c) time < m ν ( t ) > (d) FIG. 9. (color online) Time evolution of the moments (cid:104) m ν ( t ) (cid:105) averaged over 2 × realizations for systems with length N = 8 × . Different colours represent different index ν varying from ν = 0 . ν = 4 (dashed red) withstep ∆ ν = 0 .
5. Momentum excitation ( B = 2) was appliedat n = N/ α = 1 . α = 2 .
5. Displacement excitation ( A = 2) was applied at n = N/ α = 1 . α = 2 . and P n (0) ≡
0; (ii) momentum excitation with P n (0) = Bδ n,n and u n (0) ≡
0. As it is known, these classes ofinitial condition yield different asymptotic properties dueto the fact that the mode energy distribution is differentin the two cases [26]. Based on the conservative dynamicsof the system and by using a fourth-order symplectic al-gorithm [27], the set of coupled governing equations werenumerically solved. We mainly focused on lattices with α = 1 . α = 2 . N d on the considered lattices were different as N was kept fixed. Obtained results were then averaged over2 × realizations. Fig. 8 presents spread of the aver-aged energy density over the two systems after t = 3 × for both momentum and displacement excitations. In ad-dition to a faster spread of energy, a broadening appearedaround the excitation for the system with α = 1 . ξ α ( ω ) ∝ ω − α in the small-frequency regime ( ω → (cid:104) e n ( t ) (cid:105) ∝ (cid:90) ∞ ω α e −| n − n | ω α dω, (28)which asymptotically decays as | n − n | − (1+3 /α ) . For amomentum excitation, (cid:104) e n ( t ) (cid:105) can be similarly describedby (cid:104) e n ( t ) (cid:105) ∝ (cid:90) ∞ ω α e −| n − n | ω α dω, (29) with asymptotic behaviour as | n − n | − (1+1 /α ) . Theinsets in Fig. 8 show the excellent agreement betweenour numerical results and the analytical predictionsof the α -dependent power-law decay. In the case ofmomentum excitation, data for α = 1 . | n − n | − / with a goodness of 95% as shown in theinset in Fig. 8(a). For the displacement excitation, datafor α = 1 . | n − n | − with a goodness of95% which is illustrated in Fig. 8(b).Moreover, another interesting entity is time evolutionof the moments of the Hamiltonian which can be definedas [26, 28] m ν ( t ) = (cid:80) n | n − n | ν e n ( t ) H , (30)where ν is a positive number and m ν ( t ) represents the ν th moment of the Hamiltonian at time t . The sec-ond moment quantifies degree of the spreading of thewavepacket.In Fig. 9, time evolution of different moments (cid:104) m ν ( t ) (cid:105) of a wavepacket are shown and compared for the twosystems. Beyond the clear signature of the growth of thedisorder averaged moments over time, it was observedthat around t ≈ , the growth exponent suddenlychanged in the system with α = 1 .
5. At shorter times,however, we may assume that correlation effects arenot yet strong and system behaved like a short-rangecorrelated one. Note that such effect is even moreemphasized in the case of displacement excitation. Tonumerically estimate the scaling of (cid:104) m ν ( t ) (cid:105) ∝ t β ( ν ) , apower-law fit was performed on the averaged moments.The obtained results shown in Fig. 10, reveal that β ( ν )is smaller in the system with α = 1 . α = 2 . β ( ν ) is different insystems with different L´evy parameters. In an effortto determine a general analytic relation for β ( ν, α ) thatdescribes dependence on both parameters ν and α , theasymptotic forms of Eqs. (28) and (29) were insertedinto Eq. (30). Later, by defining an upper cutoff forthe summation in the numerator at the ballistic distance | n − n | = ct , it can be shown that for displacement ex-citation β ( ν, α ) = (cid:40) ν − α , αν >
30 , αν < , (31)while for momentum excitation β ( ν, α ) = (cid:40) ν − α , αν >
10 , αν < . (32)As shown in Fig. 10, there is an excellent agreement be-tween the theoretical estimates, (according to Eqs. (31)and (32)), and the numerical data obtained by power-lawfits on (cid:104) m ν ( t ) (cid:105) in the long-time regimes. ν β ( ν ) α =1.5 α =2.5analytical ν β ( ν ) α =1.5 α =2.5analytical FIG. 10. (color online) Comparison of β ( ν ) for the two sys-tems with L´evy exponent α = 1 . α = 2 . (cid:104) m ν ( t ) (cid:105) . Dotted black lines show the numerical solu-tions of analytical Eqs. (31) and (32) for each system. (a)Momentum excitation ( B = 2) (b) Displacement excitation( A = 2). V. SUMMARY AND CONCLUSION
In this paper, we have studied the localization proper-ties of classical lattice waves in the presence of power-lawcorrelated disorder. The model is mathematically sim-ple and the choice of the disorder is motivated by thesamples known as L´evy glasses [1]. In the first part, wecomputed frequency-dependent localization length andanalyzed how different characteristic parameters can af-fect it. Using theoretical arguments and numerics, wehave shown that in the regime 1 ≤ α < γ ( ω ) ∼ ω α so that γ mostly decreases as α increases and this trendremains the same even in the regime α >
2. The physicalinterpretation of this result is that, for small frequenciesand long wavelengths, the waves see an effective disorderwhose fluctuation is scale-dependent. This stems fromthe fact that the disorder is, by construction, scale-free.Large ordered regions can exist inside the system de-pending on the L´evy exponent α . As expected, weshowed that eigenfunctions of the system are localizedwhere the density of disorder is high and extended in thelarge dilute spaces inside the system. Moreover, due tothe superdiffusive nature of the transport in the range1 ≤ α <
2, we showed that the spatial growth of theroot-mean-squared phase (cid:104) φ n (cid:105) is greater in systems with smaller values of α .The anomalous localization properties reflects in theproblem of wavepacket spreading. An initially local-ized perturbation attains at large times a characteris-tic power-law decay as in the uncorrelated case [28, 29].However in the present case, the exponent of the decaydepends on the exponent α of the disorder distribution.This implies that the disorder averaged momenta (cid:104) m ν ( t ) (cid:105) must diverge with time, although the initial local energyexcitation does not spread completely. The predictionsare clearly supported by the numerical results. As a mat-ter of fact, it should be remarked that consideration of m ( t ) alone is not sufficient to conclude that the energydiffuses. Clearly, the origin of such a growth is drasti-cally different to a genuinely (possibly anomalous) dif-fusive processes observed in nonlinear oscillator chains[30–33] and just stems from the slowly decaying tails.In the present work we mostly focused on the case α >
1. As argued in section II C, the case α <
ACKNOWLEDGMENTS
Authors wish to thank Mario Mulansky for readingthe manuscript and the financial support from the Eu-ropean Research Council (FP7/2007-2013), ERC GrantNo. 291349. [1] P. Barthelemy, J. Bertolotti, and D. S. Wiersma, Nature , 495 (2008).[2] J. Bertolotti, K. Vynck, L. Pattelli, P. Barthelemy,S. Lepri, and D. S. Wiersma, Adv. Funct. Mater. ,965 (2010).[3] P. W. Anderson, Phys. Rev. , 1492 (1958).[4] P. Sheng, Introduction to wave scattering, localizationand mesoscopic phenomena , vol. 88 (Springer Verlag,2006).[5] F. Falceto and V. A. Gopar, EPL (2010).[6] A. A. Fernandez-Marin, J. A. Mendez-Bermudez, and V. A. Gopar, Phys. Rev. A , 035803 (2012).[7] P. J. Wells, J. d’Albuquerque Castro, and S. de Queiroz,Phys. Rev. B , 35102 (2008).[8] A. Iomin, Phys. Rev. E , 062102 (2009).[9] H. Matsuda and K. Ishii, Prog. Theor. Phys. Suppl. ,56 (1970).[10] E. Barkai, V. Fleurov, and J. Klafter, Phys. Rev. E ,1164 (2000).[11] C. Beenakker, C. Groth, and A. Akhmerov, Phys. Rev.B , 024204 (2009).[12] P. Barthelemy, J. Bertolotti, K. Vynck, S. Lepri, and D. S. Wiersma, Phys. Rev. E , 011101 (2010).[13] R. Burioni, L. Caniparoli, and A. Vezzani, Phys. Rev. E , 060101 (2010).[14] P. Buonsante, R. Burioni, and A. Vezzani, Phys. Rev. E , 021105 (2011).[15] P. Bernab´o, R. Burioni, S. Lepri, and A. Vezzani, Chaos,Solitons & Fractals , 11 (2014).[16] L. Devroye, Non-uniform random variate generation (Springer-Verlag, New York, 1986).[17] T. Geisel, in
L´evy flights and related topics in physics (Springer, 1995), pp. 151–173.[18] F. A. B. F. de Moura, M. Coutinho-Filho, E. Raposo,and M. Lyra, Phys. Rev. B , 012202 (2003).[19] A. Costa and F. De Moura, J. Phys.: Condens. Matter , 065101 (2011).[20] A. Crisanti, G. Paladin, and A. Vulpiani, Products ofrandom matrices (Springer, 1993).[21] B. Derrida and E. Gardner, J. Phys. France , 1283(1984).[22] F. Izrailev, S. Ruffo, and L. Tessieri, J. Phys. A: Math.Gen. , 5263 (1998).[23] F. Izrailev, A. Krokhin, and N. Makarov, Phys. Rep. , 125 (2012).[24] I. F. Herrera-Gonz´alez, F. M. Izrailev, and L. Tessieri,EPL (Europhysics Letters) , 14001 (2010).[25] S. John and M. J. Stephen, Phys. Rev. B , 6358 (1983).[26] P. K. Datta and K. Kundu, Phys. Rev. B , 6287 (1995).[27] R. I. McLachlan and P. Atela, Nonlinearity , 541 (1992).[28] S. Lepri, R. Schilling, and S. Aubry, Phys. Rev. E ,056602 (2010).[29] P. Krapivsky and J. Luck, J. Stat. Mech: Theory Exp.p. P02031 (2011).[30] G. Zavt, M. Wagner, and A. L¨utze, Phys. Rev. E ,4108 (1993).[31] P. Cipriani, S. Denisov, and A. Politi, Phys. Rev. Lett. , 244301 (2005).[32] L. Delfini, S. Denisov, S. Lepri, R. Livi, P. K. Mohanty,and A. Politi, Eur. Phys. J. Spec. Top , 21 (2007).[33] M. Mulansky and A. Pikovsky, EPL (Europhysics Let-ters) , 10015 (2010).[34] P. Gaspard and X. J. Wang, Proceedings of the NationalAcademy of Sciences85