Rocking bistable systems: Use and abuse of linear response theory

Europhys. Lett.,58 (3), pp. 342–348 (2002) EUROPHYSICS LETTERS 1 May 2002 Rocking bistable systems: Use and abuse of linear response theory J. Casado-Pascual1,2,J.G ´ omez-Ord´ o˜ nez1,M.Morillo 1and P. H¨ anggi2 1F´ısica Te´orica, Universidad de Sevilla - Apartado de Correos 1065 Sevilla 41080, Spain 2Institut f¨ur Physik, Universit¨at Augsburg - Universit¨atsstraße 1 D-86135 Augsburg, Germany (received 19 November 2001; accepted in final form 8 February 2002) PACS. 05.40.-a – Fluctuation phenomena, random processes, noise, and Brownian motion. PACS. 05.20.-y – Classical statistical mechanics. PACS. 02.50.Ey – Stochastic processes. Abstract. – The response of a nonlinear stochastic system driven by an external sinusoidal time-dependent force is studied by a variety of numerical and analytical approximations. The validity of linear response theory is put to a critical test by comparing its predictions with numerical solutions over an extended parameter regime of driving amplitudes and frequencies. The relevance of the driving frequency for the applicability of linear response theory is explored. The response of dissipative physical systems to small-amplitude external perturbations is usually described with the powerful tools of linear response theory (LRT) [1], as it is generally accepted that the effect of the perturbation can be described in terms of small deviations from the behavior of the unperturbed system. In particular, for long times and for systems which in the absence of driving reach an equilibrium distribution, LRT provides an approximate expression for the probability distribution obtained by keeping just the linear terms in a series expansion in the external amplitude. The purpose of the present letter is to point out the relevance of parameters other than the amplitude of the driving force, for the validity of LRT. We will show that for a periodic external force, the validity of LRT depends not just on the amplitude of the driving term but also crucially on its frequency. Let us consider a system characterized by a single degree of freedom, x, whose time evolution is governed by the nonlinear Langevin equation (in dimensionless form) ˙x(t)=x(t)−x3(t)+Acos Ωt+η(t),(1) where Acos Ωtrepresents an external signal and η(t) is a Gaussian white noise with zero average and η(t)η(s)=2Dδ(t−s). The corresponding linear Fokker-Planck equation (FPE) for the probability density P(x, t) reads ∂P ∂t =∂ ∂x−x+x3−Acos ΩtP+D∂2P ∂x2.(2) c EDP Sciences J. Casado-Pascual et al.:Rocking bistable systems: Use and abuse etc. 343 The unperturbed system has an equilibrium distribution of the form Peq(x)=Nexp −U0(x) D,(3) where Nis a normalization constant and U0(x) is the unperturbed potential U0(x)=−x2 2+x4 4.(4) This potential has two minima located at xm=±1 and a maximum at xM= 0, with a barrier height of 0.25. The potential U0(x)−Ax cos Ωtloses its bistable character for A≥AT= 4/27. The analysis of the dynamics is simplified by making use of two important theorems: the Htheorem, which ensures the existence of a unique long-time distribution function P∞(x, t) [2,3], and the Floquet theorem, which guarantees that P∞(x, t) is periodic in time with the same period as the external force [4]. For the system at hand, the symmetry of U0(x) implies the following properties for the long-time unique solution of the FPE: P∞(−x, t;−A)=P∞(x, t;A) and P∞(−x, t;A)=P∞(x, t +T/2; A), where T=2π/Ω and we have indicated explicitly the dependence of P∞on A. Using the Fourier expansion P∞(x, t;A)= ∞  m=−∞ Hm(x;A)eimΩt,(5) the first property leads to Hm(x;A)=Hm(−x;−A), while the second one implies that Hm(x;A)=(−1)mHm(−x;A). From both of them, we obtain Hm(x;−A)=(−1)mHm(x;A). It then follows immediately that the odd moments of the distribution, xn(t)∞,n=1,3,... can be written as Fourier series containing only odd harmonics as the even harmonics vanish due to the symmetries above. Analogously, even moments xp(t)∞,p=0,2,... contain just even harmonics in their Fourier series expansions [5]. Inserting the Fourier expansion, eq. (5), into the FPE, an infinite set of equations for the coefficients Hm(x;A) is obtained. Inspection of the set indicates that if Hm(x, A) is expanded in powers of A, it cannot contain powers smaller than A|m|. From the above general considerations, we have, in particular, for the first two moments, x(t)∞= nodd Mn(A)einΩt=2 ∞  n>0,odd |Mn(A)|cos(nΩt−φn) =2 ∞  n>0,odd |Mn(A)|(cos φncos nΩt+sinφnsin nΩt) (6) with Mn(A)=c(0) nA|n|+c(2) nA|n+2|+···,and x(t)2∞= peven Lp(A)eipΩt(7) with Lp(A)=b(0) pA|p|+b(2) pA|p+2|+···. The exact analytical expression for P∞(x, t) is unknown. LRT amounts to writing P∞(x, t)=Peq(x)+AP(1) 1(x, t),(8) 344 EUROPHYSICS LETTERS with P1(x, t) obtained from a first-order perturbation analysis of the FPE (see [4] and the Appendix of [6] for details(1)): AP(1) 1(x, t)=− ∞  n=1 A λ2 n+Ω 2[λncos Ωt+ΩsinΩt]dnϕn(x),(9) where ϕn(x) are the right eigenstates of the unperturbed FP operator and λnthe corresponding eigenvalues. The coefficients dnare ϕn|∂/∂x|ϕ0. It follows from eqs. (8) and (9) that the average value x(t)LRT ∞is given by x(t)LRT ∞=a1cos Ωt−φLRT 1.(10) The explicit calculation of the amplitude, a1, and phase lag, φLRT 1, requires the knowledge of the spectrum of the unperturbed system. For the bistable system at hand, no exact analytical expressions for the eigenfunctions and eigenvalues exist, although useful approximate expressions are known [7, 8]. Alternatively, the amplitude and phase lag [6, 9, 10] can be obtained from the response function. Using the two-mode approximation of Jung and H¨anggi [9], we write x(t)LRT ∞=b1cos(Ωt−β1)+b2cos(Ωt−β2),(11) where the first term on the right-hand side is due to the interwell hops, while the second one describes the influence of intrawell dynamics. It is convenient to cast the expression for x(t)LRT ∞as in eq. (10), and within the two-mode approximation, we get for the amplitude a1=A Dg2 1λ2 1 λ2 1+Ω 2+g2 2α2 α2+Ω 2+2g1g2λ1α(λ1α+Ω 2) (λ2 1+Ω 2)(α2+Ω 2) 1 2 ,(12) while the phase lag of the response with respect to the input signal, 0 ≤φLRT 1≤π/2, is given by φLRT 1= arctan g1λ1Ω λ2 1+Ω2+g2αΩ α2+Ω2 g1λ2 1 λ2 1+Ω2+g2α2 α2+Ω2 .(13) In the above formulas, λ1is given by [11] λ1≈√2 π1−3 2Dexp[−1/4D],(14) and α= 2. The weights, g1and g2can be obtained from the expressions g2=λ1x2eq λ1−α+x2eq −x4eq λ1−α,(15) g1=x2eq −g2.(16) To leading order in D, we can replace λ1by λK=√2/π exp[−1/4D], g1≈1andg2≈D/α. This is the limit considered in [12]. Linear response theory leads to the following predictions: the first moment x(t)∞should contain a single harmonics with the frequency of the driving force, the output amplitude should (1)Note that the plus signs in (A.23) of ref. [6] should read minus. This in turn yields a minus sign on the right-hand side in (A.28). J. Casado-Pascual et al.:Rocking bistable systems: Use and abuse etc. 345 order of harmonics 0 0.2 0.4 amplitude Ω=10−1 1 3 5 7 9 0 0.2 0.4 0.6 amplitude Ω=10−4 1 3 5 7 9 Fig. 1 0 100000 200000 t -0.5 -0.3 -0.1 0.1 0.3 0.5 < x (t) > Fig. 2 Fig. 1 – Amplitudes of the Fourier components of x(t)∞for noise strength D=0.1 and input amplitude A=0.2 and frequencies Ω = 10−4(upper panel) and Ω = 10−1(lower panel). Fig. 2 – Time evolution of x(t)for D=0.1, A=0.04 and Ω = 10−4as obtained from the numerical solution of the FPE (solid line), the adiabatic approximation (dotted line), the two-mode LRT (dashed line) and the two-mode LRT to leading order in D(dot-dashed line). behave linearly with A. Certainly, for finite values of D, if the amplitude of the driving force is infinitesimally small, the expansion procedure in Ais valid and LRT applies. The point that we want to address here is that for finite small amplitudes, A<A T, the value of Ω has to be taken into account when applying LRT. The upper limit for the values of Afor which LRT remains valid(2) depends as well on the driving frequency. The adiabatic approximation gives a description of the dynamics when Ω is small compared to any other characteristic frequency of the system. In this approach [4], the probability density is assumed to be given by Pad(x, t)=N(t)exp−U0(x)−Ax cos(Ωt) D,(17) where N(t) is the normalization constant. An analysis of the corrections to the bare adiabatic approximation has recently been presented by Talkner [14]. Even in the absence of driving, no exact explicit time-dependent analytical solution of the FPE for the model system at hand is known. We have resorted to numerical solutions of eq. (2). We follow a technique based on the use of the split propagator method of Feit et al. [15] and detailed in [16]. From the numerical solution of the FPE we can easily obtain the time dependence of x(t)∞. As this is a periodic function of time, its Fourier components can be obtained by numerical quadrature. In fig. 1, we show the amplitudes of the relevant Fourier components of the output signal for D=0.1, A=0.2 and two very different driving frequencies, Ω = 10−1andΩ=10 −4. In this figure, as well as in the subsequent ones, we have taken D=0.1. This is a typical (2)In the linear response regime, the dimensionless ratio A/D is assumed to obey A<D.Intheopposite singular limit, AD, the dynamics assumes universal weak noise spectral properties [6, 13]. 346 EUROPHYSICS LETTERS 0 200 400 t -0.4 -0.2 0 0.2 0.4 < x(t) > Fig. 3 400 425 450 t -0.4 -0.2 0 0.2 0.4 < x(t) > Fig. 4 Fig. 3 – The same as in fig. 2 but with Ω = 10−1. Fig. 4 – The same as in fig. 2 but with Ω = 1.0. Notice that, due to the large value of the driving frequency, we only plot a few cycles of the output in the asymptotic regime. value and it is adequate for the validity of the two-mode approximation leading to eqs. (12), (13). On the horizontal axis we indicate the order of the harmonics. It is clear that even for this driving amplitude, relatively large in relation to its threshold value, the response of the system at the larger frequency contains essentially the first harmonics. On the other hand, for the small driving frequency, higher-order harmonics are generated. This is an indication of the failure of LRT to describe the dynamics at these low frequencies, while LRT might still be a good description for higher frequencies. In fig. 2, we depict the time evolution of x(t)obtained from the numerical solution of the FPE for D=0.1, A=0.04 and Ω = 10−4. We also show the behaviors obtained using the adiabatic ansatz, eq. (17), LRT within the two-mode approximation, eqs. (12), (14)-(16) and LRT to leading order in D. The input signal is largely amplified at this small frequency. The adiabatic result deviates slightly from the numerical one at the peaks. The deviations from the numerical results are larger with the LRT description. Nonetheless, the two-mode LRT and the adiabatic approximations yield an acceptable description of the dynamics. This is expected within the linear response regime, where AD. In fig. 3, we show the behavior for Ω = 10−1. It is clear that the adiabatic approach yields a signal with a very large amplitude and a large phase shift compared with the numerics. The two-mode LRT still yields a very acceptable behavior. The same qualitative features are observed in fig. 4, where Ω = 1. For this large frequency, the deviations of the two-mode LRT from the numerical result are very small and they cannot be noticed in the plot. In fig. 3, we show the full time evolution including the short transient. In fig. 4, as we consider a large frequency value, we only show a few oscillations in the asymptotic regime, so that the details of a cycle can be distinguished. These last three figures show that for very small A, LRT gives a satisfactory description of the system response, with deviations from the numerics more pronounced as the external frequency assumes smaller values. To test the validity of LRT as the input amplitude is increased, we have carried out an extensive numerical analysis of the system response to input signals of increasing amplitudes and different frequencies. We evaluate the relative error eampl =|Aout −a1|/Aout, between the output amplitude, Aout, provided by the numerics and the one obtained within LRT with J. Casado-Pascual et al.:Rocking bistable systems: Use and abuse etc. 347 0 0.05 0.1 0.15 0.2 A 0 1 2 eampl 0 1 2 3 eampl Fig. 5 10-4 10-3 10-2 10-1 100 Ω 0 0.5 1 1.5 Phase lag A=0.01 A=0.05 A=0.1 A=0.2 Fig. 6 Fig. 5 – Plots of the relative error of the output amplitude, eampl =|Aout −a1|/Aout,vs. input amplitude Aand several values of the driving frequency. In the upper panel, a1is evaluated using the LRT two-mode expressions to leading order in D, while the full two-mode formulas are used in the lower panel; see eqs. (12), (14)-(16) in the main text. The noise strength is D=0.1andthe frequencies are: Ω = 1.0 (circles), Ω = 10−1(plus signs), Ω = 10−3(crosses) and Ω = 10−4(triangles). Fig. 6 – Plot of the phase lag between the average output and the driving force vs. the angular frequency Ω. With the solid line we depict φLRT 1, evaluated using the LRT two-mode expressions to leading order in D(see eq. (13) in the main text). The symbols denote the numerically determined values of the phase lag Ψ (see text), for A=0.01 (circles), A=0.05 (triangles), A=0.1 (diamonds), and A=0.2 (crosses). The noise strength is set at D=0.1. the two-mode approximation, a1in eq. (12), as a function of the input amplitude A. Our findings are shown in fig. 5. The upper panel shows the dependence of eampl on Afor several frequencies, when the LRT is evaluated to leading order in D, while in the lower panel, the full expressions, eqs. (12), (14)-(16) have been used. For relatively high frequencies, Ω = 1.0 (circles), and Ω = 10−1(plus signs), the error remains small and is practically constant, even for input amplitudes which are rather large compared to its threshold value. On the other hand, for small values of Ω, Ω = 10−3(crosses)andΩ=10 −4(triangles), the error increases drastically with the input amplitude. In particular, the explicit relative errors at D=0.1, i.e. (e1,e2,e3,e4), corresponding to the driving frequencies (Ω1=10 −4,Ω 2=10 −3, Ω3=10 −1,Ω 4=1.0), respectively, read for A=0.01: (0.028, 0.028, 0.056, 0.063); for A=0.1: (0.249, 0.249, 0.072, 0.0659); and for A=0.2: (2.539, 0.797, 0.090, 0.074). Thus, the output amplitude predicted by LRT at these small external frequencies is very much in error, even though, for the same external amplitudes and moderate-to-large frequencies, LRT predictions are still adequate. The average output lags behind the input signal with a phase shift between 0 and π/2. The value of the phase shift predicted by LRT, φLRT 1given by eqs. (10), (13), is independent of the driving amplitude, but depends on Dand Ω. It starts at 0 for very small frequencies, then reaches a local maximum, and tends to its limiting value π/2 for very large frequencies. Its behavior for the small-to-moderate frequencies considered here (Ω <1) is depicted with the solid line in fig. 6 for D=0.1. On the other hand, the phase lag of the numerical result, Ψ, depends on D,Aand Ω. The Ψ values plotted have been calculated from the difference between the instant of times within a period, at which the driving signal and the periodic output, x(t)∞, cross signs, i.e., the corresponding phase delay in crossing zero. In fig. 6, we 348 EUROPHYSICS LETTERS plot the values of Ψ for several values of the driving amplitude and frequency. For very small frequencies, the output is almost in phase with the input for all the amplitudes considered. As the frequency increases, deviations between the numerical predictions and LRT results are manifested, being larger for larger driving amplitudes. In conclusion, our analysis clearly indicates the influence of the driving frequency Ω on the validity of the LRT predictions for the amplitude, phase and number of higher harmonics of the response of the system to subthreshold input signals. As the driving frequency assumes sufficiently small values, the output amplitude significantly deviates from its linear behavior predicted by the two-mode approximation LRT, even though the driving amplitude might still be quite small in order to preserve the bistable character of the unperturbed potential. Even for subthreshold inputs, higher-order harmonics might contribute to the system response for small driving frequencies, contrary to the predictions of LRT. Although the global behavior of the phase lag indicated by LRT is qualitatively correct, as expected, its quantitative predictions are not reliable as the input amplitude increases. ∗∗∗ Support by the Direcci´on General de Ense˜nanza Superior of Spain (Project No. PB981120), the Junta de Andaluc´ıa (JC-P, JG-O, MM) and the Deutsche Forschungsgemeinschaft HA1517/13-4 (PH) is gratefully acknowledged. REFERENCES [1] Kubo R.,J. Phys. Soc. Jpn.,12 (1957) 570. [2] Lebowitz J. L. and Bergmann P. G.,Ann. Phys. (Leipzig),1(1957) 1. [3] Risken H.,The Fokker-Planck Equation (Springer-Verlag, Berlin) 1984, p. 135. [4] Jung P. and H¨ anggi P.,Phys. Rev. A,44 (1991) 8032. [5] H¨ anggi P., Jung P., Zerbe C. and Moss F.,J. Stat. Phys.,70 (1993) 25. [6] Gammaitoni L., H¨ anggi P., Jung P. and Marchesoni F.,Rev. Mod. Phys.,70 (1998) 223. [7] H¨ anggi P. and Thomas H.,Phys. Rep.,88 (1982) 207. [8] Dykman M. I., Haken H., Gang Hu, Luchinsky D. G., Mannella R., McClintock P. V. E., Ning C. Z., Stein N. D. and Stocks N. G.,Phys. Lett. A,180 (1993) 332. [9] Jung P. and H¨ anggi P.,Z. Phys. B,90 (1993) 255. [10] G´ omez-Ord´ o˜ nez J. and Morillo M.,Phys. Rev. E,49 (1994) 4919. [11] H¨ anggi P., Talkner P. and Borkovec M.,Rev. Mod. Phys.,62 (1990) 251. [12] Gang Hu, Haken H. and Ning C. Z.,Phys. Lett. A,172 (1992) 21. [13] Shneidman V. A., Jung P. and H¨ anggi P.,Europhys. Lett.,26 (1994) 571. [14] Talkner P.,New J. Phys.,1(1999) 4.1. [15] Feit M. D., Fleck J. A. jr. and Steiger A.,J. Comput. Phys.,47 (1982) 412. [16] G´ omez-Ord´ o˜ nez J. and Morillo M.,Physica A,183 (1992) 490.