Localization in one-dimensional chains with Lévy-type disorder
We study Anderson localization of the classical lattice waves in a chain with mass impurities distributed randomly through a power-law relation with as the distance between two successive impurities and . This model of disorder is long-range 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 small-frequency 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 scale-dependent. Numerical simulations show that an initially localized wavepacket attains, at large times, a characteristic inverse power-law front with an -dependent exponent which can be estimated analytically.
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 power-law distribution , where is the so-called 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évy-type disorder that include transport in quantum wires Falceto and Gopar (2010), photonic heterostructures Fernandez-Marin 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 power-law 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 power-law distribution . This model of disorder is clearly inspired by the peculiar structure of Lévy glasses. As usual, the one-dimensional case allows for a detailed study. In particular, the Lyapunov exponent , the inverse of the frequency-dependent localization length can be computed straightforwardly. Our main result is that, in the small-frequency 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 steady-state 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 power-law exponent. In contrast, for the range , depending on the type of the initial excitation, the power-law 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évy-type disordered lattice
We considered a one-dimensional harmonic disordered lattice of sites with masses . The governing equations of motion are
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 power-law 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
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 higher-order 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 power-law random numbers were generated to construct disordered lattices. For a preliminary statistical characterization, we computed the ensemble-averaged power spectrum:
where the average is over different realization of the disorder. The low-wavenumber behavior of provides information on the large-scale decay law of the disorder correlation function which will be necessary for the subsequent analysis. Fig. 2 indicates that, for , has a power-law 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)
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 almost-ordered 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 power-law singularity at small wavenumber. It was there shown that in the non-stationary 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
that, upon defining
for can be recast in the form
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)
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
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 zero-average random variable and is finite at least for . Plugging and into Eq. (5) and setting , one finds the following two-dimensional map Izrailev et al. (1998)
The following transformation relations can be introduced by using the canonical variables where and are amplitude and phase of the eigenvector at site , respectively:
Long-range correlations of lead to sub-diffusive phase fluctuations. Within the same approximation, the Lyapunov exponent is obtained by expanding the logarithm
We should get something like
Note that this formula agrees with the long-wavelength limit of Eq.(22) in Herrera-Gonzá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 short-range 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),
For long-range correlations and in the case , the fluctuation spectrum diverges at small frequencies so the second term in Eq. (18) dominates so that
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 power-law correlated in space . According to their Eq. (7.5), in 1D the Lyapunov exponent should be
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 self-consistent 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 long-range 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 self-averaging 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 power-law fitting on the small-frequency 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 sample-to-sample fluctuations become so important that averaging over a much larger ensemble of lattices would be required.
Moreover, it is interesting to look at eigenfunctions of a Lévy-type disordered system. In general, spatial part of the frequency dependent solutions of Eq. (1) can be obtained by solving
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 step-function 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
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 sub-lattices of length . For each sub-lattice, the quantity , where was computed. Averaged over all the sub-lattices, evolution of the rms-phase 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 random-walk 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 low-frequency states form a continuum with arbitrarily large localization lengths.
The starting point of this analysis is a Hamiltonian formulation of the system:
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:
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 one-dimensional 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 fourth-order 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 frequency-dependent localization length as in the small-frequency regime (), the time and disorder averaged energy in a Lévy-type disordered system for a displacement excitation can be written as
which asymptotically decays as . For a momentum excitation, can be similarly described by
with asymptotic behaviour as . The insets in Fig. 8 show the excellent agreement between our numerical results and the analytical predictions of the -dependent power-law 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).
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 short-range correlated one. Note that such effect is even more emphasized in the case of displacement excitation. To numerically estimate the scaling of , a power-law 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
while for momentum excitation
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 power-law fits on in the long-time regimes.
V Summary and conclusion
In this paper, we have studied the localization properties of classical lattice waves in the presence of power-law 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 frequency-dependent 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 scale-dependent. This stems from the fact that the disorder is, by construction, scale-free.
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 root-mean-squared 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 power-law 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 sample-to-sample 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/2007-2013), ERC Grant No. 291349.
- 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).
- Fernandez-Marin et al. (2012) A. A. Fernandez-Marin, J. A. Mendez-Bermudez, 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, Non-uniform random variate generation (Springer-Verlag, 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. Coutinho-Filho, 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).
- Herrera-González et al. (2010) I. F. Herrera-Gonzá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).