# Spreading of wave packets in disordered systems with tunable nonlinearity

###### Abstract

We study the spreading of single-site excitations in one-dimensional disordered Klein-Gordon chains with tunable nonlinearity for different values of . We perform extensive numerical simulations where wave packets are evolved a) without and, b) with dephasing in normal mode space. Subdiffusive spreading is observed with the second moment of wave packets growing as . The dependence of the numerically computed exponent on is in very good agreement with our theoretical predictions both for the evolution of the wave packet with and without dephasing (for in the latter case). We discuss evidence of the existence of a regime of strong chaos, and observe destruction of Anderson localization in the packet tails for small values of .

###### pacs:

05.45.-a, 05.60.Cd, 63.20.Pw## I Introduction

The presence of uncorrelated spatial disorder in one-dimensional linear wave equations results in the localization of their normal modes (NMs). This is the well-known phenomenon of Anderson localization A58 () and has been experimentally observed in a variety of systems, including as examples light SBFS07 (); LAPSMCS08 () and matter BJZBHLCSBA08 (); REFFFZMMI08 () waves.

Understanding the effect of nonlinearity on the localization properties of wave packets in disordered systems is a challenging task, which, has attracted the attention of many researchers in recent years. Several numerical studies of wave packet propagation in different models showed that the second moment of norm/energy distributions grows subdiffusively in time following a power law of the form M98 (); KKFA08 (); PS08 (); FKS09 (); GMS09 (); SKKF09 ().

In FKS09 (); SKKF09 () the mechanisms of spreading and localization were studied for the disordered discrete nonlinear Schrödinger equation (DNLS)

(1) |

and the quartic Klein-Gordon lattice (KG) of anharmonic oscillators with nearest neighbor coupling. The exponent was numerically found to be close to and a theoretical explanation of this value was provided.

Applying the same theoretical argumentation to a generalized DNLS (gDNLS) model

(2) |

with being a positive integer, the dependence of on was predicted in FKS09 (). The validity of this estimation has not been analyzed in detail. Mulansky M09 () presented numerical simulations of the gDNLS model for a few integer values of . In VKF09 () numerical simulations of the gDNLS model were performed for non integer values of on rather short time scales, leaving the characteristics of the asymptotic () evolution of wave packets aside and open.

The main scope of the present work is to verify the validity and the generality of the theoretical predictions presented in FKS09 () for one-dimensional disordered nonlinear chains as a function of . In particular, we choose to perform numerical simulations for the generalized KG (gKG) model instead of the gDNLS system. This choice was done for two reasons. Firstly, it allows us to test whether the estimations obtained in FKS09 () hold irrespectively of the presence of a second integral of motion (the norm for the gDNLS model). The second reason is a practical one. From the comparative study of DNLS and KG models (which correspond to ) performed in FKS09 (); SKKF09 () it was observed that the KG model requires less CPU time than the DNLS system in order to be integrated up to the same time with the same precision. Since in our study we are mainly interested in the characteristics of the asymptotic dynamical behavior of wave packets, the gKG model was preferred, as it permits long integrations of large lattice sizes within feasible CPU times.

## Ii The generalized Klein-Gordon model

The Hamiltonian of the gKG model is

(3) |

where is the lattice site index, and are respectively the generalized coordinates and momenta, defines the order of the nonlinearity, denotes the disorder strength and are chosen uniformly from the interval . The case corresponds to the standard KG model, which exhibits a similar dynamical behavior with the DNLS model FKS09 (); SKKF09 (). The equations of motion are and yield

(4) |

Hamiltonian (3) is a conservative system and its total energy is preserved and serves as a control parameter of the nonlinearity strength. Equations (4) become linear by neglecting the nonlinear term or in the limit (i.e. ) where . Then the ansatz reduces them to the linear eigenvalue problem

(5) |

The normalized eigenvectors ( are the NMs of the system with the corresponding eigenvalues . The width of the squared eigenfrequency spectrum is with .

The asymptotic spatial decay of an eigenvector is given by where is the localization length KM93 (). The spatial extend of a NM can be characterized by the average localization volume , with being its second moment D10 (). for , and for . The average spacing of squared eigenfrequencies of NMs within the range of a localization volume is for and for .

The squared frequency shift of a single-site oscillator induced by the nonlinearity is

(6) |

where is the energy of the oscillator (see Appendix A). For in Eq. (6) it follows SKKF09 ().

In our study and we follow the evolution of single site excitations for long times and for . Then we have , and . We excite site by setting , for with . In our simulations we use symplectic integration schemes of the SABA family of integrators SKKF09 (); LR01 () with some particularities particular ().

In our computations the number of lattice sites and the integration time step varied between to and to , in order to exclude finite size effects in the wave packet evolution and allow long integrations up to time units. In all our simulations the relative energy error was kept smaller than .

Following the methodology of FKS09 (); SKKF09 () we order the NMs in space by increasing value of their center-of-norm coordinate . We consider normalized energy density distributions with , where is the amplitude of the th NM and the corresponding squared eigenfrequency. In our analysis we use the second moment (), which quantifies the wave packet’s degree of spreading, the participation number , which measures the number of the strongest excited modes in , and the compactness index , which quantifies the sparseness of a wave packet, since decreases as the wave packet becomes more sparse (see SKKF09 () for more details).

## Iii Different dynamical regimes

Different dynamical regimes have been discussed in recent publications KKFA08 (); PS08 (); FKS09 (); SKKF09 (). In a recent paper one of us presented a coherent interpretation of collected numerical data within a simple framework which uses two averaged parameters of an initial wave packet as essential control parameters for the dynamical evolution: the average norm or energy density in a packet, and its typical size F10 (). According to these results, a wave packet can be selftrapped (see also KKFA08 (); SKKF09 ()) in a regime of strong nonlinearity. This happens when nonlinear frequency shifts are larger than the width of the spectrum of the linear equations. If the nonlinearity is weak enough to avoid selftrapping, the wave packet will spread either in an intermediate regime of strong chaos, followed by an asymptotic regime of weak chaos, or the strong chaos regime is skipped, and spreading starts in the regime of weak chaos and stays there. The outcome depends on the ratio of the nonlinear frequency shift of the wave packet after spreading into the full localization volume and the average spacing . Here we will adapt these arguments to the study of single site excitations. Then the three possible regimes are:

(7) | |||

(8) | |||

(9) |

Wave packets in the strong chaos regime spread faster than in the weak chaos case (see also Sect. IV.1 below).

The location of the three different dynamical regimes in the parameter space of the system’s energy and the order of the nonlinearity is shown in Fig. 1. For the regime of strong chaos is absent (in the sense that if at all, it will coexist with selftrapping). Therefore, if not selftrapped, the wave packet is expected to spread in the asymptotic regime of weak chaos. Following Anderson’s definition of localization A58 (), we measure the fraction of the wave packet energy in a localization volume around the initially excited site. For a localized state this fraction should asymptotically tend to a nonzero constant. We find that in the weak chaos regime (lower inset of Fig. 1) the fraction continuously drops down in time, indicating complete delocalization. Contrary, in the strong nonlinear regime of selftrapping (upper inset of Fig. 1) the fraction appears to tend towards a nonzero constant. These behaviors are also clearly seen in Fig. 2.

For the selftrapped regime is shifted to very large energies. At the same time for the regime of strong chaos is widening its window with decreasing . A wave packet, if not selftrapped, may then spread in the intermediate regime of strong chaos, and cross over to weak chaos at energy densities which are the smaller the smaller is. Therefore the crossover to the asymptotic regime of weak chaos may be pushed to very large times, if the initial energy is fixed, and lowered. For low enough initial energies the regime of strong chaos can be again avoided. However low initial energies imply large time scales which are needed to detect any type of spreading FKS09 (); SKKF09 (). Since for the equations become linear again, Anderson localization is restored for all times. Therefore, we expect that for fixed initial energy and , the characteristic time scales diverge as well.

Representative examples in the selftrapping and weak chaos regimes for and , are shown in Fig. 3. In the regime of weak chaos wave packets initially evolve as in the linear case: i. e. they show Anderson localization up to some time and both and remain constant. For the wave packet starts to grow with , (blue curves). The values for and for , obtained by a theoretical prediction for the weak chaos regime given in FKS09 (); F10 () (see also Sect. IV.1 below), describe quite well the numerical data. Increasing the energy shortens the time (green curves). The evolution of the compactness index for these cases is shown in the insets of Fig. 3. For both values of the compactness index eventually oscillates around some constant non-zero value, implying that the wave packet does not become more sparse when it spreads. In the regime of selftrapping a large part of the wave packet remains localized, and therefore is practically constant, while the rest spreads and the second moment increases as (red curves). The time evolution of the energy fraction contained in a localization volume around the initially excited site for in the cases of immediate subdiffusion and of selftrapping is shown in the insets of Fig. 1.

## Iv Spreading of wave packets

### iv.1 Theoretical predictions

In the regime of weak chaos, only a fraction of modes in the packet resonantly interact and facilitate randomization of phases FKS09 (); SKKF09 (); F10 (). The second moment of the wave packet increases in time and follows a power law of the form . In these cases, the participation number follows the law . The dependence of the exponent on the order of the nonlinearity was predicted to be FKS09 (); F10 ()

(10) |

The derivation is based on the result that the probability of resonance for a given packet mode depends on the energy density in the packet. For (strong chaos) it follows , while for (weak chaos) we get F10 (). The diffusion rate is conjectured to be , which in the case of weak chaos, leads to (10) together with .

In the regime of strong chaos we get and the exponent of is given by FKS09 (); F10 ()

(11) |

Since the energy density in the packet is decreasing with increasing time, the condition for strong chaos will be eventually violated, and the spreading will cross over into the regime of weak chaos F10 (). The crossover duration can be very large, and complicates the fitting analysis of numerical data.

If the phases of normal modes are randomized by an explicit routine during the evolution of the system, then the regime of strong chaos is enforced irrespectively on whether the energy density satisfies the criterion of strong or weak chaos. In that case the spreading should again follow the law (11) but now for all times FKS09 ().

The above predictions are expected to hold if a coherent transfer of energy can be neglected, and only incoherent diffusive transfer is relevant.

In previous studies with and single site excitations, spreading should start in the regime of weak chaos. Corresponding fits of the exponent yield for the KG model FKS09 (); SKKF09 (), in agreement with Eq. (10). Mulansky computed spreading exponents for the gDNLS model with single site excitations and M09 (). Since for strong chaos is avoided, the fitting of the dependence with a single power law is reasonable. The corresponding fitted exponents (), () and () agree well with the predicted weak chaos result from (10). For the initial condition has been launched in the regime of strong chaos. A single power law fit will therefore not be reasonable. Since the outcome is a mixture of first strong and later possibly weak chaos, the fitted exponent should be a number which is located between the two theoretical values and . Indeed, Mulansky reports a number . Veksler et al. VKF09 () considered short time evolutions of single site excitations (up to ) for gDNLS models. While the time window may happen to be too short for conclusive results, the observed increase of fitted exponents with increasing for is possibly also influenced by the transition from weak to strong chaos.

### iv.2 Numerical results

In our study we provide numerical evidence for the validity of both predictions (10) and (11), by performing extensive simulations of the gKG model for various values of and for energies away from the selftrapping regime. We consider not only integer but also noninteger values of .

The used energies are not too small in order to avoid long delays on the onset of spreading, and their values cross the boundary between the weak and strong chaos regimes as decreases (Fig. 1). The transition from weak to strong chaos is not an abrupt one. One should keep in mind that the curves plotted in Fig. 1 should be considered as rough indicators of the borders between different regimes, since they are influenced by the characteristics of each particular disorder realization.

#### iv.2.1 Computational techniques

For several values of the nonlinearity order we compute the evolution of single site excitations up to a large final time for an energy value which excludes selftrapping, and for 20 different disorder realizations. For each case we verify that the time evolution of and indicate that the wave packet is not selftrapped. Although we also compute we chose to analyze and present results only for , because it grows faster than allowing a more accurate determination of the exponent . Typically ranges from for small to for larger for which slower spreading is observed. For each realization the is fitted by a law, the value of the exponent is determined and the average value over the 20 realizations is computed. In practice, we perform a linear fit of the , values, instead of a nonlinear fit of the actual , values. The data used for the fitting lie in the time window whose lower end varies from up to , in order to guarantee that even the smallest window contains enough data points (for at least two orders of magnitude of ) allowing a reliable evaluation of .

If the time evolution of is well approximated by a law the averaged value should not depend on the initial time of the time window. In other words, if eventually becomes a constant function of for large enough values of then we consider that satisfactorily describes the asymptotic () behavior of . If on the other hand, does not tend to become constant as increases, the available numerical data cannot be represented reliably by a fit and no exponent can be determined. Such an example is seen in Fig. 4(a) where we see that for increases monotonically as increases. The horizontal dashed line denotes the theoretically predicted weak chaos exponent from Eq. (10) . At the largest values of the exponent is close to the strong chaos result from Eq. (11) (horizontal dotted line).

For (Fig. 4(b)) the values of , after some initial transient time, seem to saturate to a constant value implying that a law well approximates the evolution of . Actually, eventually attains values close to the theoretically predicted weak chaos exponent , denoted by a horizontal line in Fig. 4(b).

A significant dynamical quantity is the minimum time after which becomes practically constant, since it characterizes the onset of the validity of the asymptotic approximation . The value of is estimated from the numerically evaluated function as the minimum time such that for all the local (numerically estimated) rate guarantees a less than 5% change of at with respect to its current value at . This definition has a degree of arbitrariness but it manages to capture for all our simulations the time after which a fit seems to approximate quite accurately the values of . For example for (Fig. 4(b)), we find , and the corresponding exponent value is very close to the theoretically predicted value from Eq. (10).

The cases of and shown in Fig. 4 are two typical examples where an exponent for can or cannot be defined respectively. The time evolution of for three particular disorder realizations for and is plotted in Fig. 4(c). For the three curves tend to increase according to the theoretical prediction given in Eq. (10) (solid line). On the other hand, for seems to increase with an increasing rate, not showing a tendency to approach any power law, at least up to the final integration time . In this case, deviates from the weak chaos theoretical power law prediction given in Eq. (10) (dashed line). Although a reliable numerical estimation of a constant was not possible from our simulations, the values of for large values of seem to be close to the strong chaos prediction given from Eq. (11) (dotted line). The deviation from the power law predictions for can also be clearly seen by plotting the mean value of over the 20 realizations as a function of time (Fig. 4(d)). Similar results for show that the numerical results are in good agreement with the weak chaos prediction .

#### iv.2.2 Subdiffusive spreading

We performed extensive simulations for 25 different values of the nonlinearity order in the interval . For each , we followed the evolution of single site excitations for 20 different disorder realizations by considering an energy value away from the selftrapping regime, which allows the immediate subdiffusion of the wave packet. For each case we tried to determine the exponent () and the time by the above-described procedure. Our criterion for determining and was satisfied for and . For all other tested values of () a reliable value for was not obtained, since exhibited behaviors similar to the one seen for (Fig. 4). It is possible that for these cases integration up to larger final times might allow to attain its asymptotic power law behavior, and permit the estimation of , but the extremely long CPU times needed for such simulations do not make them easily feasible. As we can see from Fig. 1, all these cases belong to the strong chaos regime. The fact that exactly for the intermediate (and possibly rather long lasting) regime of strong chaos may set in, is a possible explanation for the observed difficulties.

The computed exponents for different values of are plotted in Fig. 5 (filled squares). The values of obtained in M09 () for gDNLS are also plotted (empty circles). These values are very close to our results for and , while for a reliable estimate of was not obtained from our simulations. In Fig. 5 the theoretically predicted law (10) is plotted by a dashed line and all computed exponents are in good agreement with it. The exponents for , and slightly deviate from this law, possibly due to the influence of the strong chaos regime which exists for . The time after which is well approximated by (inset of Fig. 5), increases as approaches zero, implying that integration for longer times might be needed in order to estimate for .

From Fig. 1 we expect that for small values of the range of energies for which we observe immediate subdiffusive spreading in the regime of strong chaos should increase. As an example we consider the case of , for which the exponent was obtained for (Fig. 5). In Fig. 6(a) the evolution of (average value over 20 realizations) is plotted for (green curve), (red curve) and (blue curve). For all energies increases in a similar way following a power law, after some transient initial time interval, which is well approximated by the theoretically predicted function (11) (dashed line). It is evident that the exponent , which for was estimated to be at , does not depend on the energy value. Note that the error bar is too large here to discriminate between the strong chaos and weak chaos predictions.

For very small values of the dynamics should approach the behavior of the linear system (), i. e. the wave packet should be localized. This tendency is seen in Fig. 6(b) where is plotted for a disorder realization with in the cases of , , , (curves from top to bottom). It is evident that the time increases as . In particular, for this time was not reached until the end of the simulation at .

The order of nonlinearity influences not only the spreading rate of wave packets, but also the morphology of their profiles. In Fig. 7 we plot the normalized energy distributions of initial single site excitations, for different values in NM (upper plot) and real (lower plot) space. Starting from the outer, most extended wave packet we plot distributions for (black curves), (magenta curves), (red curves), (blue curves), (green curves) and (brown curves). All wave packets were considered for the same disorder realization but at different times of their evolution when they have the same value of second moment . These times are for , for , for , for , for and for and increase for since the spreading becomes slower for larger . As we have seen in Fig. 6, when wave packets remain localized for very large time intervals before they start to spread. This is why for the second moment becomes at a larger time than in cases with and .

From the results of Fig. 7 we see that for large enough values of (), the distributions on a logarithmic scale have a chapeau-like shape consisting of a highly excited central part and exponential tails having practically the same slope. Contrarily, the distributions for and become more extended having different slopes in the tails.

A characteristic of the NM space distributions in Fig.7 for is that they exhibit very large value fluctuations (up to 10-15 orders of magnitude) in their tails, contrarily to the corresponding distributions in real space. Tail NMs are driven by the core of the wave packet, and may also interact with neighboring tail NMs. The presence of large tail amplitude fluctuations signals that neighboring tail NMs do not interact significantly (otherwise we would expect a tendency towards equipartition). Tail NMs are then excited only by the core; the further away they are, the weaker the excitation. But within a small tail volume, NMs with larger localization length will be more strongly excited than those with smaller localization length, hence the large observed fluctuations, which on a logarithmic scale are of the order of the relative variation of the localization length. Therefore Anderson localization is preserved in the tails of the distributions over very long times (essentially until the given tail volume becomes a part of the core). But the NM space distributions for and exhibit less fluctuations in their tail values with respect to the other distributions in the upper plot of Fig. 7, implying that tail NMs are now interacting with each other on comparatively short time scales and reach a visible level of local equipartition. Therefore we observe for these cases a destruction of Anderson localization even in the tails of the spreading wave packets.

#### iv.2.3 Dephasing

The assumption that all NMs in a wave packet are chaotic leads to a power law increase of whose exponent is given by Eq. (11). In order to check the validity of this prediction we enhanced the wave packet chaoticity by a periodic dephasing of its NMs. Every 100 time units on average 50% of the NMs were randomly chosen and the signs of their momenta were changed. In this way a faster spreading of the wave packet, with respect to its normal evolution, was achieved. We also note that has a similar effect on the shape of wave packets as in the case of normal evolution without dephasing, because energy distributions in NM and in real space for the same values and different values have similar profiles to the ones shown in Fig. 7.

Performing a similar numerical analysis (as in the case of normal wave packet evolution) we computed the exponent of and the time for several values of the nonlinearity order . The obtained values are plotted in Fig. 5 by filled triangles. The numerically computed exponents are in good agreement with the theoretical prediction of Eq. (11) (solid line in Fig. 5). In the case of normal wave packet evolution not all NMs in the packet are chaotic, because the exponents presented by filled squares in Fig. 5 are always smaller that the prediction of Eq. (11). The dephasing procedure increases drastically the chaotic nature of the dynamics since the corresponding exponents are quite close, but somewhat smaller, than the predicted values given by Eq. (11) (strong chaos regime).

From the results of Fig. 5 we see that exponents were determined for a larger value interval of () with respect to the normal evolution case. Only for and we were not able to estimate an exponent from the performed numerical simulations.

In addition, time (inset of Fig. 5) has always larger values with respect to the ones obtained for the normal evolution of wave packets, especially for large values of . This means that although dephasing increases the chaoticity of the wave packet, a considerably large amount of time is needed for the evolution to be characterized by .

## V Summary and discussion

We performed extensive numerical simulations of the evolution of single site excitations in the gKG model (3) for different values of the nonlinearity order . According to the analytical treatment presented in F10 (), in this case a wave packet could either: a) be selftrapped for large enough values of the nonlinearity, i. e. the total energy of the gKG system, or b) spread subdiffusively for smaller values of . Particularly for energy values not in the selftrapping regime, the single site excitation belongs either to the weak or the strong chaos regime F10 ().

In the weak chaos regime the wave packet spreads subdiffusively and its second moment grows as . The expected dependence of on in this case is given by Eq. (10). The detrapping time after which spreading stars, increases as decreases because the system is closer to a linear model where no spreading is observed due to Anderson localization. In order to be able to observe spreading for large time intervals we avoided very small energy values.

If the wave packet is launched in the strong chaos regime its subdiffusive spreading is initially characterized by an exponent given from Eq. (11), which is expected to eventually cross over to its asymptotic () value given by Eq. (10). The time at which this crossover starts (as well as its duration) could become very large, limiting our ability to observe it numerically.

According to the estimations of F10 (), if the single site excitation is not selftrapped then it belongs to the weak chaos regime for . For the existence of the strong chaos regime is possible, but not for very small energy values where the dynamics is again in the weak chaos regime (see Fig. 1).

We performed numerical simulations for various values of and for belonging both to the weak and the strong chaos regimes, having care to avoid very small energy values where we might encounter very large detrapping times (Fig. 1). The curves in Fig. 1 do not define exactly the separation between different regimes, instead they should be considered as indications of the location of border regions between them. The energy values used here cross the boundary between the weak and strong chaos regimes around the interval .

Our results verify the validity of Eq. (10) for the weak chaos regime. In particular, the numerically computed exponents are in very good agreement with the theoretical prediction of Eq. (10) for , while they exhibit an increasing deviation from this prediction for , and respectively (Fig. 5). For grows faster than the weak chaos estimation (10), but up to the final times that we performed computations we were not able to reliably fit its evolution with a power law and compute an exponent . Nevertheless, this problem, as well as the deviation of the computed exponents from the theoretical weak chaos prediction (10) for , clearly indicate that in our simulations the boundary between the weak and strong chaos regimes lies in the interval , being in agreement with the theoretical predictions of Fig. 1.

Regarding the simulations for we note that maybe longer integrations could allow for a reliable estimation of . This is a very hard computational task as it requires large CPU times. Since simulations with values around are located close to the borders of different regimes, it is possible that for different disorder realizations the wave packet’s evolution is a mixture of different dynamical behaviors, not allowing the statistical analysis to clearly determine . In addition, for small values of we observe differences in the wave packet dynamics, which could affect the determination of , since Anderson localization is destroyed also in the tails of the wave packets (Fig 7).

In order to enforce the wave packet to evolve continuously in the strong chaos regime we enhanced its chaoticity by repeatedly performing a dephasing of its NMs. When dephasing was applied, the computed exponents were in good agreement with the strong chaos theoretical prediction of Eq. (11) for almost all tested values of (Fig. 5).

Predictions (10) and (11) were derived in Refs. FKS09 (); F10 () for the gDNLS model (2) and for integer values of . Thus, as a final remark we note that our results also establish the generality of these predictions since they proved to be valid for a different dynamical system (the gKG model (3)) and for noninteger values of .

###### Acknowledgements.

We thank J. D. Bodyfelt, S. Fishman, D. O. Krimer, Y. Krivolapov, T. Lapteva, N. Li, M. Mulansky, A. Ponno and H. Veksler for useful discussions.## Appendix A Frequency shift of nonlinear oscillators

Let us consider a nonlinear oscillator described by the Hamiltonian function

(12) |

with . For a given value of the energy the oscillator’s period is

(13) |

where is the positive root of equation . The change of variable , with

(14) |

and being the positive number satisfying

(15) |

transforms Eq. (13) to

(16) |

For we get , and Eq. (12) corresponds to an harmonic oscillator with period and frequency . In order to estimate the change of the squared frequency , we note that the period of the oscillator for is , at a first order approximation in , with ‘’ denoting derivative with respect to . Differentiating Eq. (15) we find . Using this result, and differentiating Eq. (16) according to the Leibniz integral rule we get

(17) |

Then, frequency is given by

(18) |

Since , the squared frequency shift is

(19) |

## References

- (1) P. W. Anderson, Phys. Rev. 109 1492 (1958).
- (2) T. Schwartz, G. Bartal, S. Fishman and M. Segev, Nature 446 52 (2007)
- (3) Y. Lahini, A. Avidan, F. Pozzi, M. Sorel, R. Morandotti, D. N. Christodoulides and Y. Silberberg, Phys. Rev. Lett. 100 013906 (2008).
- (4) J. Billy, V. Josse, Z. Zuo, A. Bernard, B. Hambrecht, P. Lugan, D. Clement, L. Sanchez-Palencia, P. Bouyer and A. Aspect, Nature 453, 891 (2008)
- (5) G. Roati, C. D’Errico, L. Fallani, M. Fattori, C. Fort, M. Zaccanti, G. Modugno, M. Modugno and M. Inguscio, Nature 453, 895 (2008).
- (6) M. I. Molina, Phys. Rev. B 58 12547 (1998).
- (7) G. Kopidakis, S. Komineas, S. Flach and S. Aubry, Phys. Rev. Lett. 100 084103 (2008).
- (8) A. S. Pikovsky and D. L. Shepelyansky, Phys. Rev. Lett. 100 094101 (2008).
- (9) S. Flach, D. Krimer and Ch. Skokos, Phys. Rev. Lett. 102 024101 (2009); ibid. 209903 (2009).
- (10) I. García-Mata and D. L. Shepelyansky, Phys. Rev. E, 79 026205 (2009).
- (11) Ch. Skokos, D. Krimer, S. Komineas and S. Flach, Phys. Rev. E 79 056211 (2009).
- (12) M. Mulansky, Diploma thesis, Universität Potsdam (2009).
- (13) H. Veksler, Y. Krivolapov and S. Fishman, Phys. Rev. E 80 037201 (2009).
- (14) B. Kramer and A. MacKinnon, Rep. Prog. Phys. 56 1469 (1993).
- (15) D. O. Krimer, private communication.
- (16) J. Laskar and P. Robutel, Cel. Mech. Dyn. Astr. 80 39 (2001)
- (17) For the first derivative , and the second derivative , of the nonlinear term , needed for the implementation of any SABA integrator with the addition of the corrector term C, are continuous functions. In particular, for the integration of the and cases where we can substitute by , we apply the SABAC integration scheme. For all other values of , which yield noninteger powers of , the corrector term becomes cumbersome. For these cases the use of SABA integrator was preferred over SABAC.
- (18) S. Flach, Chem. Phys. in print (2010), doi:10.1016/j.chemphys.2010.02.022, ; arXiv:1001.2673v1 (2010).