Correlations of correlations: Secondary autocorrelations in finite harmonic systems

# Correlations of correlations: Secondary autocorrelations in finite harmonic systems

Dan Plyukhin Department of Computer Science, University of Toronto, Toronto, ON, Canada    Alex V. Plyukhin Department of Mathematics, Saint Anselm College, Manchester, NH, USA
July 29, 2019
###### Abstract

The momentum or velocity autocorrelation function for a tagged oscillator in a finite harmonic system decays like that of an infinite system for short times, but exhibits erratic behavior at longer time scales. We introduce the autocorrelation function of the long-time noisy tail of (“a correlation of the correlation”), which characterizes the distribution of recurrence times. Remarkably, for harmonic systems with same-mass particles this secondary correlation may coincide with the primary correlation (when both functions are normalized) either exactly, or over a significant initial time interval. When the tagged particle is heavier than the rest, the equality does not hold, correlations shows non-random long-time scale pattern, and higher order correlations converge to the lowest normal mode.

###### pacs:
05.40.Ca, 05.20.Gg

## I Introduction

The theme of fluctuations in finite systems of harmonic oscillators emerges naturally in both application and theory. From a theoretical point of view, the study of the stochastic dynamics of a tagged degree of freedom in finite harmonic systems provides a valuable illustration, and often more than that, of the role of the thermodynamic and weak-coupling limits, ergodicity, thermalization, recurrences, synchronization, and other basic concepts in nonequilibrium phenomena Zwanzig (); Mazur (); Cukier (); Kac (); Vig (); Jin (); Bend (). Another relevant area is Langevin dynamics generated by a coupling to a finite harmonic bath(s), and its application to mesoscopic systems and networks; see PS (); Hanggi (); Zhou (); Beims (); Onofrio (); Hasegawa (); Akai ().

Being nonergodic, the capability of harmonic systems to illustrate general phenomena in statistical mechanics might seem doubtful at first glance. By means of a canonical transformation a harmonic system of any size can be transformed into a collection of independent oscillators, or normal modes; since the energies of normal modes are the integrals of motion, a single isolated harmonic system does not equilibrate and is not very interesting from the point of view of statistical mechanics.

A more fruitful approach is to consider an ensemble of harmonic systems, assuming that in the past they were in contact with a larger thermal bath in equilibrium at a given temperature, and that the initial normal modes of the ensemble are distributed according to the canonical distribution. Within this framework, one evaluates statistical averages of dynamical variables over the ensemble of the system’s initial coordinates rather than over time. Such averages show the transition of the ensemble to thermal equilibrium in the limit of a large number of particles, and thus the nonergodic nature of harmonic systems does not explicitly manifest itself, and for most cases is inessential. It should however be stressed that this framework, which is standard for most works on stochastic dynamics of harmonic systems both classical and quantum, assumes a very special type of coupling between the system and the external thermal bath: this coupling justifies the initial conditions for the the system’s degrees of freedom, yet is assumed to be sufficiently weak, or completely turned off, as not to affect the system’s further dynamics.

The inequivalence for harmonic systems of ensemble- and time-averages, together with the almost exclusive exploitation in literature of the former, does not necessarily entail that the latter are inadequate. Rather, we introduce in this paper a new class of time-average correlations (we call these secondary correlations) which characterize recurrences in finite harmonic systems. For systems of same-mass particles, these correlations are shown to be very close, and under certain conditions exactly identical, to the conventional (primary) time correlations defined with ensemble averaging. This implies that for finite nonergodic systems, the use of both ensemble and time averages may give meaningful complementary descriptions, and that correlations with the two types of averaging may be related in some subtle way.

## Ii Secondary correlations

Consider the temporal autocorrelation function of a dynamical variable in a finite system of size - typically, such a function exhibits two distinctive regimes, separated by a crossover time of order , where is the speed of signal propagation in the system. For short times , the variable does not feel the presence of the boundaries, and the correlation function decays in a smooth regular way, following the same laws as for an infinitely large system. On the other hand, for longer times the dynamics of the variable are affected by signals reflected from the boundaries. For long time regimes such as this, rather than decaying smoothly the correlation functions may exhibit erratic, apparently noisy, behavior Zwanzig (); Mazur ().

We illustrate this behavior in Fig. 1 by way of the normalized momentum correlation function for the central particle in a harmonic chain with fixed ends. The Hamiltonian of the system is

 H=12mN∑i=1p2i+mω22N∑i=0(qi−qi+1)2, (1)

which describes linearly coupled particles indexed with terminal particles fixed, with displacement . Assuming is odd, the middle particle indexed

 i0=N+12 (2)

has normalized () momentum correlation function

 Ci0(t)=1⟨p2i0(0)⟩⟨pi0(0)pi0(t)⟩=2N+1N∑∑′j=1cosωjt, (3)

where the prime indicates that the summation is only over odd . In this expression (we outline its derivation in the Appendix), the terms are frequencies of normal modes

 ωj=2ωsinπj2(N+1) (4)

where is the frequency of a single oscillator, and the average is taken over the equilibrium ensemble of initial conditions. For the correlation function is very close to that of an infinite chain, given by the Bessel function

 Ci(t)≈C∞(t)=limi,N→∞Ci(t)=J0(2ωt). (5)

This can be readily justified by approximating the sum (3) with an integral, and recognizing the latter as the well-known integral representation of , see e.g. Lee (); PS ().

More interesting from the perspective of this paper is the regime in which the correlation becomes irregular, see Fig. 1. It can be shown that the function given by (3) belongs to the class of almost periodic functions: any value which the function achieves once is achieved again, infinitely many times. Traditionally, such functions are characterized by the average frequency with which they return to , or by the reciprocal, i.e. the mean recurrence time . For correlations of type (3) with large , the famous result for the recurrence time, first obtained by Kac Kac () (see also Zwanzig (); Mazur (); Vig ()),

 τ(c)∼eNc2 (6)

implies that recurrences of order or larger are exponentially rare. This result resolves, or rather (being derived for a model system) shows the direction of resolution for the paradoxes of irreversibility Zwanzig ().

In this paper we propose to characterize the irregular part of the function in another way, which is more in the spirit of nonequilibrium statistical mechanics than the mathematics of almost periodic functions. Namely, observing that for large the correlation function appears to behave like stationary noise, we are encouraged to characterize it by a new correlation function

 Di(t)=1⟨C2i(τ)⟩τ⟨Ci(τ)Ci(τ+t)⟩τ (7)

defined with the time average

 ⟨...⟩τ=limT→∞1T∫T0(...)dτ. (8)

Since we are only interested in the interval when behaves irregularly, one might prefer to set the lower integration limit in definition (8) to instead of zero. However this would only be an unnecessary complication, as the limit makes the two definitions numerically equivalent (assuming always that the integral from to converges).

We shall refer to , defined by relations (7) and (8), as the secondary correlation function, and call the primary one. We would like to promote the secondary correlation as a meaningful statistical tool for characterizing the distribution of recurrences times in a system of finite size. Such information is not contained in the Kac formula (6) for the average recurrence time , so the two functions and do not duplicate each other but describe recurrences in complementary ways.

## Iii Relation to primary correlations

Since the primary and secondary correlations and characterize recurrences at different levels and are defined using different types of averaging (over ensemble and time, respectively), the existence of any specific relation between them is perhaps a priori unexpected. Yet a simple numerical experiment with Eqs. (1-6) suggests, for the middle atom of a chain with fixed ends, the equality

 Ci0(t)=Di0(t). (9)

Closer scrutiny reveals that the equality is exact and holds for any , such that the secondary correlation completely repeats the structure of the primary one for both regular () and noisy () domains and has the same crossover time . The proof follows immediately from the relation

 ⟨cosωjτcosωj′(τ+t)⟩τ=δjj′2cosωjt (10)

which holds for an arbitrary spectrum of (nonzero) normal mode frequencies and can be verified by direct evaluation (with the help of L’Hospital’s rule). For this may further be reduced to the familiar orthogonality relation for the Fourier basis, and thus can be considered a generalized form of the latter. From (3) and (10) one obtains for the non-normalized secondary correlation

 ⟨Ci0(τ)Ci0(τ+t)⟩τ= (2N+1)2N∑∑′j,k=1⟨cosωjτcosωk(τ+t)⟩τ= 12(2N+1)2N∑∑′j=1cosωjt. (11)

Normalizing this function to unity at by dividing it by

 ⟨C2i0(τ)⟩τ=12(2N+1)2N+12=1N+1, (12)

one obtains the normalized secondary correlation

 Di0(t)=⟨Ci0(τ)Ci0(τ+t)⟩τ⟨C2i0(τ)⟩τ=2N+1N∑∑′j=1cosωjt (13)

which coincides with the primary correlation , Eq.(3).

One may observe that for the above derivation it is essential that the primary correlation takes the form of a superposition of cosines with equal weights, as in Eq. (3). In general this, of course, is not the case. For example, for a chain with fixed ends described by the Hamiltonian (1), the normalized momentum correlation function for particle with arbitrary index has the form (see Appendix)

 Ci(t)=N∑j=1A2ijcosωjt,Aij=√2N+1sinπijN+1. (14)

For the middle particle this is reduced to (3), whereas for the other particles normal modes enter the expression (14) with different amplitudes . As one can immediately verify, the exact equality of primary and secondary correlations does not hold in these cases. An important example when this equality does hold for any particle is a harmonic chain with periodic boundary conditions. In this case the momentum correlation for each particle is a superposition of equally weighted normal modes Mazur (); Zwanzig ()

 Ci(t)=1NN−1∑j=0cosωjt,ωj=2ωsin(πjN), (15)

and repetition of the above derivation leads again to the exact equality for any particle of the system.

So far, even with the above examples of its validity, the equality of primary and secondary correlations may appear as no more than a curious coincidence. However, further numerical exercises reveal that even when equality does not hold exactly, it remains a very good approximation for the initial time interval , see Fig. 2. The duration of this interval, , is found to depend non-monotonically on particle position , and for any be equal or shorter than the crossover time, . Respectively, for both primary and secondary correlations coincide with the primary correlation for the infinite chain,

 Di(t)=Ci(t)=C∞(t)=J0(2ωt),t

The proof of the approximate equality (16) can be carried out as follows. From the expression (14) for and the definition (7) for , and using the relation (10), one gets

 Di(t)=1N∑j=1A4ijN∑j=1A4ijcosωjt, (17)

or taking into account the expression (4) for normal mode frequencies

 Di(t)=1N∑j=1A4ijN∑j=1A4ijcos[2ωtsin(πj2(N+1))]. (18)

Recognizing here the generating function for Bessel functions

 cos(xsinθ)=J0(x)+2∞∑k=1J2k(x)cos(2kθ), (19)

can be written as a superposition of Bessel functions

 Di(t)=J0(2ωt)+∞∑k=1SikJ2k(2ωt), (20)

with coefficients

 Sik=2N∑j=1A4ijN∑j=1A4ijcos(πjkN+1). (21)

A simple analysis of this expression shows that given , the coefficients are nonzero only for five sets of :

 Sik=⎧⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎩2,k=2(N+1)s−4/3,k=2(N+1)s−2i−4/3,k=2(N+1)(s−1)+2i1/3,k=2(N+1)s−4i1/3,k=2(N+1)(s−1)+4i0,otherwise (22)

where Note that this expression is invariant under the transformations , reflecting the symmetry of the left and right sides of the chain. One can observe that for large and not too close to the end or to the middle of the chain the coefficients are nonzero only for large indices . For instance, for the chain with and the particle , coefficients are nonzero only for As a result, for not too large in the expression (20), the dominating contribution comes from the first term , while the corrections given by the sum involve Bessel functions of large orders which are negligibly small for a significant time interval  Stegun ().

The above consideration not only justifies the equality for , but also accounts for a curious non-monotonic dependence of on the tagged particle index , which we noticed empirically in Fig. 2. For example, according to (22), for and particles the minimal indices for which takes nonzero values () are , respectively. Then keeping only the leading and first correction terms in the exact expression (20), one gets

 D20(t) = J0(2ωt)−43J80(2ωt), D30(t) = J0(2ωt)−43J120(2ωt), D40(t) = J0(2ωt)+13J88(2ωt). (23)

One can verify that these approximations describe the initial deviation of from very well indeed (Fig. 3 shows this for particle and ). Since for small arguments decreases with order , it is clear from (23) that the second correction terms for particles involve Bessel functions of smaller orders, and thus become essential at earlier times than for particle .

If one applies a similar analysis to the primary correlations , Eq. (14), one gets a familiar approximate relation for the left side of the chain Vig (); Lee (); PS ()

 Ci(t)=J0(2ωt)−J4i(2ωt). (24)

Here, in contrast to corresponding relations (23) for , the order of the second Bessel function, which describes effects of finite size, increases monotonically (linearly) with particle index , and so does the crossover time .

In order to study the dependence of the characteristic time , during which , on in a more quantitative way, let us consider the function

 δi(t)=Ci(t)−Di(t), (25)

which is zero when the two correlations coincide for and fluctuates at longer times. For a given , let us define somewhat arbitrarily as the time at which reaches its first local minimum or maximum; see Fig. 4(a). Similarly, we can define the crossover time as the moment when the function

 Δi(t)=Ci(t)−C∞(t) (26)

has its first local extremum, recalling that is the correlation in an infinite system. Using these definitions, we record observations of and for in Fig. 4(b), as a function of particle index . Whereas increases linearly as we approach the central particle, coincides with for and linearly decreases for . As we already know, the primary and secondary correlations coincide for the middle particle , so is identically and diverges here. Somewhat unexpectedly, we find that also diverges, i.e. identically, for (and of course the symmetric case ). Therefore it would appear that diverges whenever it changes from increasing to decreasing, or vice versa. Further calculations for different show that in general the exact equality holds for particles with indices

 i0=N+12,i1=N+13,i2=2(N+1)3 (27)

provided of course that these expressions are integers. For there are three such particles (), two for (), and none for .

Let us show that this phenomenon is readily accounted for with Eqs.(20)-(22) for the secondary correlation . First, from inspecting (22) one might observe that for the minimum for which is non-zero is and comes from the set with . This yields the approximation

 Di(t)=J0(2ωt)−43J4i(2ωt), (28)

which we already used for and in (23). It differs from the corresponding approximation (24) for only by the factor in the second term. Then, from (28) and (24), the difference functions defined above by relations (25) and (26) take the form

 δi(t)=13J4i(2ωt),Δi(t)=−J4i(2ωt). (29)

Since these two functions have local extrema at the same time, by definition we have . Furthermore, since the position of the first maximum of the Bessel function increases approximately linear with  Watson (), Eq. (29) explains the equality of the characteristic times and their linear increase for in Fig. 4(b).

As gets larger still, one observes from (22) that a minimal for which is and comes from the set with . In this case for , instead of (28), we have another approximation

 Di(t) = J0(2ωt)+13Jα(2ωt), α = 4(N+1)−8i, (30)

which we already used for in (23). Since the primary correlation is still given by (24), the difference function in this case reads

 δi(t)=−13Jα(2ωt)−J4i(2ωt)≈−13Jα(2ωt). (31)

The position of its first extremum increases approximately linearly with  Watson () and, as follows from (30), linearly decreases with . This explains the behavior of for in Fig. 4(b).

The transition of from a positive to a negative slope (over the domain ) occurs at , for which (the minimal value of the set ) becomes less than or equal to (the minimal value of the set ). Then the equality gives , which is consistent with our empirical findings (27).

The exact equality for given by (27) can be readily verified using the following expression for the primary correlations

 Ci(t)=J0(2ωt)+∞∑k=1TikJ2k(2ωt), (32)

with coefficients

 Tik=2N∑j=1A2ijcos(πjkN+1). (33)

These relations are similar to (20) and (21) for and can be derived in a similar way PS (). For given by (27), one can verify directly from (33) and (21) that for any . Then the comparison of (32) and (20) gives for those values of the exact equality .

## Iv Heavy impurity problem

So far we have discussed finite harmonic systems of similar particles. If a tagged particle is heavier than the rest, it turns out that the equality of primary and secondary correlations, and , does not hold. Though structurally similar - looks like a coarse-grained copy of - the two correlations are quite distinctive on any time scale; see Fig. 5. In particular, the approximation of exponential relaxation for , while good for , is noticeably worse for . Another observation is that for both correlations, being apparently random on a short time scale, show on a larger scale a noisy yet periodically repeating pattern; see the bottom plot in Fig. 5. This feature, absent in systems of equal-mass particles, is made all the more obvious when considering higher order correlation functions , defined recursively as

 Ck+1(t)=⟨Ck(τ)Ck(τ+t)⟩τ⟨C2k(τ)⟩τ, (34)

assuming new notations for and . (In this section we use the notation with a subscript referring to the correlation order, rather than to the index of a particle). For the heavy impurity problem one finds that as the order increases the apparent randomness of correlations on the time scale quickly diminishes, and higher correlations converge to the normal mode with the lowest eigenfrequency :

 Ck(t)→cos(Ω1t); (35)

see Fig. 6. Below we outline a theoretical framework underlying these empirical observations.

Consider a cyclic chain of particles of mass and an impurity of mass described by the Hamiltonian

 H = P22M+2N∑i=1p2i2m+mω222N−1∑i=1(qi−qi+1)2 (36) + mω22[(Q−q1)2+(Q−q2N)2],

where and are the momentum and coordinates of the impurity. Using a diagonalization method similar to that described in the Appendix (see Cukier () for details), one can show that the normalized momentum correlation function for the impurity is again an almost periodic function, now of the form

 C(t)≡C1(t)=2N−1∑j=0,1,3,⋯AjcosΩjt. (37)

The amplitudes in this expression are given by

 Aj={1+2N−1∑i=1,3,⋯ϵ2i(Ω2j−ω2i)2}−1. (38)

where

 ωi = 2ωsiniπ2(2N+1), ϵi = −2μ12ω2(22N+1)12siniπ2N+1, (39)

and is the mass ratio.

Due to the system’s symmetry only the modes with zero and odd indices contribute to the superposition (37). Their frequencies ( for cannot be expressed in closed form and must be evaluated as roots of the secular equation Cukier ()

 G(z)=z2−2μω2−2N−1∑i=1,3,⋯ϵ2iz2−ω2i=0. (40)

This transcendental equation has solutions , . It can be verified that one solution is the zero frequency , which reflects the translational invariance of the system. The remaining nonzero roots lie in the interval and must be evaluated numerically.

With the set of eigenfrequencies found, one may calculate the amplitudes with (38) and evaluate the primary correlation by carrying out the summation in (37). Then, using (10), for the secondary correlation one obtains

 D(t)=C2(t)=c2{A20+122N−1∑j=1,3,⋯A2jcosΩjt} (41)

with normalization coefficient

 c2=(A20+122N−1∑j=1,3,⋯A2j)−1. (42)

Fig. 5 presents and , calculated with Eqs. (37) and (41), for and . In a similar manner, one can obtain the expression for order- correlations from (34)

 Ck(t)=ck{Aαk0+2N−1∑j=1,3,⋯(Aj√2)αkcosΩjt} (43)

with powers and normalization coefficient

 ck={Aαk0+2N−1∑j=1,3,⋯(Aj√2)αk}−1. (44)

While expression (43) is a superposition of modes, one can observe that for larger the main contribution comes from the mode with eigenfrequency , such that . This can be accounted for by noticing that the sequence of coefficients has as its maximum element and is monotonically decreasing for . For example, for and we find that (approximately). For the primary and secondary correlations involving and respectively, such an insignificant difference in values hardly plays a role. But for higher-order correlations the maximum of the set (still at ) may be orders of magnitude greater than any other element. As a result, the superposition in (43) is increasingly dominated by the term with , and higher-order correlations quickly converge to the first normal mode; .

For systems of same-mass particles the set of normal mode amplitudes, given by the second equation in (14), is a periodic function of the mode index and has no single maximum. In this case the reduction of higher-order correlations to a dominating normal mode does not occur.

## V Conclusion

Temporal autocorrelation functions are often evaluated in the thermodynamic limit, in which case they typically decrease in a regular (non-random) fashion, either monotonically or non-monotonically. In finite systems, autocorrelation functions themselves become noisy at long time scales ; this illustrates recurrences in the dynamics of the tagged variable due to reflections of sound off boundaries. In this paper we introduced and studied some properties of the secondary correlation function defined as an autocorrelation function of the primary correlation function . If it exists, the characteristic time of decay for determines the time-scale of a typical “period” for , which in turn may be associated with the typical recurrence time of the targeted variable. These “typical” times may however be ill-defined mathematically (as is indeed the case for the harmonic systems discussed above), so to be more precise the secondary correlation can be described as a function characterizing a distribution of recurrence times: for a given , a larger value for corresponds to a greater probability (density) that will return to an assigned value in time about . Comparing the secondary correlation with the mean recurrence time , Eq. (6), the latter being more prevalent in literature, one notices that the two functions give complementary descriptions: while characterizes the number of returns to an assigned value , the secondary correlation gives the distribution of return times regardless of the assigned value of .

One interesting result is the equality for systems of same-mass particles. The equality is either exact for all or a very good approximation over the initial interval whose duration depends on the tagged particle’s position non-monotonically. Although its derivation is quite simple, the equality of primary and secondary correlations may be a remarkable property, especially considering that the former is defined over the ensemble and the latter with time averaging. We restricted the discussion to the simplest case of one-dimensional harmonic systems, but an extension to higher dimensions appears to be straightforward. Whether the equality, or perhaps some other relation, between primary and secondary correlations still holds for nonlinear systems is an open question.

Like the primary correlation, for long time-scales the secondary correlations also develop noisy tails (see the insets in Fig. 3), which themselves can be characterized by correlations of higher order. In turn, this new tertiary function has the same structure as the secondary and primary functions, exhibiting regular decay over shorter times and fluctuating over longer times. Thus one can construct an infinite hierarchy of higher order correlations whose order-scaling properties are also interesting to study. Of course, for the particular cases when the equality holds exactly for all , e.g. a harmonic chain with periodic boundary conditions, all higher order correlations are identical. We have studied higher correlations in the context of the heavy impurity problem. Here the equality of primary and secondary correlations does not hold, and the primarily correlation displays a non-random idiomatic pattern on long time scales, which becomes even more visible in correlations of higher orders. Indeed the sequence of higher-order correlations converges to the lowest normal mode.

###### Acknowledgements.
We thank S. Shea and V. Dudnik for discussion and the anonymous referee for important suggestions.

## Appendix

In this appendix we outline the derivation of expressions (3) and (14) for the momentum correlation functions of -th particle in a harmonic chain with fixed ends, described by the Hamiltonian (1). Using the normal mode transformation

 qi=1√mN∑j=1AijQj,pi=√mN∑j=1AijPj

with coefficients given by (14) and taking into account the orthogonality relation , the Hamiltonian (1) is diagonalized into the form of uncoupled normal modes

 H=12N∑j=1{P2j+ω2jQ2j}

with frequencies given by (4). The normal modes are governed by the Hamiltonian equations

 ˙Pj=−∂H∂Qj=−ω2jQj,˙Qj=∂H∂Pj=Pj,

and evolve as

 Pj(t) = Pj(0)cosωjt−ωjQj(0)sinωjt, Qj(t) = Qj(0)cosωjt+ω−1jPj(0)sinωjt.

Assuming that initially the system is in equilibrium with canonical distribution function , correlations of the normal modes’ initial values are and . Then

 ⟨Pj′(0)Pj(t)⟩=⟨Pj′(0)Pj(0)⟩cosωjt=δjj′βcosωjt

and the momentum correlation of the -th particle is

 ⟨pi(0)pi(t)⟩ = mN∑j,j′=1AijAij′⟨Pj′(0)Pj(t)⟩ = mβN∑j=1A2ijcosωjt.

Division of this expression by gives the normalized correlation function (14). In the case of the middle particle (assuming is odd), for odd and zero otherwise. In this case, one obtains the normalized correlation function in the form (3). The correlation (15) corresponding to the periodic boundary condition can be derived in a similar way. For the extension to the heavy impurity problem see e.g. Cukier ().

## References

• (1) R. Zwanzig, Nonequilibrium Statistical Mechanics, Oxford University Press, New York (2001), Ch.10.
• (2) P. Mazur and E. Montroll, J. Math. Phys. 1, 70-84 (1960).
• (3) R. I. Cukier and P. Mazur, Physica 53, 157 (1971).
• (4) M. Kac, Am. J. Math. 65, 609 (1943).
• (5) J. O. Vigfusson, Physica A 98, 215 (1979).
• (6) F. Jin, T. Neuhaus, K. Michielsen, S. Miyashita, M. A. Novotny, M. I. Katsnelson, and H. De Raedt, New J. Phys. 15, 033009 (2013).
• (7) V. A. Benderskii and E. I. Kats, JETP 116, 1 (2013).
• (8) J. Florencio and M. H. Lee, Phys. Rev. A 31, 3231 (1985).
• (9) P. Hänggi and G.-L. Ingold, Chaos 15, 026105 (2005).
• (10) A.V. Plyukhin and J. Schofield, Phys. Rev. E 64, 041103 (2001).
• (11) H.-X. Zhou and R. Zwanzig, J. Phys. Chem. A 106, 7562 (2002).
• (12) J. Rosa and M. W. Beims, Phys. Rev. E 78, 031126 (2008).
• (13) Q. Wei, S. T. Smith, and R. Onofrio, Phys. Rev. E 79, 031128 (2009).
• (14) H. Hasegawa, Phys. Rev. E 83, 021104 (2011); 84, 011145 (2011).
• (15) A. Carcaterra and A. Akay, Phys. Rev. E 84, 011121 (2011).
• (16) M. Abramowitz and I. A. Stegun, Handbook of Mathematical Functions, Dover, New York (1972).
• (17) G. N. Watson, A Treatise on the Theory of Bessel Functions, Cambridge University Press, Cambridge (1995), p. 521.
You are adding the first comment!
How to quickly get a good reply:
• Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
• Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
• Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters