# Time evolution of wave-packets in quasi-1D disordered media

## 1 Introduction

Band random matrices (BRM) represent an effective model for both 1D disordered systems with long-range hopping and quasi-1D wires. The bandwidth $`b`$ plays the role of the range of the interaction in the first case, the one of the square root of the number of independent conduction channels in the second. Up to now, studies have been mostly devoted to the analysis of the stationary solutions of the Schrödinger equation and to the corresponding spectral properties of BRM’s . Much less is known about the solutions of the time dependent Schrödinger equation, a topic on which only a few studies have been performed . The partial analogy of this latter problem with the ‘dynamical localization’ phenomenon in the kicked rotor  suggests that an initial delta-like packet spreads diffusively and eventually saturates to a localized state. The width of this asymptotic packet for BRM’s is of the order of $`b^2`$ lattice sites, i.e. the same order as the localization length of all the eigenfunctions .

The theoretically predicted scaling laws for the mean square displacement $`\stackrel{~}{M}`$ were tested numerically and a comparison of the asymptotic form of the wave-packet with a theoretical formula  derived for the 1D Anderson model was attempted . More recently some new theoretical results appeared which give a formula for the time asymptotic packet in the BRM model in the large $`b`$ limit . Therefore, it became important to check numerically this formula and to both investigate how the packet reaches its time asymptotic shape and measure the size of finite $`b`$ corrections. For what the time evolution is concerned some phenomenological expressions were suggested in ref. , based on a power-law convergence of the mean square displacement to its steady state value. However, the presence of large statistical fluctuations prevented the authors of ref.  from assessing whether the time asymptotic scaling is ruled by power law corrections or by the logarithmic corrections to the $`t^1`$ dependence suggested by rigorous results obtained for the 1D Anderson model . The fluctuations $`\mathrm{\Delta }_{\stackrel{~}{M}}`$ of the width of the asymptotic wave-packet with the realization of the disorder constitute an even more controversial issue, since not even the scaling behaviour is clearly understood. Some evidence of an anomalous behaviour was presented in two previous studies of the same problem  and in the kicked rotor . In all cases the numerics was too poor to make a convincing statement about the value of the anomalous exponent.

The bottleneck of the previous simulations was the slowness of the integration scheme, a 4-th order Runge-Kutta with a small time step to obtain a good conservation of probability over a long time span. This low efficiency prevented from reaching sufficiently large values of $`b`$ and from considering a large enough number of realizations of the BRM’s. We have instead implemented a 2-nd order Cayley algorithm, which, being unitary, exactly conserves probability, although the one-step integration error is larger than the one of the Runge-Kutta scheme (a situation similar to those of symplectic algorithms in classical Hamiltonian dynamics). This has allowed us to more than double the maximum bandwidth (from $`b=12`$ to $`b=30`$) and to increase the statistics by a factor four (in the worst case).

As a result, we have been able to complete an accurate analysis of the time evolution of the mean square displacement, finding that there is no need to invoke effective formulas with a power-law time dependence even at relatively short times. We have found a clean evidence of an anomalous scaling of the relative fluctuations of the packet width, which behave as $`\mathrm{\Delta }_{\stackrel{~}{M}}/\stackrel{~}{M}b^\delta `$ with $`\delta =0.75\pm 0.03`$. In order to confirm this anomaly, we have investigated the statistics of the packet width at a specific time in the localization regime. The probability distributions at various $`b`$ values, when appropriately rescaled, superpose, and the resulting universal ($`b`$ independent) curve is definitely different from a Gaussian with an exponential tail at large $`\stackrel{~}{M}`$ values. Finally, we have compared our results with the theoretical formula for the asymptotic wave-packet  finding a convincing agreement. The finite $`b`$ corrections to the $`b\mathrm{}`$ Zhirov expression are of order $`(1/b)`$.

## 2 Model and numerical technique

We have considered the time-dependent Schrödinger equation

$$i\frac{\psi _i}{t}=\underset{j=ib}{\overset{i+b}{}}H_{ij}\psi _j$$
(1)   
where $`\psi _i`$ is the probability amplitude at site $`i`$ and the tight-binding Hamiltonian $`H_{ij}`$ is a real symmetric band random matrix. The band structure of the Hamiltonian is determined by the condition

$$\begin{array}{cc}H_{ij}=0\hfill & \text{if }|ij|>b,\hfill \end{array}$$

the parameter $`b`$ setting the band-width; the matrix elements inside the band are independent Gaussian random variables with

$$\begin{array}{ccc}H_{ij}_d=0\hfill & \text{and}\hfill & \left(H_{ij}\right)^2_d=1+\delta _{ij}\hfill \end{array}$$

where the symbol $`_d`$ stands for the average over different realizations of the disorder. In the present work we have considered the evolution of an ‘electron’ initially localized at the centre (identified with the site $`i=0`$) of an infinite lattice. To this aim we have analysed the solution of equation (1) corresponding to the initial condition

$$\psi _i\left(t=0\right)=\delta _{i0}.$$

Since the wave-packet evolves in a supposedly infinite lattice, it is necessary to avoid any spurious boundary effect due to the inevitably finite size of the vectors used in the numerical computations. This goal has been achieved by resorting to a self-expanding lattice, i.e. a lattice whose size is progressively enlarged according to the development of the wave-function. At each integration step, our program checks the probability that the electron is in the leftmost and rightmost $`b`$ sites, adding $`10b`$ new sites whenever the amplitude $`|\psi _i|`$ is larger than $`\epsilon =10^3`$ in at least one of the $`2b`$ outermost sites. We have separately verified that $`\epsilon `$ is small enough not to significantly affect the computation of the probability distribution. For instance, by lowering $`\epsilon `$ by an order of magnitude, the mean squared displacement (computed over the same disorder realizations) changes only by a few percent. Since this systematic error is not larger than the uncertainty due to statistical fluctuations, it is not convenient to reduce the cut-off as it would turn out in a slower code with a consequent reduction of the statistics.

The Schrödinger equation (1) was integrated by approximating the evolution operator $`\mathrm{exp}\left(iHt\right)`$ with the Cayley form

$$\mathrm{exp}\left(iH\delta t\right)\frac{1iH\delta t/2}{1+iH\delta t/2},$$
(2)   
which implies that the values of the wave-function at two successive time-steps are related by

$$\left(1+\frac{1}{2}iH\delta t\right)\psi (t+\delta t)=\left(1\frac{1}{2}iH\delta t\right)\psi (t).$$
(3)   
Solving the band diagonal system of equations (3) allows one to determine $`\psi (t+\delta t)`$ once $`\psi (t)`$ is known. Cayley’s algorithm is a standard tool for the computation of the solutions of the Schrödinger equation (see for instance ref. ); to the best of our knowledge, this is the first application to the specific field of random Hamiltonians with long-range hopping. Cayley’s form (2) for the evolution operator has two relevant features: it is second-order accurate in time and unitary; in addition, the corresponding integration scheme (3) is stable. Stability is essential in order to study the long time evolution of the wave-packet; as for unitarity, it ensures the conservation of probability and, together with second-order accuracy in time, allows one to choose time steps $`\delta t`$ two or three order of magnitude bigger than those used in Runge-Kutta integration schemes. Indeed, we could make use of a time step $`\delta t10^1`$, to be compared with the time step $`\delta t10^410^3`$ used for the same problem in Refs. . To ascertain how large a $`\delta t`$ could be used, we have compared the solutions obtained through Cayley’s algorithm at various $`\delta t`$ with the exact solution of the Schrödinger equation (1), computed by diagonalizing the Hamiltonian (to avoid boundary effects due to the finite size of the diagonalized matrices, we have considered sufficiently short evolution times). By this way we came to the somewhat surprising conclusion that the validity range of the approximate equality (2) extended up to time steps as big as $`\delta t1/\sqrt{b}`$ (the scaling of $`\delta t`$ with the band-width $`b`$ was necessary to compensate the opposite scaling of the energy eigenvalues with $`\sqrt{b}`$). To check this conclusion, we have computed the mean squared displacement in the localized regime for several values of $`\delta t`$ in the range $`10^2/\sqrt{b}1/\sqrt{b}`$, finding differences of a few percent, not larger than the statistical fluctuations.

This can depend on the fact that the long time evolution of the wave-packet seems to be led by the eigenstates at the band centre. Indeed, in the energy representation, the exact evolution operator and the Cayley form can be written as

$`\mathrm{exp}\left(iH\delta t\right)={\displaystyle \underset{n}{}}|ne^{iE_n\delta t}n|`$   
$`{\displaystyle \frac{1iH\delta t/2}{1+iH\delta t/2}}={\displaystyle \underset{n}{}}|ne^{i\varphi _n(\delta t)}n|,`$   
with

$$\varphi _n(\delta t)=2\mathrm{arctan}\left(E_n\delta t/2\right),$$

where $`|n`$ is the eigenvector corresponding to the energy $`E_n`$. These equations show that, for increasing $`\delta t`$, the approximate equality (2) holds true only in the subspace spanned by the eigenvectors corresponding to the band centre, while at the band edge, where the eigenvalues $`|E_n|`$ tend to $`\sqrt{b}`$, the coefficients $`\mathrm{exp}\left(iE_n\delta t\right)`$ and $`\mathrm{exp}\left(i\varphi _n\delta t\right)`$ become quickly different. Therefore, the eigenstates at the band edges appear to play a minor role in the time evolution. This is probably due to the shorter localization length of such states compared with those in the centre: in fact, for the mean square displacement, all eigenstates are weighted with their localization length .

## 3 Results

To investigate the time evolution of the wave-packet, we have computed the mean square displacement

$$\stackrel{~}{M}(b,t)=u(t)\underset{j=\mathrm{}}{\overset{\mathrm{}}{}}j^2|\psi _j(t)|^2_d.$$
(4)   
Previous studies of this problem strongly suggest that $`\stackrel{~}{M}`$ satisfies the scaling relation

$$\stackrel{~}{M}(b,t)=b^4M(\tau =t/b^{3/2}),$$
(5)   
for large enough values of the bandwidth $`b`$. Nevertheless, in Ref. , where the most detailed numerical investigation has been carried out, it was not possible to obtain a clear verification of the scaling law (5) due to the poor statistics and the small values of $`b`$.

The faster integration algorithm described in the previous section has allowed us both to average over more realizations (400 in the worst case), i.e to reduce statistical fluctuations, and to reach larger values of $`b`$ (namely, $`b=30`$ instead of $`b=12`$ as in Ref. ). The results reported in Fig. 1 for several values of $`b`$ are clean enough to show a convincing convergence from above to a limit shape. In other words, there is no possibility to interpret the deviations as a signature of a different scaling behaviour for $`\stackrel{~}{M}`$.

In order to perform a more quantitative analysis, we have proceeded in the following way: $`M(\tau ,b)`$<sup>1</sup><sup>1</sup>1We have added the variable $`b`$ to underline the residual but asymptotically irrelevant dependence on the band-with. has been averaged over the time interval $`20<\tau <30`$ to obtain the more statistically reliable quantity $`M_t(b)`$. By assuming a dependence of the type

$$M_t(b)=M_{\mathrm{}}(1+ab^\alpha ),$$
(6)   
we have fitted the three parameters $`M_{\mathrm{}}`$, $`a`$ and $`\alpha `$, finding that the convergence rate $`\alpha `$ is very close to 1 (0.95), i.e. that the finite-band corrections are of the order $`1/b`$. The fitted value of $`M_{\mathrm{}}`$ is 0.61. The results for the finite-band correction

$$\delta M=M_{\mathrm{}}M_t(b)$$
(7)   
are plotted in Fig. 2.

The good quality of our numerical data suggests also the possibility to compare the temporal behaviour with the available theoretical formulas. In particular, it has been argued in Ref.  that the existence of the so-called Mott states should imply a $`(\mathrm{ln}t)/t`$ convergence of $`M`$ to its asymptotic value. Therefore, we propose the following expression

$$M(\tau ,b)=M(\mathrm{},b)\left(1\frac{1+A\mathrm{ln}(1+\tau /t_D)}{1+\tau /t_D}\right),$$
(8)   
which is the simplest formula that we have found able to reproduce also the initially linear (i.e. diffusive) regime. For each value of $`b`$, the best fit is so close to the numerical data of Fig.1 to be almost indistinguishable from them (this is why we do not report the fits on the same figure). The meaningfulness of the above expression is further strengthened by the stability of the three free parameters $`M(\mathrm{},b)`$, $`A`$ and $`t_D`$, which allows the calculation of a $`b`$-independent diffusion constant (see below). From the values of $`M(\mathrm{},b)`$, we can extrapolate the asymptotic value $`M(\mathrm{},\mathrm{})`$ in exactly the same way as we have done for $`M_t(b)`$, finding once more a $`(1/b)`$-convergence to a value around 0.70. This result is to be compared with the theoretical prediction $`M(\mathrm{},\mathrm{})0.668`$ . The deviation of about 0.03 can be attributed to the accuracy of the integration algorithm.

Another important parameter that can be extracted from formula (8) is the diffusion coefficient. Indeed, by expanding Eq. (8) for small $`\tau `$, we find

$$\frac{M(\tau ,b)}{\tau }=\frac{M(\mathrm{},b)}{t_D}(1A)=D$$
(9)   
The diffusion constant $`D`$ turns out to be close to 0.50 for all values of $`b`$ and, what is more important, close to the value that we obtain from a quadratic fit of the initial growth rate of the packet. This is a very encouraging result, since it confirms the correctness of formula (8) for both the diffusive and the localized regimes. Let us notice that the value $`D0.50`$ is somewhat smaller than the one reported in  ($`D0.83`$). Taking into account statistical fluctuations and systematic deviations, we find that $`D=0.50\pm 0.05`$.

In past papers, a phenomenological expression involving a power-law convergence to the asymptotic value of the mean squared displacement has been proposed , arguing that it should provide an effective description of both the diffusive and localized regime

$$M(\tau ,b)=M(\mathrm{},b)\left(1\frac{1}{\left(1+\tau /t_D\right)^\beta }\right).$$
(10)   
The success of expression (8) shows that there is no need to introduce anomalous power laws to reproduce the numerical findings. However, for the sake of completeness, we have fitted our numerical data also with Eq. (10), finding an equally good agreement. Therefore, on the basis of the quality of the fit we cannot conclude which of the two expressions is better; nevertheless, it is worth recalling that the former one has the correct asymptotic behaviour and, moreover, the fitted parameters are more stable.

A much more controversial situation exists about the fluctuations of the packet width. Let us introduce the r.m.s. deviation

$$\mathrm{\Delta }_{\stackrel{~}{M}}(b,t)\sqrt{u(t)^2_du(t)_d^2}.$$
(11)   
In fact, it has not yet been clarified how the above variable scales in the large-$`b`$ limit. In particular, the correct value of the scaling exponent $`\nu `$ in the relation

$$\mathrm{\Delta }_{\stackrel{~}{M}}(b,t)=b^\nu \mathrm{\Delta }_M(\tau =t/b^{3/2},b)$$
(12)   
is still unknown; this is why the $`b`$ dependence in $`\mathrm{\Delta }_M`$ is explicitly maintained. A scaling like $`b^4`$ would imply that the packet-width is not a self-averaging quantity, since the relative size of the fluctuations would not go to zero for increasing $`b`$. Conversely, an exponent $`\nu =3`$ corresponds both to self-averaging and ‘normal’ behaviour. In fact, the number $`N_c`$ of independent channels (lattice sites) actively contributing to the localized region (i.e. the localization length) is of the order $`b^2`$. If we assume that all such contributions to the second moment $`\stackrel{~}{M}`$ are independent of one another, then we are led to conclude that the relative fluctuations should decrease as $`1/\sqrt{N_c}=1/b`$, thus yielding an absolute growth as $`b^3`$. Since previous studies  have suggested a small anomaly, i.e. $`\nu `$ slightly larger than 3, we have chosen to report the behaviour of $`\mathrm{\Delta }_M(\tau ,b)`$ for $`\nu =3`$. The data shown in Fig. 3 reveals a drastically different behaviour from what observed in Fig. 1. First of all, the curves tend to grow for increasing $`b`$; moreover, there is no obvious indication of a convergence to some finite value. Altogether, these features imply that $`\nu `$ is strictly larger than 3, qualitatively in agreement with previous simulations.

In order to perform a more quantitative analysis, we have computed the average of $`\mathrm{\Delta }_{\stackrel{~}{M}}`$ and rescaled it to the average of $`\stackrel{~}{M}`$,

$$\mathrm{\Delta }_M^a(b)\frac{\mathrm{\Delta }_{\stackrel{~}{M}}_t}{\stackrel{~}{M}_t}.$$
(13)   
($`_t`$ is again to be interpreted as the average over the time interval $`20<\tau <30`$). The advantage of this renormalization, already adopted in Ref. , is that it reduces finite-band corrections. The results reported in Fig. 4, reveal a clean power law decay with an exponent $`\delta 0.75`$. This value is slightly larger than the one found in the previous studies, but follows from a much cleaner numerics. A more global check of the scaling behaviour can be made by plotting the rescaled fluctuations

$$\mathrm{\Delta }_M^g(\tau )b^\delta \frac{\mathrm{\Delta }_{\stackrel{~}{M}}}{\stackrel{~}{M}}.$$
(14)   
for the various values of $`b`$. The optimal value of $`\delta `$ can thus be identified as that one yielding the best data collapse. The curves reported in the inset of Fig. 4 have been obtained for $`\delta =0.75`$. It is necessary to modify $`\delta `$ by at least $`\pm 0.03`$ units in order to see a significant worsening of the data collapse. Accordingly, the best estimate of the anomalous exponent is $`\delta =0.75\pm 0.03`$, so that the dependence of $`\mathrm{\Delta }_M`$ on $`b`$ in Eq. (12) is removed for $`\nu =4\delta 3.25`$.

In order to find further support for this anomalous behaviour, we have investigated the probability distribution $`P(M)`$ for the second moment at the time $`\tau =30`$ (the longest time we reached for the larger $`b`$-values), i.e. when the wave-function has almost entered the steady-state regime. The construction of reliable histograms has forced us to consider smaller values of $`b`$. In fact, we have studied the cases $`b=8`$, $`b=12`$ and $`b=22`$, using $`10^4`$ realizations of the disorder in the first two cases and $`10^3`$ in the last one. The results are reported in Fig. 5, where, following a method suggested in Ref. , we have conveniently rescaled the probability distribution $`P(M)`$. In particular, defining by $`M_{av}`$ the average value of $`M`$ at $`\tau =30`$ and $`\sigma (M)`$ the standard deviation over the ensemble of disorder realizations, we have plotted $`P^{}(M^{})=\sigma (M)P(M)`$ vs. $`M^{}=(MM_{av})/\sigma `$. After this rescaling, the distributions $`P^{}(M^{})`$, corresponding to the three $`b`$ values, have zero average and unit standard deviation. It is remarkable to notice that all curves nicely overlap indicating a striking scaling behaviour. A further important feature is the deviation from a Gaussian behaviour, especially for large values of $`M^{}`$, where a clear exponential tail is visible. The dotted line just above the three curves (corresponding to the pure exponential $`\mathrm{exp}(M^{})`$) has been added to give an idea of the decay rate which is slightly larger than 1. The results of this analysis are important in two respects: i) the exponential tail ‘explains’ the difficulties encountered in getting rid of statistical fluctuations in the estimate of $`\mathrm{\Delta }_M`$; ii) the deviations from a Gaussian behaviour provide an independent evidence of the anomalous scaling behaviour of the fluctuations. It is interesting to remark that a preliminary quantitative comparison has revealed a striking identity of the probability $`P^{}(M^{})`$ with the distributions found in Ref.  for such diverse quantities as the magnetization in the 2D XY model and the power consumption in a confined turbulent flow. This close correspondence deserves further investigations.

Finally, we want to compare our results with the theoretical predictions for the asymptotic shape of the wave-packet. In  a reasonable agreement was found between the numerical data and the formula obtained by Gogolin for strictly one-dimensional systems . Since a theoretical expression has been derived in the meantime for quasi-1D systems , it is desirable to compare our data also with this expression. In Fig. 6 we present the disorder-averaged probability profiles $`|\psi _j(t)|^2_d=\stackrel{~}{f}(j,t)`$ for large times, rescaled under the assumption 

$$f_s(x)=b^2\stackrel{~}{f}(j,\mathrm{});x=j/b^2,$$
(15)   
and compare them with Zhirov’s theoretical formula, which is denoted by the white line. No appreciable deviation is noticeable except for the extreme part of the tails, where it is reasonable to expect numerical errors due to boundary effects.

The good overlap is partly due to the (unavoidable) choice of logarithmic scales in Fig. 6. However, if we zoom the region around the maximum (with the exception of the zero channel), one can see in Fig. 7 a slow tendency of the various curves to grow towards the theoretical expectation. This is consistent with the behaviour of $`M`$ reported in Fig. 1, which reveals a convergence from above for the mean squared displacement. It is interesting to notice that all such deviations are mainly due to the finiteness of $`b`$ while the lack of asymptoticity in $`t`$ appears to be much less relevant.

Finally, we consider separately the zero channel, i.e. the return probability to the origin $`f_s(0)`$. In Fig. 8 we plot this quantity versus $`\tau `$ for different values of $`b`$. In all cases, a quite fast convergence, as compared with the behaviour of the packet-width, to the asymptotic value is clearly seen. In practice, as soon as $`\tau `$ is about 1, the average value of $`f_s(0)`$ reaches the asymptotic value. It is instructive to compare our numerical findings with the asymptotic (in time and $`b`$) analytic expression $`f_s(0)=6`$. By fitting the dependence of the time average of $`f_s`$ (in the interval $`1<\tau <30`$) on $`b`$ as in Eq. 6), we find that the asymptotic value is about 5.7 and that the convergence rate is $`1/\sqrt{b}`$. The numerical value is in a reasonable agreement with the theoretical one, considering that it is the result of an extrapolation of data already affected by errors of the order of a few percent. The non trivial part of the result is the rate of convergence of this probability, which is definitely slower than the $`1/b`$ behaviour displayed by the second (and other low order) moments. The behaviour of the return probability, however, is in agreement with the theoretical predictions made in , where the finite $`b`$ deviations from the asymptotic steady-state probability distribution were estimated to be of order $`O(1/\sqrt{b})`$ in the $`|x|1/b`$ neighbourhood of the origin.

The $`b`$-dependence of the return probability is further illustrated in Fig. 9, where we plot the deviation $`\delta f=6f_s(0)`$ from the asymptotic value $`lim_b\mathrm{}f_s(0)=6`$ as a function of $`b`$ in bilogarithmic scale. The displayed numerical values were obtained by averaging the return probability both over disorder realizations and the time interval $`20<\tau <30`$; the data were then fitted with two expressions, exhibiting deviations from the asymptotic value $`f_s(0)=6`$ of order $`O(1/b)`$ and $`O(1/b^\alpha )`$ respectively. In the second case, the exponent $`\alpha `$ was used as a fitting parameter and the best fit value was $`\alpha =0.53`$. As can be seen from Fig. 9, the power law with $`O(1/\sqrt{b})`$ corrections fits the data much better than the one with deviations of order $`O(1/b)`$.

## 4 Conclusions and perspectives

We have studied the time evolution of an initial $`\delta `$-like wave-packet in a 1D disordered lattice with long-range hopping. The main results of this paper are the following. We have confirmed with clean numerics the scaling law (5) for the mean square displacement $`\stackrel{~}{M}`$, first proposed and studied in Ref. . This scaling law is valid in the large $`b`$ limit; here we have found that finite $`b`$ corrections are of the order $`1/b`$. We have proposed formula (8) for fitting the time evolution of $`\stackrel{~}{M}`$ towards its steady state value; this formula contains the logarithmic corrections suggested by the existence of Mott states. We confirm the presence of an anomaly in the scaling law of the relative fluctuations $`\mathrm{\Delta }_{\stackrel{~}{M}}/\stackrel{~}{M}`$ of the mean square displacement, finding that they vanish for large $`b`$ as $`b^{0.75}`$. We have linked this anomaly to the presence of non-Gaussian fluctuations of the mean square displacement. In fact, the probability distribution of $`\stackrel{~}{M}`$ displays an exponential tail for large values of $`\stackrel{~}{M}`$. The conveniently rescaled probability strikingly coincides with the distributions obtained in Ref.  for such diverse quantities as the magnetization in the 2D XY model and the power consumption in a confined turbulent flow. The degree of universality of such distribution deserves further investigations. Finally, we have compared the numerical results on the steady state probability profile with the theoretical formula proposed by Zhirov for large $`b`$, finding a good agreement. We have computed for the first time finite $`b`$ corrections, obtaining $`O(1/b)`$ deviations for the moments of the probability profile and $`O(1/\sqrt{b})`$ corrections for the return probability to the origin.

## 5 Acknowledgements

We thank B. Chirikov for having suggested the numerical verification of Zhirov’s formula. We thank O. Zhirov for having exchanged with us his ideas and for his comments on our numerical results. We finally thank F.M. Izrailev, H. Kantz and B. Mehlig for useful discussions.