Force-linearization closure for non-Markovian Langevin systems with time delay
FForce-linearization closure for non-Markovian Langevin systems with time delay
Sarah A. M. Loos and Sabine H. L. Klapp
Institut f¨ur Theoretische Physik, Hardenbergstr. 36,Technische Universit¨at Berlin, D-10623 Berlin, Germany (Dated: July 3, 2018)This paper is concerned with the Fokker-Planck (FP) description of classical stochastic systemswith discrete time delay. The non-Markovian character of the corresponding Langevin dynamicsnaturally leads to a coupled infinite hierarchy of FP equations for the various n -time joint distri-bution functions. Here we present a novel approach to close the hierarchy at the one-time levelbased on a linearization of the deterministic forces in all members of the hierarchy starting from thesecond one. This leads to a closed equation for the one-time probability density in the steady state.Considering two generic nonlinear systems, a colloidal particle in a sinusoidal or bistable potentialsupplemented by a linear delay force, we demonstrate that our approach yields a very accurate rep-resentation of the density as compared to quasi-exact numerical results from direct solution of theLangevin equation. Moreover, the results are significantly improved against those from a small-delayapproximation and a perturbation-theoretical approach. We also discuss the possibility of accessingtransport-related quantities, such as escape times, based on an additional Kramers approximation.Our approach applies to a wide class of models with nonlinear deterministic forces. PACS numbers: 05.40.Jc, 02.30.Yy, 05.10.Gg, 02.50.Ey
I. INTRODUCTION
Many (non-)equilibrium systems from the macroscopicworld down to the quantum level are governed by dy-namical equations which involve both, noise and delay intime [1–3]. Noise due to imperfections or environmentalinfluences is essentially omnipresent in many real-worldsystems and experimental setups (e.g., at the cantileverof an Atomic Force Microscope [4]). On the mathemat-ical level it often also results from the presence of hid-den degrees of freedom, which have disappeared in thecontext of a coarse-graining procedure. Similarly, timedelay can have various sources. One example are finiteinformation processing and reaction times, which occurfor example in social systems (e.g., financial markets [5]or economic processes [6]). Such latencies as well ap-pear in neural systems due to finite signal transmissionand refractory times (which, e.g., become apparent inthe human pupil reflex [7], or in stick balancing experi-ments [8]). In biological systems, delays can be caused bymaturation times (as in population growth [9] or prey-predator dynamics [10]), or, on the biomolecular level,by biochemical reactions kinetics (such as the transcrip-tional and translational delay in gene networks [11]). Inlaser dynamics with optical feedback, delay is induced bythe traveling time of laser light on its round trip in thecavity [12]. On an even smaller scale, time delays oc-cur in quantum-optical systems coupled to a structuredphotonic reservoir [13].Beyond these intrinsic delays, time delay is an impor-tant issue in experimental setups with feedback control,where the delay is due to the finite time to proceed theoutput (plus the signal transmission times to obtain theinformation and feed it back into the system). Indeed,feedback control has emerged as an important tool tomanipulate small systems. Important representatives for such systems are colloids in a thermal bath under timedelayed feedback control [14–16]. All these examples in-volve noise and delay, illustrating that due to the gen-erality and omnipresence of such features, noisy systemswith time delays occur on all scales. In many of thesecases, the interplay of noise and delay moreover leads tointriguing dynamical behavior such as multistability andstochastic switches [12, 17].However, even in the classical case the mathematicaldescription of noisy systems with time delay continues tobe a major challenge. In presence of a discrete time de-lay τ , the standard mathematical equation for stochasticmotion, that is, the (under- or overdamped) Langevinequation (LE), becomes a stochastic delay differentialequation (SDDE), or, in a more physical language, a non-Markovian Langevin equation. The usual concepts andsolutions for Markovian LEs then do not apply and in-deed, the development of a general solution method forSDDEs is yet an unresolved problem [18]. Moreover, con-trary to the Markovian systems, the route from the LEto a Fokker-Planck equation (FPE) for the correspond-ing probability density is more involved [18–21]. Thetime delay alone leads to an infinite hierarchy of cou-pled FP equations for the n -time (joint) probability den-sities [18–21], for which as well no general solution hasbeen found to date. Similar hierarchy problems also oc-cur in other contexts in statistical physics. A well knownexample is the Bogoliubov–Born–Green–Kirkwood–Yvon(BBGKY) hierarchy, in which the time evolution of theone-particle density depends upon the two-particle den-sity and so on. For this problem, various closing strate-gies have been proposed, such as the simple mean-field(factorization) approximation and the more sophisticateddynamical density functional theory (involving an adia-batic approximation) [31]. Another closure example isthe mode coupling theory for glassy systems [32]. We a r X i v : . [ c ond - m a t . s t a t - m ec h ] J u l stress, however, that in all these examples the hierar-chical structure emerges due to (conservative) particle-particle interactions. In contrast, the hierarchy of FPEsappearing in time-delayed systems – although having asimilar structure – arises due to delay-induced temporal correlations.These aspects are also of major relevance in the emerg-ing field of stochastic thermodynamics [22–24]. Indeed,concepts like entropy production, fluctuation theoremsand information exchange have been widely formulated(and tested experimentally [25, 26]) for Markovian noisysystems, for which the correspondence between the var-ious levels of description (LE, FPE, path integrals) canreadily be utilized. For example, the entropy productionfrom the system can be directly calculated via the prob-ability density and the corresponding probability current[22]. For non-Markovian systems, all these concepts haveto be revisited (see [18, 27] for recent work in this direc-tion).The only solvable class of delayed noisy systems arethose with linear forces, where all the n -time probabil-ity densities are given by multivariate Gaussian distribu-tions. Here, the system can be solved exactly on the LE[28] and the FPE level [21], or by utilizing both levelsof description [29]. In the more general case of nonlin-ear systems, one has to rely on approximations [1]. Forcontinuous systems, the two most established strategiesare a Taylor expansion of the LE in τ [19, 30] and a per-turbation theory on the level of the FPE [20], where theentire delay force is treated as a small deviation from theMarkovian dynamics.In this work, we propose a novel FPE-based approach,which we call “force linearization closure” (FLC). Themain idea is to linearize the deterministic force in allmembers of the FPE hierarchy starting from the secondone. We then utilize the (Gaussian) solution for the lin-ear, delayed higher-order FPEs to obtain a closure ofthe FPE for the one-time probability density in the non-linear case. Our strategy can be applied whenever thelinearized system has a stable steady state.We apply this concept to a classical system, namely,an (overdamped) colloidal particle subject to a nonlinearstatic force and a time-delayed force representing a feed-back control. The control “target” is the particle positionwhich, due to the typical (micron-scale) size of a colloid,is readily accessible in real space experiments and simula-tions [14–16]. Feedback control can be implemented, e.g.,as a co-moving “optical tweezer” [24, 25, 33–36], whichcan be well approximated by a co-moving quadratic po-tential [25, 33, 37]. Depending on the experimental setup,a time delay naturally arises during the position mea-surement and the adjustment of the tweezer, or it can beintentionally implemented as a feature. As a result, thesystem is subject to a linear time-delayed feedback force.Regarding the nonlinear static force (or potential, re-spectively) we consider two paradigmatic examples: first,a periodic sinusoidal (“washboard”) potential, which isfrequently used to model diffusion and transport on rough surfaces [38–40] and can be realized experimen-tally e.g. by optical landscapes [39]. The impact of timedelay on diffusion in washboard potentials has alreadybeen studied on the basis of a heuristic delayed FPE[24, 35, 36] and by simulations [16]. Our second exampleis a colloidal particle in a double-well potential. Bistablenoisy systems often serve to study escape problems [41–44]. Bistable systems with delay have been studied the-oretically, e.g., regarding the Kramers rate [17, 30, 45]and in the context of coherence resonance (see the studyof Tsimring and Pikovsky [46], who provided a solutionbased on a discretization procedure, and subsequent work[17, 46–48]). Experimental realizations have been sug-gested in Refs. [45, 49], which involve the polarizationdynamics of lasers with optical feedback (where τ is as-sociated with the external cavity length).The remainder of this paper is organized as follows. Westart by briefly reviewing the Langevin and the Fokker-Planck description of noisy overdamped systems with dis-crete time delay, as well as the two main earlier approxi-mations for the one-time probability density ρ (Sec. II).We then introduce the force linearization closure and pro-vide the main steps to derive its key equation, a closedFPE for ρ (Sec. II C). In the second part of the pa-per, we apply the FLC to the two examples mentionedabove and compare the results with quasi-exact numeri-cal results obtained from Brownian dynamics (BD) sim-ulations, i.e., direct simulations of the delayed Langevinequation. We demonstrate that the FLC renders a verygood approximation of the one-time probability density,which is, moreover, more accurate than the predictionsfrom the two main earlier approaches. Furthermore, webriefly discuss the possibility to estimate escape times(which are rigorously connected to two-time probabilitydensities) within the novel approach (Sec. III D). We sum-marize and conclude in Sec. V. The paper contains fourappendices with additional theoretical results (e.g., thesecond member of the FPE hierarchy), and technical as-pects of the calculations. II. THEORETICAL FRAMEWORKA. Delayed Langevin and Fokker-Planck equation
We consider an overdamped Brownian particle at thetime-dependent stochastic position χ ( t ) which representsthe dynamical quantity of interest. The particle movesin a static external potential V s ( x ), with spatial variable x ∈ Ω, where Ω denotes the spatial domain of the sys-tem. The particle is further subject to a “delay force” F d , which depends on the instantaneous and on the de-layed particles position, i.e., at a given time t on χ ( t )and χ ( t − τ ), with τ > F ( x, x τ ) = F s ( x ) + F d ( x, x τ ) , (1)with x τ ∈ Ω being a second spatial variable needed dueto the two involved particle positions in F d , and F s ( x ) = − ∂ x V s ( x ) being the negative of the spatial derivative of V s . The delayed Langevin equation (LE) readsd χ ( t ) / d t = γ − F [ χ ( t ) , χ ( t − τ )] + (cid:112) D Γ( t ) , (2)where γ is the friction coefficient of the surroundingmedium at temperature T and Γ( t ) denotes the Langevinforce. The latter introduces additive Gaussian whitenoise to the system, i.e., (cid:104) Γ( t ) (cid:105) = 0 ∀ t and (cid:104) Γ( t )Γ( t (cid:48) ) (cid:105) = δ ( t − t (cid:48) ) ∀ t, t (cid:48) , with (cid:104) ... (cid:105) denoting the ensemble averageand δ ( t ) the Delta distribution. Further, D is the short-time diffusion coefficient satisfying D = k B T /γ (with k B being the Boltzmann constant) [43]. We measure thetime in units of the “Brownian” time scale τ B = σ /D ,where σ is the particle diameter. This is the time in whicha free Brownian particle (no external forces) travels overa distance equal to its own size.As usual for delay equations, the evolution of the dy-namical variable (here χ ) depends on a given historyfunction φ , which serves as an initial condition χ ( t ) = φ ( t ) , ∀ t ∈ [ − τ, ρ ( x, t ) = (cid:104) δ ( x − χ ( t )) (cid:105) ina similar manner to the Markovian case. We particu-larly refer to Ref. [20], which presents a derivation basedupon Novikov’s theorem [50]. Novikov’s theorem gener-ally links the variational derivative of an arbitrary func-tional of the Langevin force, Λ[Γ], to the correlations be-tween that functional and the Langevin force. For Gaus-sian white noise it reads [20, 50]: (cid:10) Λ[Γ]Γ( t ) (cid:11) = (cid:42) δ Λ[Γ] δ Γ( t ) (cid:43) . (3)A central step in the derivation in [20] is the us-age of Novikov’s theorem to express the correlations (cid:10) { δ ( x − χ ( t )) δ ( x τ − χ ( t − τ )) } Γ( t ) (cid:11) by the related vari-ational derivatives. For the case of additive noise, thedelayed FPE for ρ is given by [19, 20, 29] ∂ t ρ ( x, t ) = − γ − ∂ x (cid:104) (cid:98) F ( x, t ) ρ ( x, t ) (cid:105) + D ∂ xx ρ ( x, t ) , (4a)where (cid:98) F ( x, t ) = F s ( x ) + (cid:90) Ω F d ( x, x τ ) ρ c ( x τ , t − τ | x, t )d x τ . (4b) For the sake of a shorter notation, we drop here andin the following the dependency on the history func-tion, i.e., ρ ( x, t ) ≡ ρ ( x, t | φ ) and ρ c ( x τ , t − τ | x, t ) ≡ ρ c ( x τ , t − τ | x, t ; φ ). Comparing Eq. (4a) to an ordinaryFPE for Markovian systems, one notes the appearanceof a delay-averaged drift term (cid:98) F instead of the usualdrift term, which only depends on instantaneous quan-tities. Specifically, the delay-averaged drift involves theconditional probability ρ c , which is related to the two-time (joint) probability density via ρ ( x, t ; x τ , t − τ ) = ρ c ( x τ , t − τ | x, t ) ρ ( x, t ) [with (cid:82) Ω ρ ( x, t ; x τ , t − τ )d x τ = ρ ( x, t )]. Equation (4) is hence not self-sufficient. Thedelayed FPE for ρ , on the other hand, involves ρ (as explicitly shown in the Appendix A), and so forth.Thus, an infinite hierarchy of equations emerges, whose n +1st member depends on the n -time probability den-sity ρ n ( x, t ; ... ; x nτ , t − nτ ). By finding suitable approxi-mations for ρ c , Eqs. (4) can be closed and the hierarchytruncated. This is our objective in the present work. B. Earlier approaches to approximate the one-timeprobability density
One such approach, involving a first-orderperturbation-theoretical (PT) ansatz for the den-sity on level of the FPE, was previously introduced inRef. [20]. Within PT, the delay force F d is regarded asa small perturbation to the non-delayed dynamics, i.e., | F d | (cid:28) | F s | . The resulting first-order equation has thesame form as Eq. (4a), but the (cid:98) F -term [Eq. (4b)] containsthe conditional probability density with respect to theunperturbed, i.e., non-delayed, system, which we willrefer to as ρ F d ≡ c . The latter follows a closed (Markovian)FPE, and consequentially, the combination of the twoFPEs is self-sufficient. From a practical perspective, thePT approach is particularly appealing for systems wherean analytical expression for ρ F d ≡ is available. However,this is the case only for very few nonlinear static forces(even in the absence of any additional (delay) forces).For this reason, the application of the PT approachoften requires additional approximations. In the presentstudy we use two different approaches, which we willspecify in Sec. III A.Yet another approach involves a small delay expansionon the level of the LE [1, 19, 30]. More specifically, theauthors of Ref. [19, 30] suggest a Taylor expansion upto linear order of the total deterministic force (and ofthe noise intensity when delayed noise is considered) inpowers of τ . The Taylor expansion is somewhat prob-lematic, since it involves the derivative of the position χ ,which is a stochastic variable and hence not continuouslydifferentiable. Performing this Taylor expansion ad hoc hence implies a certain inaccuracy, which may becomeespecially apparent in the range of large noise intensities(large D ), see [20] for a discussion.The novel approach introduced in the present paperworks on the basis of the FPE, like the perturbation the-ory. As we will demonstrate, our ansatz provides an im-proved approximation of the one-time probability densityfor the considered nonlinear static forces. C. Force linearization closure
In the following, we focus on non-equilibrium steadystates (NESS). Indeed, many Langevin systems de-scribed by Eq. (2) automatically reach a NESS, definedby ∂ t ρ , ss ( x, t ) = 0. Steady states generally allow for moreanalytical treatment than transients. Throughout thiswork we denote steady state quantities with the subscript“ss”. Since the steady state conditional probability onlydepends on the time difference, and not on the instancesof time themselves, we will further use the shortened no-tation: ρ c , ss ( x τ | x ; τ ) = ρ c , ss ( x τ , t − τ | x, t ).The basis of our approach is that, within a NESS, theFPE hierarchy can be solved exactly when all determin-istic forces are linear in x and x τ [28, 29]. A further, yetless restrictive requirement is that the system obeys nat-ural boundary conditions, i.e., lim x →±∞ ρ ( x, t ) = 0 andΩ = R . Given this background, the main idea of ourapproach is to achieve a closed approximate FPE for ρ by linearizing the deterministic forces in all members ofthe infinite FPE hierarchy apart from the first one. Dueto the involved linearization procedure, which we outlinein detail below, we call our approach “force linearizationclosure” (FLC).
1. Linearization of the deterministic forces
We start by considering the (time-dependent) energylandscape resulting from the total force (1). As a simpleestimate for the total steady state energy landscape, weassume at this step the system to be at rest, i.e., χ ( t − τ ) ≡ χ ( t ), such that the total (static) potential is given by V STAT ( x ) = − (cid:90) F ( x, x τ = x )d x. (5)In general, V STAT may have multiple (local) minima (de-fined by F = 0 and F (cid:48) > x i with an integer i ∈ I . The index set I hence containsone element for each local minimum. To formulate theFLC, we split the spatial domain Ω into non-overlappingintervals Ω i , such that both boundaries of Ω and each lo-cal maximum of V STAT represents a bound of an interval,and every Ω i contains exactly one minimum x i .By this procedure, we obtain a splitting of Ω into sub-domains Ω i , whose union is again the entire spatial do-main: Ω = (cid:84) i ∈I Ω i . For static potentials with a sin-gle minimum x j , this means Ω j = Ω. For each i , wethen perform a Taylor expansion of the entire deter-ministic force F = F s + F d in both spatial variables x and x τ . More specifically, we expand around the devia-tions with the enclosed minimum, i.e., ∆ x i = x − x i and ∆ x iτ = x τ − x i . Neglecting all terms of quadratic orders O (∆ x i ), O (∆ x i ∆ x iτ ), O (∆ x i τ ), or higher, we obtain F lin ,i ( x, x τ ) = − α i ∆ x i − β i ∆ x iτ , (6)where − α i and − β i are the first-order derivatives of F with respect to x and x τ , respectively, evaluated at theminimum x i . [Note that since we expand around theminima of V STAT (5), we always have F ( x i , x i ) = 0, suchthat the constant terms vanish.]This linearization procedure yields an approximationof the steady state energy landscape composed of a se-quence of quadratic polynomials with time-dependentcenters, each subdomain Ω i reaching from one local max-imum to the following one (or to a bound of the entirespatial domain).
2. Analytical solution for linearized forces
Now, we will turn back to the full system which in-volves thermal noise. Without further reasoning, onewould expect that the linearization of the determinis-tic forces renders a good approximation of the stochas-tic dynamics, whenever the probability density close tothe minima of the (approximate) total potential in thesteady state [Eq. (5)] is large. This is, for instance, thecase, when the potential barriers are high compared tothermal fluctuations.The following steps are performed separately for eachsubdomain Ω i . We first apply the linearization F ≈ F lin ,i [see Eq. (6)] to the force terms in all membersof the infinite FPE hierarchy starting from the second,i.e., the equation for the two-time steady state density ρ , ss ( x, t ; x τ , t − τ ) (see Appendix A for the general fromof this FPE). Assuming natural boundary conditions atevery subdomain bound, the FPEs for all ρ n , ss can besolved by multivariate Gaussian distributions [28]. Thesecond member of the FPE hierarchy in the linear caseand the corresponding analytical solutions ρ lin2 , ss and ρ lin3 , ss are given in the Appendices B and C. For each subdo-main, we thus have access to the function ρ , ss ≡ ρ lin ,i , ss (the explicit formula for ρ lin ,i , ss can also be found in [51]).The corresponding conditional probability density reads ρ lin ,i c , ss ( x τ | x ; τ ) = (cid:114) π K i d i (cid:16) − d i (cid:17) exp (cid:34) ∆ x i K i (cid:35) × exp (cid:104) d i (cid:16) x i ∆ x iτ d i − ∆ x i − ∆ x iτ (cid:17)(cid:105) (7)with the coefficients ω i ≡ ω ( α i , β i ) = (cid:112) ( α i ) − ( β i ) /γ ∈ C , (8a) K i ≡ K ( α i , β i , τ ) = D γ + ( β i /ω i ) sinh( τ ω i ) α i + β i cosh( τ ω i ) , (8b) d i ≡ d ( α i , β i , τ ) = ( β i ) K i / β i K i ) − ( D γ − α i K i ) , (8c) d i ≡ d ( α i , β i , τ ) = ( D γ − α i K i ) / ( β i K i ) . (8d)We recall ∆ x i = x − x i , where x i is the enclosed min-imum of V STAT [Eq. (5)]. Note that ω i becomes imag-inary if | α i | < β i , such that the hyperbolic functionsin K ( α i , β i , τ ) [see Eq. (8b)] convert to trigonometricones. In this case, there exist critical τ values τ c ω =arccos( − α i /β i ) + 2 πκ, ∀ κ ∈ Z , for which K i diverges,and no stable NESS exists. There is also no stable state,when − β i ≤ α i , or τ α i ≤ −
1. Only if − α i < β i < α i , asteady state is approached for all τ . A derivation anddiscussion of the steady state conditions is presented in[28].With Eq. (7), the first member of the approximate FPEhierarchy (where the force linearization is applied to thesecond and all higher members) accordingly reads γD ∂ xx ρ FLC1 , ss ( x ) = ∂ x [ (cid:98) F FLC ( x ) ρ FLC1 , ss ( x )] , (9a)where (cid:98) F FLC ( x ) = (cid:90) Ω F ( x, x τ ) ρ lin ,i c , ss ( x τ | x ; τ )d x τ , ∀ x ∈ Ω i . (9b)In Eq. (9a), ρ FLC1 , ss is the steady-state one-time probabilitydensity obtained within the FLC approach. Notice thatthe deterministic force F appearing in the integral inEqs. (9b) is not linearized.Equations (7-9b) [together with the linearization rule(6)] form a closed set of equations. The here presentedclosure of the FPE hierarchy is the cornerstone of ourapproach. Since we have a general expression for ρ lin ,i c , ss [Eq. (7)] with which (cid:98) F FLC can be calculated separately,the FLC formally converts the delayed FPE for ρ , ss itselfinto a closed, quasi-Markovian one [Eq. (9a)]. The factthat the delay is still present in our approximate FPE, be-comes apparent by considering the special case of a linearforce F . Then, Eqs. (9) coincide with the correspondingexact delayed FPE. Furthermore, one can easily see fromEqs. (9) that the usual FPE is recovered when the delaytime vanishes, or when F becomes independent of x τ : inboth cases (cid:98) F FLC = F .
3. Vanishing steady state probability current
The key equations of the FLC [Eqs. (7-9b)] are validfor steady states of, in principle, arbitrary systems whichcan be meaningfully linearized according to the proce-dure described in Sec. II C 1. A particularly simple situ-ation arises, if additionally the steady-state probabilitycurrent given by J = [ (cid:98) F FLC ( x ) − γD ∂ x ] ρ FLC1 , ss ( x ) vanishes.If J = 0, the formal solution of Eq. (9a) takes the simpleBoltzmann-like form ρ FLC1 , ss ( x ) = Z − exp (cid:2) − V FLCeff ( x ) / ( k B T ) (cid:3) (10)with the effective static potential V FLCeff ( x ) = − (cid:90) x ˆ x d x (cid:48) (cid:98) F FLC ( x (cid:48) ) , (11) where (cid:98) F FLC is given in Eq. (9b), and ˆ x ∈ Ω is arbitrarybut fixed. Here and in the following, we denote the nor-malization constant by Z .We stress that the assumption J = 0 simplifies the anal-ysis, but it is not a necessary condition. Also for non-zerocurrents, Eq. (9a) can be treated using standard tech-niques for (Markovian) Fokker-Planck equations [43, 44]. III. APPLICATIONS
We now apply the FLC to two generic examples, in-volving a multistable and a bistable static potential com-bined with a linear delay force. In the subsequent sec-tion III A, we first define this force and provide someresults which apply to both model systems. In Secs. III Band III C we then present results from the FLC and com-pare its performance to the two main approaches knownfrom the literature, i.e., the perturbation theory (PT)and the small delay expansion. As a test of the differ-ent approximations, we provide numerical results fromBrownian dynamics (BD) simulations. Details about thenumerical methods are given in the Appendix D.
A. Linear delay force, general results
We consider a linear delay force, with amplitude k ≥ F d ( x, x τ ) = − k ( x − x τ ) , (12)which vanishes for k → τ →
0. The particular forcein Eq. (12) can be associated with the delayed confiningpotential V d ( x, x τ ) = ( k/ x − x τ ] . Such quadratic feed-back potentials are commonly used to model optical traps[25, 33, 37], which are implemented in many experimen-tal setups to control colloidal particles [25, 33, 34]. Froma more general perspective, the delay force in Eq. (12) isof Pyragas type [52], and has been extensively studiedin the context of chaos control [3]. It is important tonote that the FLC generally also applies to systems withnonlinear delay forces.For both exemplary systems, we use natural boundaryconditions. Therewith, the steady state probability cur-rent vanishes irrespective of the particular form of thestatic potential V s , and solution (10) readily applies. Wecan give a general expression for ρ FLC1 , ss as follows. Thedelay force is already linear, and yields β i = − k, ∀ i ∈ I ,where β i and k are the coefficients appearing in Eqs. (6)and (12), respectively. After performing several Gaussianintegrals, Eq. (9b) yields (cid:98) F FLC ( x ) = − ∂ x V s ( x ) − k ∆ x i + kxd i (cid:113) K i d i × (cid:113) − d i exp (cid:110)(cid:104) d i d i − d i + 1 / (2 K i ) (cid:105) ∆ x i (cid:111) . (13) (cid:98) F FLC can be further simplified by using the identity d i d i − d i + 1 / (2 K i ) = 0 , (14)where the quantities d i and d i are given in Eqs. (8c) and(8d). This yields the piecewise defined effective potential V FLCeff ( x ) = V s ( x ) + k/ (cid:0) − | d i | (cid:1) ∆ x i , (15)for x ∈ Ω i . Inserting Eq. (15) into Eq. (10) one obtains ρ FLC1 , ss . Furthermore, since F d ( x, x ) = 0 ∀ x ∈ Ω, our es-timate for the steady state energy landscape coincideswith the static potential V STAT = V s , and hence x i andthe bounds of Ω i are readily determined by the extremaof V s . Thus, for a given V s , one only needs to calcu-late α i by linearizing F s , and therewith the coefficient d i [Eq. (8d)].Before specifying V s , we review some results from thetwo main earlier approaches, for the case of the lineardelay force [Eq. (12)] (and natural boundary conditions).Just like the FLC [Eq. (10)], the small delay expansiongives rise to a Boltzmann-distributed one-time probabil-ity density ρ approx1 , ss ( x ) = Z − exp [ − V approxeff ( x ) / ( k B T )] . (16)Here, the effective potential V approxeff ≡ V small τ eff reads [20] V small τ eff ( x ) = V s ( x ) / (1 − kτ ) . (17)The application of the perturbation-theoretical approach,on the other hand, is not straightforward. The closedapproximate FPE from the PT involves the conditionalprobability of the corresponding unperturbed ( F d ≡ ρ F d ≡ , ss ( x τ | x ; τ ). Since for most static potentials(including the ones we will consider in Secs. III B andIII C) no analytical expression for this quantity is avail-able, further approximations become inevitable. As sug-gested in Ref. [20], one can for this purpose utilize the“short time propagator” [43] ρ c ( x τ | x ; τ ) =1 √ πD τ exp (cid:34) − { x τ − x − F s ( x ) τ /γ } D τ (cid:35) , (18)which is derived and discussed in Ref. [43] on p. 73 andprior. The usage of the short time propagator impliesa first-order approximation in τ [20, 43]. The resultingone-time probability density is also the Boltzmann dis-tribution Eq. (16) with V approxeff ≡ V PTeff given by V PTeff ( x ) = (1 + kτ ) V s ( x ) . (19)Interestingly, when the short time propagator is used for ρ F d ≡ , ss , the PT approach obviously yields qualitativelyvery similar results as the small τ expansion. In particu-lar, both approaches render effective potentials that areproportional to the static potential [Eqs. (17) and (19)].This means, also the approximate densities have the samefunctional form as in the case F d ≡
0. In fact, the onlyremaining effect of the delay force, according to these ap-proximations, is a constant ( x -independent) factor within the exponent of ρ approx1 , ss . Moreover, the small τ expan-sion and the PT with short time propagator even becomeequivalent for small τ k , since then these factors coincide1 / (1 − kτ ) = (1 + kτ ) + O ([ τ k ] ).Alternatively, we propose to approximate ρ F d ≡ , ss bythe conditional probability from the corresponding un-perturbed ( F d ≡ V PTeff ( x ) = V s ( x ) + ( k/ × (cid:8) − exp [ − ( α i − k ) τ /γ ] (cid:9) ∆ x i , ∀ x ∈ Ω i . (20)As opposed to the result obtained with the short-timepropagator [Eq. (19)], the effective potential in Eq. (20)has a different functional form than V s ( x ). More specif-ically, the delay force now effectively adds to the staticpotential a fixed quadratic potential around the minimaof V s . Since the linear delay force [Eq. (12)] is indeedexpected to trap the particle in a quadratic confining po-tential (with history-dependent center), this effective po-tential appears to be more realistic than Eq. (19). How-ever, the history-dependency of the position of the “trap”imposed by F d is, of course, also not captured by this ap-proximations.As a first consistency check, one readily sees fromEqs. (15-17),(19),(20) that all four approximations giverise to Boltzmann-distributed density profiles, which, inthe limits of vanishing delay force ( τ → k → ρ , ss ∝ exp[ − V s ( x ) /k B T ]. How-ever, they yield different effective potentials for finite de-lay forces, which we will compare in the following.Finally, we note that considering the force (12) ratherthan a force with two different amplitudes F d ( x, x τ ) = − k (cid:48) x − kx τ , does not imply a loss of generality concerningthe results in this section. A minor adjustment requiredif k (cid:48) (cid:54) = k , is that x i are then the minima of V STAT (cid:54) = V s . B. Periodic static potential
We start by considering the periodic “washboard” po-tential V s ( x ) = − (∆ V/
2) cos( x/σ ) (21)with barriers of height ∆ V and minima at x i = 2 πσi and i ∈ I = Z . Particles in sinusoidal potentials rep-resent a well-studied paradigm to model transport onrough surfaces [39]. Linearizing the static potential yields α i = ∆ V/ (2 σ ) + k (and β i = − k ) for all i ∈ I . On all sub-domains Ω i = [(2 i − πσ, (2 i + 1) πσ ], the density is hencegiven by the Boltzmann distribution (10) with effectivepotential (15), where V s is given in Eq. (21) and d i = γω cosh( τ ω ) − [∆ V/ (2 σ ) + k ] sinh( τ ω ) γω − k sinh( τ ω ) (22) x/σ ρ , ss ( x ) FLCPT OUPT s.t.small τF d ≡ BD FIG. 1. (Color online) One peak of the periodic steadystate probability density ρ , ss in the “washboard” with de-lay force. Red symbols: numerical data (BD), thick red line:force linearization closure (FLC), blue dashed line: pertur-bation theory with Ornstein–Uhlenbeck approximation (PTOU), blue thin line: PT with short time propagator (PTs.t.), dashed-dotted line: small delay expansion (small τ ),and solid black line: density without delay force ρ , ss ( x ) ∝ exp {− V s ( x ) / ( k B T ) } ( F d ≡ k = 8 k B T /σ and τ = τ B , respectively.Here and in the following plots the barrier height is set to∆ V = 8 k B T . with γωσ = (cid:112) (∆ V/ σ ) + k ∆ V and ∆ x i = x − x i .Figure 1 shows the density profile in one potential wellgenerated with BD simulations of the delayed LE (2), onthe one hand, and the corresponding FLC approxima-tion (22), on the other hand, for an exemplary parameterchoice. Here and throughout the entire paper, we set thebarrier heights to ∆ V = 8 k B T . One clearly sees that theFLC generates a very good approximation of the proba-bility density distribution. Similar convincing results arefound for all other considered values of τ and k in thetested range τ /τ B ∈ [10 − ,
10] and k/ ( σ − k B T ) ∈ [0 , µ n of the dis-tributions within one potential well (for x values be-tween two maxima of V s ), e.g., of ρ , ss ( − πσ ≤ x ≤ πσ ).Due to the spatial inversion symmetry of the total po-tential around the enclosed minimum, all odd momentsvanish. We hence consider the second central moment µ = (cid:104) ( χ − µ ) (cid:105) ss (with µ = (cid:104) χ (cid:105) ss ) of one peak of thedistribution. Figure 2 shows µ vs. the delay time τ fortwo exemplary values of delay force amplitude. We findthat in the presence of the delay force (12), the variance ofthe distribution is reduced. This is not surprising, sincethe delay force arises from a quadratic potential which“traps” the particle and hence enhances the confining ef-fect of V s . Within the considered parameter regime, thevariance is seen to decrease with increasing τ , until a sat-uration value is approached at τ Wsat ≈ . τ B (due to thestrongly increasing simulation times, we did not simulatemuch larger τ ). On the level of the probability densities and the FPE,the saturation of µ at large τ indicates that in the con-sidered range of delay times, the delay-averaged drift (cid:98) F , and hence ρ c , ss ( x τ | x ; τ ), are essentially constant for τ > τ Wsat . In other words, the conditional probabilityfor the time difference τ remains essentially unchanged,when τ is further increased. This suggests some kindof relaxation mechanism within the valley, where thehistory of the stochastic process becomes “washed out”on the level of ensemble averaged quantities, like theprobability densities. In this context it is interesting tonote that, for the case without delay force ( k = 0), therelaxation time within a potential well is of the order τ Wir ≈ γ/V (cid:48)(cid:48) s ( x min ) = 2 τ B k B T / ∆ V (see Ref. [44], p. 348).In the present case, ∆ V = 8 k B T , that is τ Wir ≈ . τ B .The saturation of µ thus sets in after the density relax-ation within a potential well, consistent with our expec-tation.A different situation occurs for much larger τ valuesthan the ones considered in Fig. 2, in particular, when τ gets into the range of the mean escape times (numericalresults for the escape times are provided in Sec. III D).Then we expect again a τ -dependency of ρ c , ss and there-with of µ , since the transport between the potential val-leys becomes important.We now compare the density from the FLC approachwith corresponding data from the small delay (Taylor)expansion [Eq. (16) and (17)] and from the PT approach.Within the latter, we either use the short time propagator[see Eq. (19)], or, the corresponding Ornstein–Uhlenbeckapproximation [see Eq. (20) and above]. As visible inFig. 1 and 2, the PT generally overestimates the height ofthe density peak. Approximating the conditional proba-bility with the short time propagator yields worse resultsthan using ρ c , ss of the corresponding Ornstein–Uhlenbeckprocess. This observation matches our expectations, seeSec. III A. At least when the Ornstein–Uhlenbeck approx-imation is used, the PT reproduces the quantitative be-havior of the function µ ( τ ), see Fig. 2. On the contrary,the small delay expansion fails completely for larger delaytimes.We conclude that the FLC generates the best approx-imation of the steady-state density distribution in theperiodic potential, especially in the regime of large τ orlarge k . This is the regime where the “perturbation”,i.e. | F d | , is not small compared to the force applied bythe static potential, such that the basic assumption ofthe perturbation theory is not fulfilled. The small de-lay expansion, on the other hand, involves a truncatedTaylor expansion in τ , making it plausible that also thisapproach fails for large delay times. -2 -1 0 1 log ( τ/τ B ) s e c o n d m o m e n t µ k =2 k B T σ − k =4 k B T σ − FLCPT OUPT s.t.small τ BD FIG. 2. (Color online) Second moment µ of the steady stateprobability density within one well of the “washboard” poten-tial, e.g., of ρ , ss ( − πσ ≤ x ≤ πσ ). Symbols show results fromBrownian dynamics (BD), and thick solid lines from force lin-earization closure (FLC). Cyan color is k = 2 k B T /σ . Blackis k = 4 k B T /σ , with thin solid line: perturbation theorywith short time propagator (PT s.t.), dashed line: perturba-tion theory with Ornstein–Uhlenbeck ρ F d ≡ , ss (PT OU), anddot-dashed line: small τ Taylor expansion (small τ ). C. Bistable static potential
As a second example, we consider the double-well po-tential V s ( x ) = ∆ V [( x/σ ) − x/σ ) ] , (23)with a barrier of height ∆ V (set to ∆ V = 8 k B T ) at x/σ = 0and two minima at x i = iσ with i ∈ I = {− , } .A Brownian particle in a double-well potential is ageneric model for a bistable noisy system [43, 44], whichin recent years has been considered also under the im-pact of a linear delay force [17, 46–49]. As opposed tothe corresponding expression given in Eq. (12) (whichcorresponds to a moving “optical tweezer”), these ear-lier theoretical studies consider a slightly different delayforce involving solely the delayed particle position, i.e., F d ( x, x τ ) = − kx τ . Our approach can, however, as wellbe applied to such a delay force.The force linearization (according to the procedure de-scribed in Sec. II C 1) yields α i = 8 ∆ V/σ + k for both i (and β i = − k ). For x/σ ∈ [ −∞ ,
0] and x i = − = − σ , or x/σ ∈ [0 , ∞ ] and x i =10 = σ , the effective potential V FLCeff is given by Eq. (15) with d i = γω cosh( τ ω ) − [(8∆ V/σ ) + k ] sinh( τ ω ) γω − k sinh( τ ω ) (24)with γωσ = (cid:112) (8∆ V/σ ) + k ∆ V .Also for this model, the FLC clearly renders a verygood approximation of the one-time probability density,as shown in Fig. 3 for an exemplary parameter choice.Because the second moment is barely affected by the delay force, we here consider the third moment of onepeak of the bimodal distribution, i.e., ρ , ss ( x > µ x> = (cid:82) x> [( x − µ ) / √ µ ] ρ , ss ( x )d x , i.e., the thirdcentral moment normalized with the standard deviation.The comparison with the numerical data in Fig. 4 revealsthat the FLC approach renders reasonably good predic-tions of µ x> . We moreover find convincing results forall other tested parameters within the considered ranges τ /τ B ∈ [10 − ,
5] and k/ ( σ − k B T ) ∈ [0 , k val-ues so high that there is no second minimum of the totalpotential, when χ ( t − τ ) rests in the first minimum. Forall these qualitatively different cases, the agreement be-tween FLC and numerical results is good. However, itis not as accurate as that for the second moment in thecase of the periodic potential (see Fig. 2). One reason forthat might be the fact that the approximate potential is,by construction, symmetric with respect to the enclosedminimum (within each subdomain). In this sense, thedouble-well potential, which is asymmetric around eachminimum, is not as well approximated as the symmetric“washboard” potential. This manifests in a larger magni-tude of the higher order Taylor terms that are neglectedwithin the linearization procedure (6). We also note thatthe third moment is per se very sensitive against inac-curacies, since it involves cubic terms in the relative po-sition. The deviations between the FLC and the exactresult are hence expected to be larger.Regarding the impact of increasing k and τ , we seefrom Fig. 4 that the skew µ x> becomes smaller, whichmeans that the probability distribution becomes moresymmetric. This is not surprising since the pure delay po-tential is quadratic (i.e., symmetric) in the system statevariable. Moreover, similar to the second moment in thecase of the “washboard” potential (see Fig. 2), µ x> ap-proaches a saturation value (within the considered pa-rameter range). The corresponding delay time is aboutone order of magnitude smaller than in the “washboard”case: τ Dsat ≈ . τ B for ∆ V = 8 k B T . This value of τ Dsat is, as in the case of the “washboard”, significantly largerthan that for the relaxation time τ Dir within a well. Forthe bistable system at k = 0 one finds [44] (see p. 348) τ Dir /τ B ≈ / k B T / ∆ V ) = 1 /
64 at ∆ V = 8 k B T . Thus,the saturation of µ can be explained by the same argu-ments as in the case of the “washboard” potential.Figure 3 and 4 also show the results according to thePT [Eqs. (19) and (20)], and from small delay expansion,see Eq. (17). Very similar to the case of the “washboard”,only the PT with Ornstein–Uhlenbeck approximation iscapable of reproducing the quantitative behavior (espe-cially for large τ ), and the FLC clearly provides the bestapproximation of the one-time probability. x/σ ρ , ss ( x ) FLCPT s.t.PT OUsmall τF d ≡ BD FIG. 3. (Color online) One peak of the bimodal steady stateprobability density ρ , ss in the double-well potential. Colorcode as in Fig. 1. The delay force amplitude and delay timeare set to k = 18 k B T /σ and τ = 0 . τ B , respectively. D. Escape times
1. The Kramers-FLC estimate
The escape time τ K denotes the average time χ ( t )spends in the vicinity of a potential minimum, until itleaves the valley by jumping to an adjacent one. Byconstruction, mean escape times involve two locationsand two instances in time, and hence are connected totwo-time probabilities. Within our FLC approach, thehigher order probabilities ρ FLC(n > , ss are multivariate Gaus-sian densities. These belong to a quadratic static poten-tial with only one valley and no escape process. There-fore, there is no direct route to calculate the mean escapetime on the basis of the FLC approach.However, in the limit of small noise intensities, wecan find an approximation for τ K by using Kramerstheory for (Markovian) systems characterized by a staticpotential landscape U . When the potential barriers∆ U are large compared to the thermal energy, i.e.,when γD (cid:28) ∆ U , the Kramers theory provides anArrhenius formula [35, 42, 44] for the escape rate, r K = (cid:112) U (cid:48)(cid:48) ( x min ) | U (cid:48)(cid:48) ( x max ) | / (2 πγ ) exp ( − ∆ U /γD ), (with U (cid:48)(cid:48) ( x ex ) = ∂ xx U ( x ) | x = x ex , with x ex ∈ { x min , x max } ). Thisestimate for r K relies on an equilibrium approximationfor ρ and does not involve ρ . Due to the two symmetricways to leave a valley of the “washboard” potential,the resulting mean escape time can be approximated by τ K = 1 / (2 r K ) [42] and for the double-well potential withthe unique direction to exit each valley, by τ K = 1 /r K .To use the Kramers theory in the context of the FLCapproach, we note that the first member of the FPE hi-erarchy [Eq. (9a)] is formally identical to a Markovian(non-delayed) FPE for ρ , ss with a static effective po-tential V FLCeff given by Eq. (11). Thus, we can directlyapply the Arrhenius formula with U ≡ V FLCeff . For thelinear delay force in our examples, one just needs to sub- -3 -2 -1 0 log ( τ/τ B ) t h i r d m o m e n t µ x > k =2 k B T σ − k =4 k B T σ − FLCPT s.t.PT OUsmall τ BD FIG. 4. (Color online) Third moment µ x> of the steadystate probability density within the right well of the double-well potential, i.e., of ρ , ss ( x > stitute ∆ U = ∆ V + ( k/ x max − x min ) and U (cid:48)(cid:48) ( x ex ) = V (cid:48)(cid:48) s ( x ex ) + k (1 − | d i | ). Taken altogether, our estimateof τ K involves two separate approximations: the FLC toobtain V FLCeff and the application of the Arrhenius for-mula to this non-Markovian system (and implicitly, allsimplifications made within the Kramers theory).
2. Escape times in the delayed “washboard” potential
We now apply the Kramers-FLC approximation to cal-culate the escape times for the delayed “washboard” po-tential. Exemplary results for τ K are plotted in Fig. 5together with numerical data. For F d ≡
0, the BD sim-ulation results and the Kramers-FLC estimate coincide.In the regarded parameter range, τ K generally increaseswith k and τ , and the FLC results are roughly in agree-ment with the numerical results. This is rather remark-able, due to the crudeness of applying the Kramers the-ory to the non-Markovian system. At large values of τ , one observes a qualitatively different behavior: whilethe escape times resulting from the Kramers-FLC esti-mate saturate, the corresponding BD data continue toincrease with τ . This difference becomes particularlyprominent for large k . We propose that the discrep-ancy can be explained as follows. As already discussedin Sec. III B, the conditional probability saturates at afinite value τ sat , which is related to the intra-well relax-ation time τ ir . Accordingly, also V FLCeff , and therewith τ K from the Kramers-FLC approximation, saturate at τ sat .In the true non-Markovian system, the escape times notonly depend on ρ , ss , but also on higher n -time probabil-ities. Therefore they don’t necessarily saturate togetherwith ρ , ss . Indeed, our numerical investigations revealthat the saturation of τ K sets in at much higher τ values,see Fig. 5. For larger τ , another, not yet discussed timescale becomes increasingly important, that is the jump0 τ/τ B l og ( τ K / τ B ) k =0 . k B T σ − k =0 . k B T σ − BDKramers-FLC
FIG. 5. (Color online) Logarithm of the mean escape time τ K for the “washboard” as a function of the delay time τ .Cyan color: delay force amplitude k = 0 . k B T /σ , black color: k = 0 . k B T /σ . Symbols: BD, lines: Kramers-FLC estimate. duration time. In the present system, this time is roughly2 τ B (according to BD). For τ in the range of the jump du-ration times, the temporal changes of the total potentialand the significant changes of the system state variable χ occur on similar time scales. The interplay between bothmotions then leads to new (inter-well) dynamical behav-ior not captured by our approach. It hence becomes lessjustified to treat the delayed system as a quasi-Markovianone, whose stochastic (non-delayed) dynamics evolve ina static (effective) potential.One example for such new dynamical behavior arequasi-regular oscillations of χ between two valleys. Infact, we have observed analogous dynamics in the de-layed bistable system. For both static potentials, thedelay-induced oscillations have a mean period of about τ and start at random times. We further observe that theyoccasionally pause and randomly set in again. Figure 6shows BD results for the distribution of waiting timesbetween sequential jumps in the delayed double-well po-tential for an exemplary parameter setting. One sees thatmost waiting times lie within an interval ∆ t jump ∈ [0 , τ ]around a single, yet broad peak at about τ /
2. This il-lustrates the stochastic character of the oscillations withmean period of about τ . Similar quasi-regular oscilla-tions can been found in the double-well potential withlinear delay force F d ( x, x τ ) ∝ x τ as reported in [46] andfurther discussed in [17]. Interestingly, this version ofdelayed bistable system has moreover been shown to ex-hibit coherence resonance, i.e., that the regularity of theoscillations is maximal at a certain finite noise intensity.Preliminary numerical studies suggest that our modelsystems also display coherence resonance. Further in-vestigations on the delay-induced oscillatory behavior inour systems are in progress. τ τ τ τ τ τ τ ∆ t jumps P ( ∆ t j u m p s ) › ∆ t jumps fi k =8 k B T /σ , τ =5 . τ B FIG. 6. (Color online) Waiting time distribution obtainedfrom a numerically (BD) generated (normalized) histogramof time between sequential jumps ∆ t jumps , for the delayeddouble-well potential with k = 8 k B T /σ and τ = 5 . τ B . Thedashed black line marks the mean waiting time (cid:104) ∆ t jumps (cid:105) . IV. SPATIAL AUTOCORRELATIONFUNCTION
In this final section, we aim to discuss a fundamentaldifference between the FLC and the PT, which has notbeen in the focus so far. This point also provides an ex-planation why the FLC performs better on the exemplarysystems considered in the present work.We recall that both approaches yield closed FPEs for ρ via an approximation of the conditional probability ρ c , ss ( x τ , t − τ | x, t ), and therewith of the two-time prob-ability density ρ , ss ( x, t ; x τ , t − τ ) of the delayed system.The essential difference is that in the framework of theFLC, ρ c , ss is approximated by the corresponding functionof the linearized delayed system, while in the PT ap-proach, ρ c , ss stems from the corresponding non-delayedsystem. The FLC hence involves a ρ related to a non-Markovian system, while the PT utilizes a “Markovian ρ ”. On the level of the two-time probability density,the difference between a Markovian and a non-Markoviansystem manifests itself, e.g., via the spatial autocorrela-tion C ( τ ) = (cid:104) χ ( t ) χ ( t − τ ) (cid:105) ss − (cid:104) χ ( t ) (cid:105) obtained by inte-grating ρ times xx τ over both spatial variables and sub-tracting the squared first moment, see Eq. (B7). Whilein the Markovian case C ( τ ) always vanishes for large τ ,reflecting the loss of memory, it stays finite for typicalnon-Markovian systems, revealing that the system statesat the times t and t − τ are correlated for large τ . Thisis illustrated in Fig. 7, where we have plotted C ( τ ) fortwo of the considered non-Markovian systems (the “wash-board” potential and the quadratic potential received bythe force linearization, both with linear delay force), andfor the corresponding Markovian systems (obtained bysetting the delay force zero). The fact that only the FLCinvolves a non-Markovian quantity to approximate ρ c , ss ,which accounts for these finite correlations, gives an ex-1 -1 -0.5 0 0.5 1 log ( τ/τ B ) › χ ( t ) χ ( t − τ ) fi ss F d =0 F d linear F s non-linear F s FIG. 7. (Color online) Steady state spatial autocorrelationfunction C ( τ ) = (cid:104) χ ( t ) χ ( t − τ ) (cid:105) ss as a function of the delay time τ within the “washboard” potential and the force-linearizedversion of it. Symbols: BD results for the “washboard”, solidlines: exact results [Eq. (B7)] for the force linearized systems,black color: k = 8 k B T /σ , magenta color: k = 0. planation why it yields better results for large delay timesthan the perturbation theory. V. CONCLUSION
In this work, we introduce the FLC as a novel approachto close the FPE for the steady state probability densityof delayed systems. The new closure is achieved by lin-earizing the total deterministic force in all members of theFPE hierarchy starting from the second. We have spec-ified on a linear delay force, which can experimentallybe implemented, for instance, by an optical laser tweezeracting on a polarizable colloidal particle [33]. However,the novel approach generally also applies to nonlinear de-lay forces.By considering two exemplary model systems, wepresent numerical evidence that the novel approach gen-erates indeed a very good approximation of the one-timeprobability density. In fact, the FLC performs much bet-ter here compared to the small delay (Taylor) expansionintroduced in Ref. [19], and to the perturbation theoryon the level of the FPE from Ref. [20]. We also note thatthe FLC is somewhat more general than the PT in so faras it does not rely (contrary to the PT) on an analyticalexpression for the conditional probability for the corre-sponding system without delay force. And even if thisquantity is available, the respective quantity of the cor-responding linearized , delayed system ρ linc , ss (which is thekey ingredient within the FLC), often renders a betterapproximation, especially for large delay times. This isbecause ρ linc , ss is a non-Markovian quantity which can ac-count for the finite correlations between the system statesat the times t and t − τ . These correlations are not rep-resented in the PT. Moreover, the FLC is not restricted to small noise levels, as the small delay expansion, or tosmall feedback forces, as the PT approach.However, the FLC also has its limitations. It onlyworks in the steady state, in contrast to the earlier ap-proaches. Additionally, also the corresponding force lin-earized system must have steady state. Another issueoccurs when our estimate for the total potential in thesteady state V STAT [Eq. (5)] has no minima giving riseto an ambiguity regarding the points around which thelinearization should be carried out. However, one caneasily adjust the ansatz of the FLC, by simply usingother (e.g., random) centers of the linearization, prefer-ring those where the probability density is high (in orderto minimize the error in the final approximation). Forinstance, the minima of the static potential can be rea-sonable choices.Having in mind these restrictions, the FLC can beapplied to a much wider class of delayed systems thanthe ones considered in the present work. One impor-tant extension are systems with inertia, for which theFokker-Planck description has been discussed extensivelyin Ref. [18]. Such systems have been shown to exhibitcomplex dynamical behavior including chaos; an exam-ple is a delayed double-well potential with inertia andperiodic forcing [53]. Furthermore, it is in principle feasi-ble to generalize the FLC towards systems with multiplediscrete delays, which, for instance, can play a crucialrole in the stochastic process of gene regulation [11] orin the collective noisy motion of animals [54]. The force-linearized case of such systems is still solvable [21, 54],with the n -time probabilities being again given by mul-tivariate Gaussian distributions. Also for the more gen-eral class of noisy systems with distributed delays, theforce linearized case can be handled analytically for ar-bitrary memory kernels [21, 54–56]. This generally pavesthe path to adapt the basic idea of the FLC. However,to the best of our knowledge, the derivation of the FPEwith multiple or distributed delays and nonlinear forceshas not been carried out to date, at least not in thegeneral case. A lack of this equation inherently limitsthe applicability of the FLC. Another important gener-alization would be the case of multiplicative noise (anddelay), for which so far not many analytical results areknown. It is clear that the situation will substantiallydepend on the specific functional form of the noise term.While at least the delayed Fokker-Planck equation hasbeen worked out for the general case (independent of thespecific multiplicative noise term) [19, 57], there is to thebest of our knowledge no general solution for the linearforce case at hand. For some specific multiplicative noiseterms, however, the solution for the linear force modelis indeed known [29, 56, 57], making the FLC in factfeasible (with non-Gaussian conditional probability dis-tributions). Generally speaking, for all types of delayedstochastic systems, the applicability of the FLC dependson the availability of both, the general FPE and the so-lution of the force-linearized case.A further aspect touched in the present work concerns2the calculation of transport-related quantities such as es-cape times or waiting time distributions. In general, theyrequire the two-time probability ρ or even higher n -timedensities. In the present form, the FLC does not provideaccess to ρ n ≥ . Still, we have demonstrated that a rea-sonable estimate of the escape times is possible by com-bining the FLC and the Kramers theory. This estimate,however, breaks down when the delay times are in therange of (or large compared to) the jump duration times.Under such conditions, the interplay between the dynam-ics of the system state variable χ and the time-dependentenergy landscape causes oscillatory motion, whose de-scription is clearly beyond the rather crude FLC-Kramersapproximation. For this reason, an extension of the the-ory towards higher members of the FPE hierarchy wouldbe worthwhile. The analytical results concerning the sec-ond member of the FPE hierarchy and ρ for the linearcase, given in the Appendices A-C, are essential steps inthis direction.A different way to access transport-related quantitieswas presented in [46], where the autocorrelation func-tion in a very similar model (a double-well potentialsupplemented with linear delay force) was successfullyapproximated by that of an appropriate discrete model(where an entire potential valley is regarded as one dis-crete state). In this case the Master equation can analyti-cally be solved. We aim to stress that, in its present form,the FLC does not compete with this approach, but yieldsa description on a different time and length scale. WhileRef. [46] yields a convincing approximation of the inter-well dynamics, our approach rather targets the intra-welldynamics. Hence, to some extent, both approaches com-plement each other.Finally, we want to briefly comment on the possiblerelevance of the FLC for stochastic thermodynamics. In-deed the thermodynamics of delayed systems is a subjectof strong current interest [18, 27]. Of particular interestare the information and entropy change rates, which relyon two-time probabilities. Time-delayed feedback forceshave been shown to yield nontrivial contributions to thesequantities [18, 48]. In its present form, the framework ofthe FLC does not provide access to the involved prob-abilities. This is another motivation to extend our ap-proach towards higher n -time probability densities. Onecould then use the FLC predictions in order to estimatethese intriguing thermodynamic quantities far from equi-librium. ACKNOWLEDGMENTS
This work was supported by the Deutsche Forschungs-gemeinschaft through SFB 910 (project B2).
Appendix A: SECOND MEMBER OF THEDELAYED FPE, GENERAL CASE
In the following we present a derivation of the delayedFPE for ρ , which, to the best of our knowledge, has notbeen reported elsewhere so far. We extend the deriva-tion presented in Ref. [20] [on which we comment beforeEqs. (4)], towards the second member of the infinite FPEhierarchy. Accordingly, we start by considering the for-mal time derivative of ρ ( x, t ; x τ , t − τ ) = (cid:10) δ [ x − χ ( t )] δ [ x τ − χ ( t − τ )] (cid:11) , and use the delayed LE (2) to substitute thederivatives ∂χ ( t ) /∂t and ∂χ ( t − τ ) /∂t . We then introducethe functional Λ[Γ] ≡ δ ( x − χ ( t )) δ ( x τ − χ ( t − τ )), justlike in Ref. [20]. Using Novikov’s theorem [see Eq. (3)], wecan express the emerging correlations (cid:10) Λ[Γ]Γ( t ) (cid:11) by func-tional derivatives. We now use basic variational deriva-tive rules: δ Λ[Γ] δ Γ( t ) = δ Λ[Γ] δχ ( t ) δχ ( t ) δ Γ( t ) + δ Λ[Γ] δχ ( t − τ ) δχ ( t − τ ) δ Γ( t ) (cid:124) (cid:123)(cid:122) (cid:125) → , (A1)and analogously for δ Λ[Γ] /δ Γ( t − τ ). The last term inEq. (A1) vanishes, since it is noncausal: The particle po-sition at time t − τ cannot depend on the noise at the latertime t , since Γ is a random force with vanishing temporalcorrelation [51]. In the case of additive noise, the varia-tional derivatives δχ ( t ) /δ Γ( t ) and δχ ( t ) /δ Γ( t − τ ) do notexplicitly depend on Γ, so we can evaluate the ensembleaverages and obtain the FPE for ρ ∂∂t ρ ( x, t ; x τ , t − τ ) = − ∂∂x (cid:26) F ( x, x τ ) γ ρ ( x, t ; x τ , t − τ ) (cid:27) − ∂∂x τ (cid:26)(cid:90) Ω F ( x τ , x τ ) γ ρ ( x, t ; x τ , t − τ ; x τ , t − τ )d x τ (cid:27) + √ D (cid:32)(cid:20) ∂ ∂x + ∂ ∂x τ (cid:21)(cid:26) δχ ( t ) δ Γ( t ) (cid:12)(cid:12)(cid:12) χ ( t )= xχ ( t − τ )= xτ ρ ( x, t ; x τ , t − τ ) (cid:27) + ∂ ∂x τ ∂x (cid:26) δχ ( t ) δ Γ( t − τ ) (cid:12)(cid:12)(cid:12) χ ( t )= xχ ( t − τ )= xτ ρ ( x, t ; x τ , t − τ ) (cid:27)(cid:33) . (A2)This is the second member of the Fokker-Planck hier-archy for the general case of a system with delay force F ( x, x τ ). As expected, it contains the next higher-orderprobability distribution function, that is ρ . We also notethat, in contrast to the first member of the FPE hierarchy[see Eq. (4)], Eq. (A2) contains the variational derivatives δχ ( t ) /δ Γ( t ) and δχ ( t ) /δ Γ( t − τ ). These need to be spec-ified for the specific F under consideration. Appendix B: SECOND MEMBER OF FPE FORLINEAR FORCES
Now we specialize the delayed FPE for ρ [see Eq. (A2)]for forces that are linear in x and x τ , i.e., forces of the3form F L ( x, x τ ) = − αx − βx τ . The corresponding delayedLE can be solved iteratively by using the method of steps[51]: to solve the equation in a time interval t ∈ [ nτ, ( n +1) τ ] , ∀ n ∈ N , one substitutes the solution of the precedinginterval [( n − τ, nτ ] to get rid of the delay term (for n = 0 the history function is used). One finds that oneach interval, the formal solution has the same structure Y n +1 ( t ) = Y n ( nτ ) e − α ( t − nτ ) / ( γσ ) + (cid:90) tnτ e − α ( t − t (cid:48) ) / ( γσ ) [ (cid:112) D Γ( t (cid:48) ) − β Y n ( t (cid:48) − τ )]d t (cid:48) , (B1)with χ ( t ) = Y n +1 ( t ), for t ∈ [ nτ, ( n +1) τ ]. From Eq. (B1)we find the functional derivatives δχ ( t ) δ Γ( t ) = δχ ( t − τ ) δ Γ( t − τ ) = (cid:114) D , (B2) δχ ( t ) δ Γ( t − τ ) = (cid:112) D e − ατ/ ( γσ ) . (B3)Inserting Eqs. (B2) and (B3) into Eq. (A2) we obtain forthe case of linear forces F L ( x, x τ ) = − αx − βx τ : ∂∂t ρ ( x, t ; x τ , t − τ ) = − ∂∂x (cid:26) F L ( x, x τ ) γ ρ ( x, t ; x τ , t − τ ) (cid:27) − ∂∂x τ (cid:90) Ω F L ( x τ , x τ ) γ ρ ( x, t ; x τ , t − τ ; x τ , t − τ ) d x τ + D (cid:20) ∂ ∂x + ∂ ∂x τ (cid:21) ρ ( x, t ; x τ , t − τ )+ 2 D e − ατ/ ( γσ ) ∂ ∂x∂x τ ρ ( x, t ; x τ , t − τ ) . (B4)In the steady state, Eqs. (B4) and the FPE for ρ [Eq. (4)]are solved by the respective (multivariate) Gaussian dis-tributions [28] ρ n , ss = 1 (cid:112) (2 π ) n det { D n } e − (1 / x −(cid:104) x (cid:105) ) D n − ( x −(cid:104) x (cid:105) ) , (B5)with x = ( x, x τ , .., x ( n − τ ) T , D = C (0), and D = (cid:18) C (0) C ( τ ) C ( τ ) C (0) (cid:19) , D = C (0) C ( τ ) C (2 τ ) C ( τ ) C (0) C ( τ ) C (2 τ ) C ( τ ) C (0) . (B6) The covariance matrices D n involve the steady state spa-tial autocorrelation function between the systems state attime t and time t + z : C ( z ) := (cid:104) χ ( t ) χ ( t + z ) (cid:105) ss − (cid:104) χ ( t ) (cid:105) . (B7)Explicit expressions for C (0), C ( τ ), and C (2 τ ) (whichoccur in D k , k= 1 , ,
3) are given in the next section.
Appendix C: SPATIAL AUTOCORRELATIONFUNCTION
Due to the symmetry property C ( z ) = C ( − z ), we onlyconsider non-negative arguments z in the following. Asshown in Refs. [19, 51], the spatial autocorrelation func-tion on the interval z ∈ [0 , τ ] is given by C (0) = D γ + ( β/ω ) sinh( ωτ ) α + β cosh( ωτ ) , (C1) C ( z ) = C (0) cosh( ωz ) − ( D /ω ) sinh( ω | z | ) (C2)with ω = (cid:112) α − β /γ ∈ C . In order to obtain C ( z ) for z >
0, one can directly deduce from the delayed LE [51]:d C ( z )d z = (cid:42) χ ( t ) d χ ( u )d u (cid:12)(cid:12)(cid:12) u = t + z (cid:43) ss = − αγ C ( z ) − βγ C ( z − τ ) + (cid:112) D (cid:104) χ ( t )Γ( t + z ) (cid:105) . (C3)For z >
0, the last term vanishes, since a non-zero cor-relation between the system state χ and the future noisewould be noncausal. The resulting differential equationfor z > C ( z ) / d z = − ( α/γ ) C ( z ) − ( β/γ ) C ( z − τ ) , (C4)can be solved iteratively using the method of steps (seetext at the beginning of the Appendix B). Using the so-lution for the interval z ∈ [0 , τ ] in Eqs. (C1) and (C2), weobtain the correlation function on the interval z ∈ [ τ, τ ], C ( z ) = D e α ( τ −| z | ) [ β cosh( ωτ ) + α ] + [ ... ] ∗ β [ α + β cosh( ωτ )] , (C5)where[ ... ] ∗ = − α cosh[ ω ( τ − z )] − β cosh[ ω (2 τ − z )]+ [( β − α ) /ω ] sinh[ ω ( τ − | z | )] − ( αβ/ω ) sinh[ ω (2 τ − | z | )] . (C6)Inspection of Eq. (C5) shows that there exist critical τ -values for which C ( z ) diverges, see also text after Eq. (8d)(with α = α i and β = β i in the discussion there).We would like to note that in Ref. [51], the solution(C2) is erroneously stated to be valid for all z ∈ R . Alongtheir derivation, Eq. (18) of Ref. [51] is actually not validfor | z | > τ , although they claim otherwise. A correlationfunction (C2) ∀ z ∈ R would in fact be unphysical. Con-sider for example β < α , where the system is spatiallyconfined. Then, it is not reasonable that the autocorrela-tion function grows without bounds for large time differ-ences, as Eq. (C2) would predict. The iterative solutionobtained from Eq. (C4) by using the method of steps, isindeed bounded, consisting with the physical intuition.4 Appendix D: BROWNIAN DYNAMICSSIMULATIONS
In order to test our theoretical results against quasi-exact data, we perform Brownian dynamics (BD) simula-tions of the delayed LE (2). We use the Euler-Maruyamaintegration scheme [58, 59], with a varying temporal dis-cretization of ∆ t ∈ [10 − , − ] τ B , such that 1000 ∆ t ≤ τ .We perform each simulation multiple times (with differ-ent random number seed), and take the average over atleast 10 realizations. The (pseudo) random numbers aregenerated with the algorithm “Mersenne Twister” [60]and the Box-M¨uller method [59]. As initial condition, we use an equilibrium configura-tion in the respective static potential without delay force.Before we measure the steady state properties, we let100 τ B times pass in the presence of the delay force tolet the system reach a steady state. The simulation timethereafter is more than 200 τ B .The density profiles from the BD simulations are ob-tained from histograms with a spatial resolution of ∆ x =0 . σ . These histograms are also used to calculate themoments µ n . With the same procedure the momentsof the analytical distributions resulting from the FLC,PT and small delay expansion are calculated (we use thesame bin sizes and positions as in the simulations). [1] A. Longtin, Complex time-delay systems: theory and ap-plications , edited by F. M. Atay (Springer-Verlag BerlinHeidelberg, 2010) pp. 177–195.[2] E. Sch¨oll, S. H. L. Klapp, and P. H¨ovel, eds.,
Control ofself-organizing nonlinear systems (Springer, 2016).[3] E. Sch¨oll and H. G. Schuster, eds.,
Handbook of chaoscontrol (John Wiley & Sons, 2008).[4] M. Montinaro, A. Mehlin, H. Solanki, P. Peddibhotla,S. Mack, D. Awschalom, and M. Poggio, Appl. Phys.Lett. , 133104 (2012).[5] G. Stoica, Proc. Am. Math. Soc. , 1837 (2005).[6] H. U. Voss and J. Kurths,
Modelling and forecasting fi-nancial data: techniques of nonlinear dynamics , editedby A. S. Soofi and L. Cao, Vol. 2 (Springer Science &Business Media, 2002) pp. 327–349.[7] A. Longtin, J. G. Milton, J. E. Bos, and M. C. Mackey,Phys. Rev. A , 6992 (1990).[8] J. G. Milton, A. Fuerte, C. B´elair, J. Lippai,A. Kamimura, and T. Ohira, Nonlinear Theory andIts Applications, IEICE , 129 (2013); J. G. Milton,T. Ohira, J. L. Cabrera, R. M. Fraiser, J. B. Gyorffy,F. K. Ruiz, M. A. Strauss, E. C. Balch, P. J. Marin, andJ. L. Alexander, PLoS ONE , e7427 (2009).[9] N. S. Goel, S. C. Maitra, and E. W. Montroll, Rev. Mod.Phys. , 231 (1971).[10] K. Das, M. Srinivas, M. Srinivas, and N. Gazi, C. R.Biol. , 503–513 (2012).[11] K. Parmar, K. B. Blyuss, Y. N. Kyrychko, and S. J.Hogan, Comput. Math. Meth. M. , 347273 (2015);D. Bratsun, D. Volfson, L. S. Tsimring, and J. Hasty,Proc. Natl. Acad. Sci. U. S. A. , 14593 (2005);C. Gupta, J. M. L´opez, W. Ott, K. Josi´c, and M. R.Bennett, Phys. Rev. Lett. , 058104 (2013).[12] C. Masoller, Phys. Rev. Lett. , 034102 (2002).[13] S. M. Hein, F. Schulze, A. Carmele, and A. Knorr, Phys.Rev. A , 052321 (2015); Y. Lu, N. L. Naumann, J. Cer-rillo, Q. Zhao, A. Knorr, and A. Carmele, arXiv preprintarXiv:1703.10028 (2017).[14] M. Braun and F. Cichos, ACS Nano , 11200–11208(2013); M. Braun, A. P. Bregulla, K. G¨unther, M. Mer-tig, and F. Cichos, Nano Lett. , 5499–5505 (2015);D. F. B. Haeufle, T. Bauerle, J. Steiner, L. Bremicker,S. Schmitt, and C. Bechinger, Phys. Rev. E , 012617(2016).[15] K. Lichtner, A. Pototsky, and S. H. L. Klapp, Phys. Rev. E , 051405 (2012); K. Lichtner and S. H. L. Klapp,EPL , 40007 (2010).[16] D. Hennig, Phys. Rev. E , 041114 (2009).[17] C. Masoller, Phys. Rev. Lett. , 020601 (2003).[18] M. L. Rosinberg, T. Munakata, and G. Tarjus, Phys.Rev. E , 042114 (2015).[19] S. Guillouzic, I. L’Heureux, and A. Longtin, Phys. Rev.E , 3970 (1999).[20] T. D. Frank, Phys. Rev. E , 031106 (2005); Phys. Rev.E , 011112 (2005).[21] L. Giuggioli, T. J. McKetterick, V. M. Kenkre, andM. Chase, J. Phys. A , 384002 (2016).[22] U. Seifert, Rep. Prog. Phys. , 126001 (2012).[23] T. Sagawa and M. Ueda, Phys. Rev. E , 021104 (2012);D. Abreu and U. Seifert, Phys. Rev. Lett. , 030601(2012); A. C. Barato and U. Seifert, , 090601 (2014);D. Mandal and C. Jarzynski, Proc. Natl. Acad. Sci.U.S.A. , 11641 (2012); J. M. R. Parrondo, J. M.Horowitz, and T. Sagawa, Nat. Phys. , 131–139(2015); A. Kutvonen, T. Sagawa, and T. Ala-Nissila,Phys. Rev. E , 032147 (2016).[24] S. A. M. Loos, R. Gernert, and S. H. L. Klapp, Phys.Rev. E , 052136 (2014).[25] K. H. Kim and H. Qian, Phys. Rev. E , 022102 (2007).[26] J. V. Koski, V. F. Maisi, T. Sagawa, and J. P. Pekola,Phys. Rev. Lett. , 030601 (2014).[27] H. Jiang, T. Xiao, and Z. Hou, Phys. Rev. E , 061144(2011); T. Munakata, S. Iwama, and M. Kimizuka, Phys.Rev. E , 031104 (2009); T. Munakata and M. L. Ros-inberg, Phys. Rev. Lett. , 180601 (2014).[28] U. K¨uchler and B. Mensch, Stoch. Stoch. Rep. , 23–42(1992).[29] T. D. Frank and P. J. Beek, Phys. Rev. E , 021917(2001).[30] S. Guillouzic, I. L’Heureux, and A. Longtin, Phys. Rev.E , 4906 (2000).[31] U. M. B. Marconi and P. Tarazona, J. Chem. Phys. ,8032 (1999); A. M. Menzel, A. Saha, C. Hoell, andH. L¨owen, J. Chem. Phys. , 024115 (2016).[32] S. P. Das, Rev. Mod. Phys. , 785 (2004); L. M. C.Janssen and D. R. Reichman, Phys. Rev. Lett. ,205701 (2015).[33] J. Kotar, M. Leoni, B. Bassetti, M. C. Lagomarsino, andP. Cicuta, Proc. Natl. Acad. Sci. U.S.A. , 7669 (2010).[34] B. Qian, D. Montiel, A. Bregulla, F. Cichos, and H. Yang, Chem. Sci. , 1420 (2013); A. Balijepalli, J. J.Gorman, S. K. Gupta, and T. W. LeBrun, Nano Lett. , 2347 (2012).[35] R. Gernert, C. Emary, and S. H. L. Klapp, Phys. Rev.E , 062115 (2014).[36] R. Gernert and S. H. L. Klapp, Phys. Rev. E , 022132(2015).[37] M. W¨ordemann, Structured Light Fields: Applicationsin Optical Trapping, Manipulation, and Organisation (Springer Science & Business Media, 2012).[38] P. Reimann, C. Van den Broeck, H. Linke, P. H¨anggi,J. M. Rubi, and A. P´erez-Madrid, Phys. Rev. Lett. ,010602 (2001).[39] M. P. N. Juniper, A. V. Straube, D. G. A. L. Aarts, andR. P. A. Dullens, Phys. Rev. E , 012608 (2016).[40] C. Emary, R. Gernert, and S. H. L. Klapp, Phys. Rev.E , 061135 (2012).[41] H. A. Kramers, Physica , 284 (1940).[42] P. H¨anggi, P. Talkner, and M. Borkovec, Rev. Mod.Phys. , 251 (1990).[43] H. Risken, The Fokker-Planck Equation (Springer, 1984).[44] C. W. Gardiner,
Handbook of Stochastic Methods , 2nded., 6th printing (Springer, Berlin–Heidelberg–New York,2002).[45] D. Goulding, S. Melnik, D. Curtin, T. Piwonski, J. Houli-han, J. P. Gleeson, and G. Huyet, Phys. Rev. E , 031128 (2007).[46] L. S. Tsimring and A. Pikovsky, Phys. Rev. Lett. ,250602 (2001).[47] L. Du and D. Mei, Indian J. Phys. , 267 (2015).[48] T. Xiao, Phys. Rev. E , 052109 (2016).[49] T. Piwonski, J. Houlihan, T. Busch, and G. Huyet, Phys.Rev. Lett. , 040601 (2005).[50] E. A. Novikov, Sov. Phys. JETP , 1290 (1965).[51] T. D. Frank, P. J. Beek, and R. Friedrich, Phys. Rev. E , 021912 (2003).[52] K. Pyragas, Phys. Lett. A , 421 (1992).[53] B. N. Nbendjo, R. Tchoukuegno, and P. Woafo, Chaos,Solitons & Fractals , 345 (2003).[54] T. J. McKetterick and L. Giuggioli, Phys. Rev. E ,042135 (2014).[55] P. H¨anggi and P. Talkner, Phys. Lett. A , 9 (1978);P. H¨anggi and H. Thomas, Phys. Rep. , 207 (1982);A. Hern´andez-Machado, J. Sancho, M. San Miguel, andL. Pesquera, EPJ B , 335 (1983).[56] A. A. Budini and M. O. C´aceres, J. Phys. A , 5959(2004).[57] T. D. Frank, Phys. Rev. E , 061104 (2004).[58] E. Buckwar, J. Comput. Appl. Math. , 297 (2000).[59] P. Kloeden and E. Platen, Numerical solution of stochas-tic differential equations (Springer–Verlag, 1992).[60] (June 2016), weblink: