Localization in onedimensional chains with Lévytype disorder
Abstract
We study Anderson localization of the classical lattice waves in a chain with mass impurities distributed randomly through a powerlaw relation with as the distance between two successive impurities and . This model of disorder is longrange correlated and is inspired by the peculiar structure of the complex optical systems known as Lévy glasses. Using theoretical arguments and numerics, we show that in the regime in which the average distance between impurities is finite with infinite variance, the smallfrequency behaviour of the localization length is . The physical interpretation of this result is that, for small frequencies and long wavelengths, the waves feel an effective disorder whose fluctuations are scaledependent. Numerical simulations show that an initially localized wavepacket attains, at large times, a characteristic inverse powerlaw front with an dependent exponent which can be estimated analytically.
pacs:
05.45.a,05.60.k,42.25.Dd,63.20.PwI Introduction
The spatial distribution of disorder plays a key role on transport properties in disordered media, allowing anomalous laws like superdiffusion to arise. This interesting topic has been basis of numerous theoretical and experimental studies in a wide variety of complex systems. Recently, one of the studies reports on realization of some engineered materials named Lévy glasses in which light rays propagate through an assembly of transparent microspheres embedded in a scattering medium Barthelemy et al. (2008). If the diameter of the microspheres is designed to have a powerlaw distribution , where is the socalled Lévy exponent defining degree of the heterogeneity of the system, light can indeed perform superdiffusion. Although a random walk model is rich and precise to describe light transport through Lévy glasses and successfully explains the experimental observations Bertolotti et al. (2010), it cannot address wave properties such as polarization and interference. Therefore, an open question to understand is how interference of waves can affect the propagation and what might be the possible role of Anderson localization Anderson (1958); Sheng (2006) in such materials.
Although our main motivation stems from the above described setup, it is worth mentioning some related studies of systems with Lévytype disorder that include transport in quantum wires Falceto and Gopar (2010), photonic heterostructures FernandezMarin et al. (2012), and disordered electronic systems Wells et al. (2008); Iomin (2009). The purpose of this work is to provide a framework to investigate localization in powerlaw correlated disordered systems and to illustrate how the localization length is affected by different characteristic features of the system such as frequency, degree of heterogeneity, disorder strength, etc. To address the above issues we consider a very simple model, a harmonic chain of coupled oscillators with random impurities separated by random distances having a powerlaw distribution . This model of disorder is clearly inspired by the peculiar structure of Lévy glasses. As usual, the onedimensional case allows for a detailed study. In particular, the Lyapunov exponent , the inverse of the frequencydependent localization length can be computed straightforwardly. Our main result is that, in the smallfrequency regime, the scaling relation is numerically estimated for the Lyapunov exponent in the range , i.e. when the variance of the distances is infinite. Instead, for , the usual scaling , typical of random uncorrelated disorder Matsuda and Ishii (1970) is recovered.
The model and some theoretical arguments, based on the Hamiltonian map formulation of the transfer method, will be presented in section II. We then present and discuss the complete numerical steadystate analysis in section III of this paper. For a more comprehensive understanding of the transport in such disordered media, we investigate and report in section IV the time evolution of an initially localized wave packet. In particular, we consider the time and disorder averaged energy profile . Here, we show that in the range , the asymptotic tails of the energy profile decay with an dependent powerlaw exponent. In contrast, for the range , depending on the type of the initial excitation, the powerlaw exponent is independent of . Moreover, we pay close attention to the time evolution of the different moments of the and discuss the dependent properties in the range . Finally, we summarize our results in the concluding section V.
Ii Theoretical description
ii.1 Model: 1D Lévytype disordered lattice
We considered a onedimensional harmonic disordered lattice of sites with masses . The governing equations of motion are
(1) 
where is spring constant (), is displacement of the th site of the lattice from its equilibrium position and .
In the present work, we consider dichotomic disorder whereby assumes two possible values, either as the background mass of the lattice or to represent mass of the defects. Illustrated in Fig. 1, disordered lattices of length were constructed by arrangement of defects through a powerlaw distribution , where was the random number of lattice sites with mass between two successive defects. As a direct consequence of such choice of , the mass of the entire lattice was the parameter with dependent statistical properties. Obviously, its average density can be written as
(2) 
where is the fraction of defects on the lattice. In the range , is finite. On the other hand, for , vanishes in the thermodynamic limit since the average distance between consecutive defects diverges in this regime. In this paper, we mainly focus on the case (in which is finite) and will comment only briefly on the case which is somehow more peculiar.
Before proceeding, we recall that random walks on such class of Lévy structures have been thoroughly studied in a series of recent papers Barkai et al. (2000); Beenakker et al. (2009); Barthelemy et al. (2010); Burioni et al. (2010); Buonsante et al. (2011); Bernabó et al. (2014) as a minimal model that includes quenched disorder and anomalous diffusion. Depending on the value of , fundamentally different regimes of transport are achieved: indeed walkers superdiffuse for . Several distinguished features of the quenched nature of disorder like the importance of initial conditions and the consequences on higherorder statistics of the diffusive process are discussed in the above mentioned papers. It can be thus envisaged that also localization properties may display unusual features.
ii.2 Statistical properties of the disorder
Using the algorithm described in Devroye (1986), positive integer powerlaw random numbers were generated to construct disordered lattices. For a preliminary statistical characterization, we computed the ensembleaveraged power spectrum:
(3) 
where the average is over different realization of the disorder. The lowwavenumber behavior of provides information on the largescale decay law of the disorder correlation function which will be necessary for the subsequent analysis. Fig. 2 indicates that, for , has a powerlaw singularity for small . This implies that the disorder correlation is indeed decaying as an inverse power of the relative distance with an exponent . The exponent will be important in the subsequent analysis: fitting of the numerical data suggests the following relation with the Lévy parameter (see the Inset of Fig. 2)
(4) 
A theoretical justification of the above relation can be given based on the similarity of this process with the correlation of a Lévy walk as discussed by Geisel in Ref. Geisel (1995). Note that in all the examined cases , so the power spectrum is integrable and the associated process can be considered as stationary for any .
For what concerns the dependence on the lattice size (data not shown), by comparing data for different lengths, we found that for the spectra are independent, while for the spectra decrease with increasing roughly like . This is because the density of defects vanishes in the thermodynamic limit and almostordered realizations dominate the statistical averages.
To conclude this subsection, let us mention that a similar model has been studied in Ref. de Moura et al. (2003) (see also the related work Costa and De Moura (2011)). The main differences with our work is that the sequence of masses is generated as a trace of a fractional Brownian motion that, by construction, has also a a powerlaw singularity at small wavenumber. It was there shown that in the nonstationary case ( in our notation), a mobility edge can exist at a finite frequency value.
ii.3 Transfer matrix approach
As it is well known, the localization properties can be studied by the transfer matrix method Crisanti et al. (1993). Assuming that oscillates harmonically in time at an angular frequency , so that , Eq. (1) can be written as an eigenvalue problem
(5) 
that, upon defining
for can be recast in the form
(6) 
where is a transfer matrix of the th site on the lattice. By iteration of the transfer matrix over the entire lattice and applying the appropriate boundary conditions of the system, solutions of Eq. (5) as a function of frequency can be obtained.
The central quantity to be computed is the Lyapunov exponent which gives the inverse of the localization length . It is convenient to define a new parameter . Inserting into Eq. (5), one obtains the following recursive relation Derrida and Gardner (1984)
(7) 
Eq. (7) can be interpreted as a ”discrete time” stochastic equation. The mass plays the role of a noise source (with bias) whose strength is gauged by the frequency . In the present model, Lyapunov exponent as the inverse of localization length can be computed by
(8) 
Also the integrated density of states follows from node counting arguments, i.e. , where is the fraction of negative values.
For deeper analysis of the systems we are interested in, let us rewrite Eq. (5) in a different notation, introducing the new variables and , the transfer map can be reformulated as follows. Letting , and is a zeroaverage random variable and is finite at least for . Plugging and into Eq. (5) and setting , one finds the following twodimensional map Izrailev et al. (1998)
(9)  
(10) 
The following transformation relations can be introduced by using the canonical variables where and are amplitude and phase of the eigenvector at site , respectively:
(11) 
Substituting Eqs. (11) back into Eqs. (9) and (10) and neglecting terms of the order and higher, one then obtains the following map
(12)  
(13) 
Eq. (13) describes evolution of the phase as it is perturbed by disorder and can be approximated as Izrailev et al. (2012)
(14) 
Longrange correlations of lead to subdiffusive phase fluctuations. Within the same approximation, the Lyapunov exponent is obtained by expanding the logarithm
(15) 
The main issue is to estimate the correlators in Eq. (15). Following the same path as in ref. Izrailev et al. (2012) (see in particular section (5.2)), it can be shown that
(16)  
(17) 
We should get something like
(18) 
Note that this formula agrees with the longwavelength limit of Eq.(22) in HerreraGonzález et al. (2010) which was obtained for weak disorder (where ). The sum is basically the Fourier transform of the correlation function, i.e. the power spectrum of the sequence of masses studied above. We thus expect three cases:

For shortrange correlations as in the case , the second term in Eq. (18) is finite and independent of the frequency. So it just provides a proportionality constant. The Lyapunov exponent is thus expected to scale, up to some constant prefactor as in the standard uncorrelated case Matsuda and Ishii (1970),
(19) 
For longrange correlations and in the case , the fluctuation spectrum diverges at small frequencies so the second term in Eq. (18) dominates so that
(20) Using the result Eq. (4) that the spectrum goes as and recalling the definition of this means that the Lyapunov exponent should be proportional to .

Furthermore, if we assume that the same argument holds also for , we may conclude that , which means that vanishes for . This argument may not be entirely correct since in this case, also vanishes.
To conclude this section, we mention that John and Stephen John and Stephen (1983) studied a related model of a classical wave in the continuum limit where the mass fluctuations were quenched random variables but powerlaw correlated in space . According to their Eq. (7.5), in 1D the Lyapunov exponent should be
(21) 
Although the models are pretty similar at large scales, our result is different from this last estimate. The origin of the discrepancy is not clear at this stage, however it is true that the approach in ref. John and Stephen (1983) relies on some approximate selfconsistent calculations while the calculations reported here are somehow exact.
Iii Numerical results
In this section, numerical results of the model described in the previous section are presented. We have carefully studied the dependence of the Lyapunov exponent on different characteristic parameters of the systems. Data analysis was performed mainly for systems with , however, some computations were done for systems with to allow comparison of the behaviours in two different transport regimes. Throughout the rest of the paper, the spring constant was assumed to be and, unless otherwise stated, the mass ratio was fixed to .
iii.1 Lyapunov exponent
According to the definition described in Eq. (8), the Lyapunov exponent is the exponential growth or decay rate of a vector in the limit .
Due to the longrange spatial correlations (especially in the regime ), the convergence of the numerical value of could vary from a single realization to another and usually requires a long series of iterations. In order to ensure that our obtained results were independent of the disorder realization (i.e. that selfaveraging holds), a series of lattices with different number of defects were studied and compared. Results led us to the conclusion that for , fluctuations of versus are relatively small (better than ) in the systems with different parameters. Clearly, as approaches to , even smaller values of can yield a fast convergence . However, to improve stability of our results even further, defects were used to construct our desired disordered systems. We pursue this number consistently in the rest of this paper for the analysis of very large lattices (asymptotic limit). In all the studied cases, same initial displacement must be assigned to start iteration of the transfer matrix. Note that Lyapunov exponent has this intrinsic property to be independent of and our numerical results were in excellent agreement with that.
Since we mostly focus on the behaviour of the Lyapunov exponent in the small frequency regime, we firstly checked that in all the presented cases the integrated density of states is linearly changing with the frequency, . The physical interpretation of this fact is simply that the spectrum is effectively equivalent to an ordered chain with a renormalized wave speed.
In Fig. 3, variations of the Lyapunov exponent versus frequency for systems with different parameter are shown in doubly logarithmic scale. In each system, as the frequency was increased from to , greater values were obtained for Lyapunov exponent and accordingly, localization length was decreased. As illustrated in the inset of Fig. 3, the scaling relation predicted in the previous section is in excellent agreement with the data obtained by powerlaw fitting on the smallfrequency region of the curves. In some cases, we also checked that the fit that includes the leading order correction follows the data in the whole range, giving further support to the theoretical arguments above.
Another interesting physical effect that should be assessed is the dependence of the Lyapunov exponent on disorder strength or equivalently mass ratio in our model. In an intuitive picture, one expects that higher mass ratio leads to stronger disorder, shorter localization length and thereby greater Lyapunov exponents. In addition to that, the scaling of intrinsically depends on transport properties of the system and hence should be independent of the mass ratio. In other words, the mass ratio should at most change the prefactor in the above scaling relation. To confirm this prediction, four different mass ratios were considered in a system with and length composed of 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 yield an increase in the Lyapunov exponent by orders of magnitude, indicating that the prefactor changes rapidly with .
Fig. 5 shows and compares dependence of the Lyapunov exponent on the Lévy exponent at two fixed frequencies and . Upon decreasing , was initially increased. This is in agreement with the intuition that the localization is enhanced by increasing fluctuations of the distances between defects. However, by approaching , started decreasing and a peak was observed in both curves. We can therefore conclude that in the range , a transition occurs in the behaviour of the Lyapunov exponent. The value of , at which the transition happens, depends slightly on the frequency . According to our numerical data, the transition appeared at for and at for (see the inset in Fig. 5, showing the Lyapunov exponents normalized by their maximums). Moreover, no appreciable dependence on the lattice length is observed in this regime. Although it seems plausible that vanishes below , the origin of the peak is not clear. A possibility is that sampletosample fluctuations become so important that averaging over a much larger ensemble of lattices would be required.
iii.2 Eigenfunctions
Moreover, it is interesting to look at eigenfunctions of a Lévytype disordered system. In general, spatial part of the frequency dependent solutions of Eq. (1) can be obtained by solving
(22) 
where is the diagonalized mass matrix diag and is the matrix of spring constants expressed by
Stationary eigenfunctions of disordered systems can be obtained by numerically solving Eq. (22). For a system with and length consisted of defects, eigenfunctions at three different frequencies are presented in Fig. 6. Top panel shows distribution of mass on the entire system which is a stepfunction of or . Eigenfunctions of the systems at three different frequencies are depicted in the other panels. According to our results, at and , eigenfunctions are almost extended in a long space between two successive defects. At , however, eigenfunctions are localized in a small region inside the system. Note that because of the computational difficulties to diagonalize the system matrix, for this analysis we could only study finite systems with limited number of defects.
iii.3 Phase dynamics
Next we present our numerical results for the spatial evolution of the phase that, using Eq. (11) and the definition , can be computed during the iteration of the transfer matrix by the following relation
(23) 
The computed phase has to be unwrapped to the real axis. Top panel in Fig. 7 illustrates variations of in systems of equal length and different Lévy exponents at .
For the statistical analysis, the average drift rate of the phase () was numerically computed and removed by defining . In order to study evolution of over the entire system, the lattice was divided into sublattices of length . For each sublattice, the quantity , where was computed. Averaged over all the sublattices, evolution of the rmsphase was consequently obtained as shown in the bottom panel of Fig. 7. According to the data, growth of is faster in systems with smaller values of . This can be explained by the fact that disorder fluctuations become stronger which is in agreement with the superdiffusive behaviour of the corresponding randomwalk problem Burioni et al. (2010).
Iv Evolution of wavepackets
In this section, we investigate the consequences of the above results on the spreading of an initially localized perturbation on an infinite lattice. Although every single eigenstate is exponentially localized, the wave propagation can be non trivial as lowfrequency states form a continuum with arbitrarily large localization lengths.
The starting point of this analysis is a Hamiltonian formulation of the system:
(24) 
where and are displacement and momentum of the mass , respectively. From the Hamiltonian, time evolution of any excitation of interest can be derived. Evidently, local energy at each site of the lattice at time can be written as:
(25)  
(26)  
(27) 
where and represent kinetic and potential energies, respectively. Together with the equation of motion Eq. (5), these relations constitute the set of governing equations of onedimensional linear discrete systems.
In our analysis, disordered lattices of length with an excitation at the center were studied. Two different types of excitation were applied: (i) displacement excitation with and ; (ii) momentum excitation with and . As it is known, these classes of initial condition yield different asymptotic properties due to the fact that the mode energy distribution is different in the two cases Datta and Kundu (1995). Based on the conservative dynamics of the system and by using a fourthorder symplectic algorithm McLachlan and Atela (1992), the set of coupled governing equations were numerically solved. We mainly focused on lattices with and to enable comparison of the two different regimes of transport. Note that the number of defects on the considered lattices were different as was kept fixed. Obtained results were then averaged over realizations. Fig. 8 presents spread of the averaged energy density over the two systems after for both momentum and displacement excitations. In addition to a faster spread of energy, a broadening appeared around the excitation for the system with .
Following the analytical approach of ref. Lepri et al. (2010) and assuming the frequencydependent localization length as in the smallfrequency regime (), the time and disorder averaged energy in a Lévytype disordered system for a displacement excitation can be written as
(28) 
which asymptotically decays as . For a momentum excitation, can be similarly described by
(29) 
with asymptotic behaviour as . The insets in Fig. 8 show the excellent agreement between our numerical results and the analytical predictions of the dependent powerlaw decay. In the case of momentum excitation, data for are fitted by with a goodness of as shown in the inset in Fig. 8(a). For the displacement excitation, data for are fitted by with a goodness of which is illustrated in Fig. 8(b).
Moreover, another interesting entity is time evolution of the moments of the Hamiltonian which can be defined as Datta and Kundu (1995); Lepri et al. (2010)
(30) 
where is a positive number and represents the th moment of the Hamiltonian at time . The second moment quantifies degree of the spreading of the wavepacket.
In Fig. 9, time evolution of different moments of a wavepacket are shown and compared for the two systems. Beyond the clear signature of the growth of the disorder averaged moments over time, it was observed that around , the growth exponent suddenly changed in the system with . At shorter times, however, we may assume that correlation effects are not yet strong and system behaved like a shortrange correlated one. Note that such effect is even more emphasized in the case of displacement excitation. To numerically estimate the scaling of , a powerlaw fit was performed on the averaged moments. The obtained results shown in Fig. 10, reveal that is smaller in the system with in comparison with .
As already shown and discussed, is different in systems with different Lévy parameters. In an effort to determine a general analytic relation for that describes dependence on both parameters and , the asymptotic forms of Eqs. (28) and (29) were inserted into Eq. (30). Later, by defining an upper cutoff for the summation in the numerator at the ballistic distance , it can be shown that for displacement excitation
(31) 
while for momentum excitation
(32) 
As shown in Fig. 10, there is an excellent agreement between the theoretical estimates, (according to Eqs. (31) and (32)), and the numerical data obtained by powerlaw fits on in the longtime regimes.
V Summary and conclusion
In this paper, we have studied the localization properties of classical lattice waves in the presence of powerlaw correlated disorder. The model is mathematically simple and the choice of the disorder is motivated by the samples known as Lévy glasses Barthelemy et al. (2008). In the first part, we computed frequencydependent localization length and analyzed how different characteristic parameters can affect it. Using theoretical arguments and numerics, we have shown that in the regime , so that mostly decreases as increases and this trend remains the same even in the regime . The physical interpretation of this result is that, for small frequencies and long wavelengths, the waves see an effective disorder whose fluctuation is scaledependent. This stems from the fact that the disorder is, by construction, scalefree.
Large ordered regions can exist inside the system depending on the Lévy exponent . As expected, we showed that eigenfunctions of the system are localized where the density of disorder is high and extended in the large dilute spaces inside the system. Moreover, due to the superdiffusive nature of the transport in the range , we showed that the spatial growth of the rootmeansquared phase is greater in systems with smaller values of .
The anomalous localization properties reflects in the problem of wavepacket spreading. An initially localized perturbation attains at large times a characteristic powerlaw decay as in the uncorrelated case Lepri et al. (2010); Krapivsky and Luck (2011). However in the present case, the exponent of the decay depends on the exponent of the disorder distribution. This implies that the disorder averaged momenta must diverge with time, although the initial local energy excitation does not spread completely. The predictions are clearly supported by the numerical results. As a matter of fact, it should be remarked that consideration of alone is not sufficient to conclude that the energy diffuses. Clearly, the origin of such a growth is drastically different to a genuinely (possibly anomalous) diffusive processes observed in nonlinear oscillator chains Zavt et al. (1993); Cipriani et al. (2005); Delfini et al. (2007); Mulansky and Pikovsky (2010) and just stems from the slowly decaying tails.
In the present work we mostly focused on the case . As argued in section II.3, the case should probably correspond to a vanishing Lyapunov exponent. This may indicate that the localization in this regime is not standard. This is further supported by the similarity of the problem with the perturbation growth in strongly intermittent dynamical systems which indeed yields a vanishing Lyapunov exponent and a subexponential growth of perturbations Gaspard and Wang (1988). Moreover, the divergence of the variance of the distances should lead to huge sampletosample fluctuations. Some preliminary data confirm this expectation suggesting that this case deserves a more detailed future study.
Finally, a natural question regards the extension to higher dimensions. Provided one is able to extend the definition of disorder, we expect the problem to be considerably more difficult in view of the nature of the predicted essential singularity of the localization length at vanishing frequencies John and Stephen (1983).
Acknowledgements.
Authors wish to thank Mario Mulansky for reading the manuscript and the financial support from the European Research Council (FP7/20072013), ERC Grant No. 291349.References
 Barthelemy et al. (2008) P. Barthelemy, J. Bertolotti, and D. S. Wiersma, Nature 453, 495 (2008).
 Bertolotti et al. (2010) J. Bertolotti, K. Vynck, L. Pattelli, P. Barthelemy, S. Lepri, and D. S. Wiersma, Adv. Funct. Mater. 20, 965 (2010).
 Anderson (1958) P. W. Anderson, Phys. Rev. 109, 1492 (1958).
 Sheng (2006) P. Sheng, Introduction to wave scattering, localization and mesoscopic phenomena, vol. 88 (Springer Verlag, 2006).
 Falceto and Gopar (2010) F. Falceto and V. A. Gopar, EPL 92 (2010).
 FernandezMarin et al. (2012) A. A. FernandezMarin, J. A. MendezBermudez, and V. A. Gopar, Phys. Rev. A 85, 035803 (2012).
 Wells et al. (2008) P. J. Wells, J. d’Albuquerque Castro, and S. de Queiroz, Phys. Rev. B 78, 35102 (2008).
 Iomin (2009) A. Iomin, Phys. Rev. E 79, 062102 (2009).
 Matsuda and Ishii (1970) H. Matsuda and K. Ishii, Prog. Theor. Phys. Suppl. 45, 56 (1970).
 Barkai et al. (2000) E. Barkai, V. Fleurov, and J. Klafter, Phys. Rev. E 61, 1164 (2000).
 Beenakker et al. (2009) C. Beenakker, C. Groth, and A. Akhmerov, Phys. Rev. B 79, 024204 (2009).
 Barthelemy et al. (2010) P. Barthelemy, J. Bertolotti, K. Vynck, S. Lepri, and D. S. Wiersma, Phys. Rev. E 82, 011101 (2010).
 Burioni et al. (2010) R. Burioni, L. Caniparoli, and A. Vezzani, Phys. Rev. E 81, 060101 (2010).
 Buonsante et al. (2011) P. Buonsante, R. Burioni, and A. Vezzani, Phys. Rev. E 84, 021105 (2011).
 Bernabó et al. (2014) P. Bernabó, R. Burioni, S. Lepri, and A. Vezzani, Chaos, Solitons & Fractals 67, 11 (2014).
 Devroye (1986) L. Devroye, Nonuniform random variate generation (SpringerVerlag, New York, 1986).
 Geisel (1995) T. Geisel, in Lévy flights and related topics in physics (Springer, 1995), pp. 151–173.
 de Moura et al. (2003) F. A. B. F. de Moura, M. CoutinhoFilho, E. Raposo, and M. Lyra, Phys. Rev. B 68, 012202 (2003).
 Costa and De Moura (2011) A. Costa and F. De Moura, J. Phys.: Condens. Matter 23, 065101 (2011).
 Crisanti et al. (1993) A. Crisanti, G. Paladin, and A. Vulpiani, Products of random matrices (Springer, 1993).
 Derrida and Gardner (1984) B. Derrida and E. Gardner, J. Phys. France 45, 1283 (1984).
 Izrailev et al. (1998) F. Izrailev, S. Ruffo, and L. Tessieri, J. Phys. A: Math. Gen. 31, 5263 (1998).
 Izrailev et al. (2012) F. Izrailev, A. Krokhin, and N. Makarov, Phys. Rep. 512, 125 (2012).
 HerreraGonzález et al. (2010) I. F. HerreraGonzález, F. M. Izrailev, and L. Tessieri, EPL (Europhysics Letters) 90, 14001 (2010).
 John and Stephen (1983) S. John and M. J. Stephen, Phys. Rev. B 28, 6358 (1983).
 Datta and Kundu (1995) P. K. Datta and K. Kundu, Phys. Rev. B 51, 6287 (1995).
 McLachlan and Atela (1992) R. I. McLachlan and P. Atela, Nonlinearity 5, 541 (1992).
 Lepri et al. (2010) S. Lepri, R. Schilling, and S. Aubry, Phys. Rev. E 82, 056602 (2010).
 Krapivsky and Luck (2011) P. Krapivsky and J. Luck, J. Stat. Mech: Theory Exp. p. P02031 (2011).
 Zavt et al. (1993) G. Zavt, M. Wagner, and A. Lütze, Phys. Rev. E 47, 4108 (1993).
 Cipriani et al. (2005) P. Cipriani, S. Denisov, and A. Politi, Phys. Rev. Lett. 94, 244301 (2005).
 Delfini et al. (2007) L. Delfini, S. Denisov, S. Lepri, R. Livi, P. K. Mohanty, and A. Politi, Eur. Phys. J. Spec. Top 146, 21 (2007).
 Mulansky and Pikovsky (2010) M. Mulansky and A. Pikovsky, EPL (Europhysics Letters) 90, 10015 (2010).
 Gaspard and Wang (1988) P. Gaspard and X. J. Wang, Proceedings of the National Academy of Sciences 85, 4591 (1988).