Comment on “Ratchet universality in the presence of thermal noise”

PHYSICAL REVIEW E 88, 066101 (2013) Comment on “Ratchet universality in the presence of thermal noise” Niurka R. Quintero,1,2,*Renato Alvarez-Nodarse,1,3,†and Jos´ e A. Cuesta4,5,‡ 1Instituto de Matem´ aticas de la Universidad de Sevilla, Universidad de Sevilla, 41012 Sevilla, Spain 2Departamento de F´ ısica Aplicada I, EPS, Universidad de Sevilla, Virgen de ´ Africa 7, 41011 Sevilla, Spain 3Departamento de An´ alisis Matem´ atico, Universidad de Sevilla, Apartado 1160, 41080 Sevilla, Spain 4Grupo Interdisciplinar de Sistemas Complejos, Departamento de Matem´ aticas, Universidad Carlos III de Madrid, Avenida de la Universidad 30, 28911 Legan´ es, Spain 5Instituto de Biocomputaci´ on y F´ ısica de Sistemas Complejos, Universidad de Zaragoza, 50009 Zaragoza, Spain (Received 4 September 2013; published 10 December 2013) A recent paper [P. J. Mart´ ınez and R. Chac´ on, Phys.Rev.E87, 062114 (2013)] presents numerical simulations on a system exhibiting directed ratchet transport of a driven overdamped Brownian particle subjected to a spatially periodic, symmetric potential. The authors claim that their simulations prove the existence of a universal waveform of the external force that optimally enhances directed transport, hence confirming the validity of a previous conjecture put forth by one of them in the limit of vanishing noise intensity. With minor corrections due to noise, the conjecture holds even in the presence of noise, according to the authors. On the basis of their results the authors claim that all previous theories, which predict a different optimal force waveform, are incorrect. In this Comment we provide sufficient numerical evidence showing that there is no such universal force waveform and that the evidence obtained by the authors otherwise is due to their particular choice of parameters. Our simulations also suggest that previous theories correctly predict the shape of the optimal waveform within their validity regime, namely, when the forcing is weak. On the contrary, the aforementioned conjecture does not hold. DOI: 10.1103/PhysRevE.88.066101 PACS number(s): 05.60.Cd, 05.40.−a, 05.70.Ln, 07.10.Cm The authors of Ref. [1] (see also the erratum [2]) simulate the equation ˙ x+sin x=√σξ(t)+γF bihar(t), (1) Fbihar(t)=ηsin(ωt)+2(1 −η)sin(2ωt +φ), where γis the global amplitude of the force; 0 ⩽η⩽1 and φ account for the relative amplitude and initial phase difference of the two harmonics, respectively; ξ(t) is a Gaussian white noise with zero mean and ξ(t)ξ(t+s)=δ(s); and σis proportional to the temperature of the system. This system exhibits ratchet transport if the external force breaks both a time-shift symmetry, namely, if Fbihar(t)=−Fbihar(t+T/2) (Tbeing the period of Fbihar), and a time-reversal symmetry, i.e., Fbihar(t)=−Fbihar(−t). This happens for all 0 <η<1 and all φ= 0,π. If initially the particle starts at x(t0)=x0, the ratchet current can be obtained as v=lim t→∞ x(t)−x0 t−t0 ,(2) where · represents an ensemble average over all trajectories satisfying the same initial condition. Obviously the ratchet current vwill be a function of the parameters of the system, in particular of those that define the external force. Since for η=0,1orφ=0,π the force breaks neither the time-shift symmetry nor the time-reversal symmetry (hence v=0), it is easily foreseen that for a certain combination of the parameters of the force vmust be maximal (in absolute value). *[email protected] †[email protected] ‡[email protected] Based on a conjecture proposed by one of the authors [3], vshould be optimal when the force maximally breaks the symmetries. For σ=0 this happens for η=4/5 irrespective of the values of γand φ(as long as φ= 0,π)[1]. An argument based on an affine transformation of the force leads the authors to conclude that this optimal shape of the force will hold even for σ>0, albeit some deviations are to be expected. This result is universal in the sense that is independent of γand φ. Figure 1(a) of [2] confirms that this is an accurate prediction even for the high intensities of the noise they use in their simulations (σ=2,3,4). The other parameters are set to ω=0.08π,φ=π/2, and γ=2 throughout their paper. They go on to claim that since all previous theories [4–9] predict a form of the ratchet current given by [2] v∝γ3η2(1 −η),(3) they all predict that vis optimal for η=2/3, a value certainly far away from the simulation results. Accordingly, the two main conclusions of this work are that (i) the conjecture of a universal force waveform that optimizes the current is confirmed even in the presence of strong noise, albeit with some deviations, and (ii) all previous theories must be incorrect because they incorrectly predict this form. We have carried out extensive simulations of the same system (1) and with the same parameter as the authors of [1,2], but instead of limiting ourselves to the single value of the global amplitude γ=2 used in their simulations we have covered a wider range of values, from γ=6downtoγ=0.8. Below this value simulations are prohibitively long because the high values of the noise intensity demand a very large number of realizations to achieve reliable results. The outcome of these simulations is summarized in Fig. 1, which represents the value of η(henceforth ηopt) that optimizes vas a function of the global amplitude of the external force γ. 066101-1 1539-3755/2013/88(6)/066101(3) ©2013 American Physical Society COMMENTS PHYSICAL REVIEW E 88, 066101 (2013) 0123456 γ 0.66 0.68 0.70 0.72 0.74 ηopt FIG. 1. (Color online) Values of the parameter ηdefining the relative amplitudes of the two harmonics of the external force Fbihar(t) [see Eq. (1)], for which the ratchet velocity vreaches its maximum absolute value, plotted as a function of the global amplitude of the external force γ(bullets). The remaining parameters are chosen as in Fig. 1of [1] (actually of the Erratum [2]): ω=0.08π,σ=2, and φ=π/2. The red diamond at (0,2/3) represents the theoretical prediction of previous theories [4–9], which should hold in the limit γ→0. There are two main conclusions that we can extract from this figure. First of all, our results are not compatible with the existence of an optimal force waveform. The values of ηopt range from near 0.69 up to near 0.75. The predicted universal value ηopt =4/5 is reached at no value of γand the closest it gets to it is at γ≈2, precisely the value used in the simulations of Refs. [1,2]. We need to make clear at this point that by setting γ=2 in our simulations our results reproduce accurately the plots of Fig. 1(a) of this work. The apparent (approximate) confirmation of the prediction ηopt =4/5 seems to arise from 0.6 0.65 0.7 0.75 0.8 η -0.07 -0.06 -0.05 -0.04 -0.03 v FIG. 2. (Color online) Value of the ratchet current vas a function of the parameter ηdefining the relative amplitudes of the two harmonics of the external force Fbihar(t)[seeEq.(1)]. The parameters are γ=2, ω=0.08π,σ=2, and φ=π/2. Bullets are the results from simulations averaged over 5000 realizations of the noise. The yellow curve is a fourth-order polynomial fit to these results (v= −3.7188 +22.068η−48.018η2+44.984η3−15.366η4and correlation coefficient r=0.99). This fit is used to determine ηopt. Notice that the minimum of the fit (at η=0.742) does not coincide with the value of ηfor which the largest absolute value of voccurs because of the fluctuations in the ratchet current. the specific choice of simulation parameters made by the authors. Second, although we cannot decrease γbelow 0.8 without introducing too much uncertainty, the figure clearly illustrates that the trend of the value of ηthat optimizes vis toward the value 2/3, which all theories predict in their range of validity, i.e., in the limit of weak external forces. On the basis of this evidence we conclude that the conjecture put forth in [3] is not at all supported by the numerical simulations. The reasoning leading from maximum symmetry breaking of the external force to a maximum response of the system is of a linear response style and does not apply to the kind of nonlinear behavior that ratchet current generation represents. For the sake of reproducibility we provide the details of the numerical procedure we have followed to obtain Fig. 1.Simulations of the stochastic differential equation (1) have been performed using the second-order weak predictor-corrector method [10] with time step t =0.01, initial condition x(0) = 0, and a final integration time tf=200π/ω. The ratchet velocity vhas been computed using formula (2) averaging over 5000 realizations of the noise. For each γin Fig. 1we have obtained an entire curve v(η) for 100 values of ηin the interval [0,1]. Despite the average over such a large number of realizations, the resulting curves are still quite noisy, too much to 0 π/2 π3π/2 2π φ -0.06 -0.04 -0.02 0 0.02 0.04 0.06 v (a) 0 0.2 0.4 0.6 0.8 1 η -0.07 -0.06 -0.05 -0.04 -0.03 -0.02 -0.01 0 v (b) FIG. 3. (Color online) Plots of (a) vas a function of φfor η=2/3 and (b) ratchet current vas a function of η,forφ=π/2. Bullets are the numerical simulations of Eq. (1) with parameters ω=0.08π, σ=2, and γ=2. The yellow solid line in (a) is a fit to a sinusoidal function (v=−0.062 sin φ). The yellow solid line in (b) is a fit to −η2(1 −η)A(η), where A(η)isthe(2,2) Pad´ e approximant A(η)= (0.025 −0.113η+0.212η2)(1 −2.662η+2.011η2)−1. 066101-2 COMMENTS PHYSICAL REVIEW E 88, 066101 (2013) 0 π/2 π3π/2 2π φ -0.06 -0.04 -0.02 0 0.02 0.04 0.06 v FIG. 4. (Color online) Plot of vas a function of φfor η=2/3, γ=6, ω=0.08π,andσ=2. Bullets are the numerical simulations of Eq. (1). The yellow solid line is a fit to the two first harmonics of Eq. (4) (v=−0.079 sin φ−0.011 sin 3φ). reliably determine the value ηopt. For this reason we have recalculated the curves v(η) for another 100 values of ηin a narrower interval that clearly contains ηopt and have fitted a fourth-degree polynomial to the results (see Fig. 2for an example). The value of ηopt is obtained by optimizing this polynomial. This is how the points of Fig. 1have been obtained. As for the second conclusion of Ref. [1], aside from the evidence provided by Fig. 1that the numerical results are consistent with the ηopt =2/3 prediction of the theories in the limit of weak external forces, we can actually go further and show that a recent extension of the theory developed in Ref. [9], valid for arbitrarily large forces [11], fits perfectly with the results presented in [1,2]. For the case of harmonic mixing represented by Eq. (1), the theory predicts that vis given by the harmonic expansion v=∞  n=0 An(γ,η)γ6n+3η4n+2(1 −η)2n+1sin[(2n+1)φ],(4) where the coefficients An(γ,η) are functions of the squares of the amplitudes of the forcing harmonics, i.e., of γ2η2and γ2(1 −η)2. This implies that if the current is well described by one sinusoidal function, then v=A0(γ,η)γ3η2(1 −η)sinφ+O(γ9)(5) and A0(γ,η) should be well described by a bivariate quadratic polynomial, or any other approximant of an equivalent order, in γ2η2and γ2(1 −η)2. Figure 3(a) shows a fit of a sinusoidal function to the simulation data for vas a function of φobtained from Eq. (1) for η=2/3 and the other parameters as in Fig. 1 of [2]. It clearly shows that retaining only the first harmonic in (4) is enough to accurately reproduce the data. Thus v should conform to (5). Accordingly, we set φ=π/2 and fit the simulation results of vvs ηto a function of the form −η2(1 −η)A(η), where we take for A(η)a(2,2) Pad´ e approximant.1The result is plotted in Fig. 3(b) to show that this fit is a very accurate description of v(η) and therefore correctly predicts the deviation of ηopt from its weak force approximation ηopt =2/3. Finally, we would like to point out that the idea of an optimal shape of the external force is very difficult to reconcile with the current shape given by Eq. (4) because as soon as the amplitude of the force γbecomes sufficiently large, new harmonics will start modulating the shape of the current (Fig. 4clearly illustrates this effect). In this regime, only a very specific dependence of the coefficients Anwith γ, which does not occur in the case of Eq. (1), would yield the universality claimed in [1–3]. We thank Angel S´ anchez for his critical reading of this Comment. We acknowledge financial support through Grants No. MTM2012-36732-C03-03 (R.A.N.), No. FIS2011-24540 (N.R.Q.), and No. PRODIEVO (J.A.C.) from Ministerio de Econom´ ıa y Competitividad (Spain); Grants No. FQM262 (R.A.-N.), No. FQM207 (N.R.Q.), No. FQM-7276, and No. P09-FQM-4643 (N.R.Q. and R.A.-N.) from Junta de Andaluc´ ıa (Spain); Project MODELICO-CM (J.A.C.) from Comunidad de Madrid (Spain); and a grant from the Humboldt Foundation through Research Fellowship for Experienced Researchers No. SPA 1146358 STP (N.R.Q.). 1The only reason to use a Pad´ e approximant instead of a polynomial is that rational approximants are less prone to introduce spurious oscillations than high-degree polynomials. The choice of a (2,2) Pad´ e approximant is dictated by its having as many unknowns as a fourth-degree polynomial, so they both are approximants of the same order. [1] P. J. Mart´ ınez and R. Chac´ on, Phys. Rev. E 87, 062114 (2013). [2] P.J.Mart ´ ınez and R. Chac´ on, Phys.Rev.E88, 019902(E) (2013). [3] R. Chac´ on, J. Phys. A: Math. Gen. 43, 322001 (2010). [4] P. H¨ anggi and F. Marchesoni, Rev. Mod. Phys. 81, 387 (2009). [5] S. Flach, O. Yevtushenko, and Y. Zolotaryuk, Phys.Rev.Lett. 84, 2358 (2000). [6] S. Denisov, S. Flach, A. A. Ovchinnikov, O. Yevtushenko, and Y. Zolotaryuk, Phys.Rev.E66, 041104 (2002). [7] W. Wonneberger and H. J. Breymayer, Z. Phys. B 43, 329 (1981). [8] F. Marchesoni, Phys. Lett. A 119, 221 (1986). [9] N. R. Quintero, J. A. Cuesta, and R. Alvarez-Nodarse, Phys. Rev. E 81, 030102 (2010). [10] P. E. Kloeden and E. Platen, Numerical Solution of Stochastic Differential Equations (Springer, Berlin, 1995). [11] J. A. Cuesta, N. R. Quintero, and R. Alvarez-Nodarse, Phys. Rev. X 3, 041014 (2013). 066101-3