Force-linearization closure for non-Markovian Langevin systems with time delay

TL;DR

This paper introduces a linearization-based closure method for non-Markovian Langevin systems with delay, providing accurate steady-state probability densities and improving upon existing approximations.

Contribution

A novel linearization closure approach for the Fokker-Planck hierarchy in delayed Langevin systems, enabling accurate steady-state density calculations.

Findings

Accurately predicts steady-state densities for delayed nonlinear systems.

Outperforms small-delay and perturbation approximations.

Applicable to a wide class of nonlinear force models.

Abstract

This paper is concerned with the Fokker-Planck (FP) description of classical stochastic systems with discrete time delay. The non-Markovian character of the corresponding Langevin dynamics naturally leads to a coupled infinite hierarchy of FP equations for the various $n$-time joint distribution functions. Here we present a novel approach to close the hierarchy at the one-time level based on a linearization of the deterministic forces in all members of the hierarchy starting from the second 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 potential supplemented by a linear delay force, we demonstrate that our approach yields a very accurate representation of the density as compared to quasi-exact numerical results from direct solution of the Langevin equation. Moreover, the results are significantly improved against those from a small-delay approximation and a perturbation-theoretical approach. We also discuss the possibility of accessing transport-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.

Figures

Equations

$$F(x,x_{\tau })=F_{\mathrm{s}}(x)+F_{\mathrm{d}}(x,x_{\tau }),$$

$$\mathrm{d}\chi (t)/\mathrm{d}t={\gamma }^{-1}F\left[\chi (t),\chi (t-\tau )\right]+\sqrt{2D_0}\Gamma (t),$$

$$\big{\langle}\Lambda[\Gamma]\Gamma(t)\big{\rangle}=\Bigg{\langle}\frac{\delta\Lambda[\Gamma]}{\delta\Gamma(t)}\Bigg{\rangle}.$$

$${\partial }_t{\rho }_1(x,t)=$$

$$+D_0{\partial }_{xx}{\rho }_1(x,t),$$

$$\hat{F}(x,t)=$$

$$V_{\mathrm{STAT}}(x)=-\int F(x,x_{\tau }=x)\mathrm{d}x.$$

$$F^{\mathrm{lin},i}(x,x_{\tau })=-{\alpha }^i\Delta x^i-{\beta }^i\Delta x_{\tau }^i,$$

$${\rho }_{\mathrm{c},\mathrm{ss}}^{\mathrm{lin},i}(x_{\tau }|x;\tau )=\sqrt{\frac{2}{\pi }{K}^i{d_1^i}^2\left(1-{d_2^i}^2\right)}\exp \left[\frac{{\Delta x^i}^2}{2K^i}\right]\times$$

$$\exp \left[d_1^i\left(2\Delta x^i\Delta x_{\tau }^i{d}_2^i-{\Delta x^i}^2-{\Delta x_{\tau }^i}^2\right)\right]$$

$${\omega }^i\equiv$$

$$K^i\equiv$$

$$d_1^i\equiv$$

$$d_2^i\equiv$$

$$\gamma D_0{\partial }_{xx}{\rho }_{1,\mathrm{ss}}^{\mathrm{FLC}}(x)={\partial }_x[{\hat{F}}^{\mathrm{FLC}}(x){\rho }_{1,\mathrm{ss}}^{\mathrm{FLC}}(x)],$$

$${\hat{F}}^{\mathrm{FLC}}(x)={\int }_{\Omega }F(x,x_{\tau }){\rho }_{\mathrm{c},\mathrm{ss}}^{\mathrm{lin},i}(x_{\tau }|x;\tau )\mathrm{d}{x}_{\tau },\forall x\in {\Omega }^i.$$

$${\rho }_{1,\mathrm{ss}}^{\mathrm{FLC}}(x)=Z^{-1}\exp \left[-V_{\mathrm{eff}}^{\mathrm{FLC}}(x)/(k_{\mathrm{B}}\mathcal{T})\right]$$

$$V_{\mathrm{eff}}^{\mathrm{FLC}}(x)=-{\int }_{\hat{x}}^x\mathrm{d}{x}^{\prime }{\hat{F}}^{\mathrm{FLC}}(x^{\prime }),$$

$$F_{\mathrm{d}}(x,x_{\tau })=-k(x-x_{\tau }),$$

$${\hat{F}}^{\mathrm{FLC}}(x)=-{\partial }_x{V}_{\mathrm{s}}(x)-k\Delta x^i+kx{d}_2^i\sqrt{2K^i{d}_1^i}\times$$

$$\sqrt{1-{d_2^i}^2}\exp \left\{\left[d_1^i{d_2^i}^2-d_1^i+1/(2K^i)\right]{\Delta x^i}^2\right\}.$$

$$d_1^i{d_2^i}^2-d_1^i+1/(2K^i)=0,$$

$$\displaystyle V_{\mathrm{eff}}^{\mathrm{FLC}}(x)=V_{\mathrm{s}}(x)+{\color[rgb]{0,0,0}{k}/{2}\!\left(1-|d^{i}_{2}|\right)}{\Delta x^{i}}^{2},$$

$${\rho }_{1,\mathrm{ss}}^{\mathrm{approx}}(x)=Z^{-1}\exp \left[-V_{\mathrm{eff}}^{\mathrm{approx}}(x)/(k_{\mathrm{B}}\mathcal{T})\right].$$

$$V_{\mathrm{eff}}^{\mathrm{small}\tau }(x)=V_{\mathrm{s}}(x)/(1-k\tau ).$$

$${\rho }_{\mathrm{c}}(x_{\tau }$$

$$\frac{1}{\sqrt{4\pi D_0\tau }}\exp \left[-\frac{{\{x_{\tau }-x-F_{\mathrm{s}}(x)\tau /\gamma \}}^2}{4D_0\tau }\right],$$

$$V_{\mathrm{eff}}^{\mathrm{PT}}(x)=(1+k\tau )V_{\mathrm{s}}(x).$$

$$V_{\mathrm{eff}}^{\mathrm{PT}}(x)=$$

$$\left\{1-\exp [-({\alpha }^i-k)\tau /\gamma ]\right\}{\Delta x^i}^2,\forall x\in {\Omega }^i.$$

$$V_{\mathrm{s}}(x)=-(\Delta V_{}/2)\cos (x/\sigma )$$

$$d_2^i=\frac{\gamma \omega \cosh (\tau \omega )-[\Delta V_{}/(2{\sigma }^2)+k]\sinh (\tau \omega )}{\gamma \omega -k\sinh (\tau \omega )}$$

$$V_{\mathrm{s}}(x)=\Delta V_{}[{\left(x/\sigma \right)}^4-2{\left(x/\sigma \right)}^2],$$

$$d_2^i=$$

$$\frac{\delta \Lambda [\Gamma ]}{\delta \Gamma (t)}=\frac{\delta \Lambda [\Gamma ]}{\delta \chi (t)}\frac{\delta \chi (t)}{\delta \Gamma (t)}+\frac{\delta \Lambda [\Gamma ]}{\delta \chi (t-\tau )}\underset{\to 0}{\underbrace{\frac{\delta \chi (t-\tau )}{\delta \Gamma (t)}}},$$

Full text

I Introduction

Many (non-)equilibrium systems from the macroscopic world down to the quantum level are governed by dynamical equations which involve both, noise and delay in time Longtin (2010); Schöll et al. (2016); Schöll and Schuster (2008). Noise due to imperfections or environmental influences is essentially omnipresent in many real-world systems and experimental setups (e.g., at the cantilever of an Atomic Force Microscope Montinaro et al. (2012)). On the mathematical level it often also results from the presence of hidden degrees of freedom, which have disappeared in the context of a coarse-graining procedure. Similarly, time delay can have various sources. One example are finite information processing and reaction times, which occur for example in social systems (e.g., financial markets Stoica (2005) or economic processes Voss and Kurths (2002)). Such latencies as well appear in neural systems due to finite signal transmission and refractory times (which, e.g., become apparent in the human pupil reflex Longtin et al. (1990), or in stick balancing experiments Milton et al. (2013); *Milton2009a). In biological systems, delays can be caused by maturation times (as in population growth Goel et al. (1971) or prey-predator dynamics Das et al. (2012)), or, on the biomolecular level, by biochemical reactions kinetics (such as the transcriptional and translational delay in gene networks Parmar et al. (2015); *Bratsun2005; *Gupta2013). In laser dynamics with optical feedback, delay is induced by the traveling time of laser light on its round trip in the cavity Masoller (2002). On an even smaller scale, time delays occur in quantum-optical systems coupled to a structured photonic reservoir Hein et al. (2015); *Lu2017.

Beyond these intrinsic delays, time delay is an important issue in experimental setups with feedback control, where the delay is due to the finite time to proceed the output (plus the signal transmission times to obtain the information and feed it back into the system). Indeed, feedback control has emerged as an important tool to manipulate small systems. Important representatives for such systems are colloids in a thermal bath under time delayed feedback control Braun and Cichos (2013); *Braun2015; *Haeufle2016; Lichtner et al. (2012); *Lichtner2010; Hennig (2009). All these examples involve noise and delay, illustrating that due to the generality and omnipresence of such features, noisy systems with time delays occur on all scales. In many of these cases, the interplay of noise and delay moreover leads to intriguing dynamical behavior such as multistability and stochastic switches Masoller (2002, 2003).

However, even in the classical case the mathematical description of noisy systems with time delay continues to be a major challenge. In presence of a discrete time delay $\tau$, the standard mathematical equation for stochastic motion, that is, the (under- or overdamped) Langevin equation (LE), becomes a stochastic delay differential equation (SDDE), or, in a more physical language, a non-Markovian Langevin equation. The usual concepts and solutions for Markovian LEs then do not apply and indeed, the development of a general solution method for SDDEs is yet an unresolved problem Rosinberg et al. (2015). Moreover, contrary to the Markovian systems, the route from the LE to a Fokker-Planck equation (FPE) for the corresponding probability density is more involved Rosinberg et al. (2015); Guillouzic et al. (1999); Frank (2005a); *Frank2005a; Giuggioli et al. (2016). The time delay alone leads to an infinite hierarchy of coupled FP equations for the $n$-time (joint) probability densities Rosinberg et al. (2015); Guillouzic et al. (1999); Frank (2005a); *Frank2005a; Giuggioli et al. (2016), for which as well no general solution has been found to date. Similar hierarchy problems also occur in other contexts in statistical physics. A well known example is the Bogoliubov–Born–Green–Kirkwood–Yvon (BBGKY) hierarchy, in which the time evolution of the one-particle density depends upon the two-particle density and so on. For this problem, various closing strategies have been proposed, such as the simple mean-field (factorization) approximation and the more sophisticated dynamical density functional theory (involving an adiabatic approximation) Marconi and Tarazona (1999); *Menzel2016. Another closure example is the mode coupling theory for glassy systems Das (2004); *Janssen2015. We stress, however, that in all these examples the hierarchical structure emerges due to (conservative) particle-particle interactions. In contrast, the hierarchy of FPEs appearing in time-delayed systems – although having a similar structure – arises due to delay-induced temporal correlations.

These aspects are also of major relevance in the emerging field of stochastic thermodynamics Seifert (2012); Sagawa and Ueda (2012); *Abreu2012; *Barato2014; *Mandal2012; *Parrondo2015; *Kutvonen2016; Loos et al. (2014). Indeed, concepts like entropy production, fluctuation theorems and information exchange have been widely formulated (and tested experimentally Kim and Qian (2007); Koski et al. (2014)) for Markovian noisy systems, for which the correspondence between the various levels of description (LE, FPE, path integrals) can readily be utilized. For example, the entropy production from the system can be directly calculated via the probability density and the corresponding probability current Seifert (2012). For non-Markovian systems, all these concepts have to be revisited (see Rosinberg et al. (2015); Jiang et al. (2011); *Munakata2009; *Munakata2014 for recent work in this direction).

The only solvable class of delayed noisy systems are those with linear forces, where all the $n$-time probability densities are given by multivariate Gaussian distributions. Here, the system can be solved exactly on the LE Küchler and Mensch (1992) and the FPE level Giuggioli et al. (2016), or by utilizing both levels of description Frank and Beek (2001). In the more general case of nonlinear systems, one has to rely on approximations Longtin (2010). For continuous systems, the two most established strategies are a Taylor expansion of the LE in $\tau$ Guillouzic et al. (1999, 2000) and a perturbation theory on the level of the FPE Frank (2005a); *Frank2005a, where the entire delay force is treated as a small deviation from the Markovian dynamics.

In this work, we propose a novel FPE-based approach, which we call “force linearization closure” (FLC). The main idea is to linearize the deterministic force in all members of the FPE hierarchy starting from the second one. We then utilize the (Gaussian) solution for the linear, delayed higher-order FPEs to obtain a closure of the FPE for the one-time probability density in the nonlinear case. Our strategy can be applied whenever the linearized system has a stable steady state.

We apply this concept to a classical system, namely, an (overdamped) colloidal particle subject to a nonlinear static force and a time-delayed force representing a feedback control. The control “target” is the particle position which, due to the typical (micron-scale) size of a colloid, is readily accessible in real space experiments and simulations Braun and Cichos (2013); *Braun2015; *Haeufle2016; Lichtner et al. (2012); *Lichtner2010; Hennig (2009). Feedback control can be implemented, e.g., as a co-moving “optical tweezer” Loos et al. (2014); Kotar et al. (2010); Kim and Qian (2007); Qian et al. (2013); *Balijepalli2012; Gernert et al. (2014); Gernert and Klapp (2015), which can be well approximated by a co-moving quadratic potential Wördemann (2012); Kotar et al. (2010); Kim and Qian (2007). Depending on the experimental setup, a time delay naturally arises during the position measurement and the adjustment of the tweezer, or it can be intentionally implemented as a feature. As a result, the system is subject to a linear time-delayed feedback force.

Regarding the nonlinear static force (or potential, respectively) we consider two paradigmatic examples: first, a periodic sinusoidal (“washboard”) potential, which is frequently used to model diffusion and transport on rough surfaces Reimann et al. (2001); Juniper et al. (2016); Emary et al. (2012) and can be realized experimentally e.g. by optical landscapes Juniper et al. (2016). The impact of time delay on diffusion in washboard potentials has already been studied on the basis of a heuristic delayed FPE Loos et al. (2014); Gernert et al. (2014); Gernert and Klapp (2015) and by simulations Hennig (2009). Our second example is a colloidal particle in a double-well potential. Bistable noisy systems often serve to study escape problems Kramers (1940); Hänggi et al. (1990); Risken (1984); Gardiner (2002). Bistable systems with delay have been studied theoretically, e.g., regarding the Kramers rate Guillouzic et al. (2000); Goulding et al. (2007); Masoller (2003) and in the context of coherence resonance (see the study of Tsimring and Pikovsky Tsimring and Pikovsky (2001), who provided a solution based on a discretization procedure, and subsequent work Tsimring and Pikovsky (2001); Masoller (2003); Du and Mei (2015); Xiao (2016)). Experimental realizations have been suggested in Refs. Goulding et al. (2007); Piwonski et al. (2005), which involve the polarization dynamics of lasers with optical feedback (where $\tau$ is associated with the external cavity length).

The remainder of this paper is organized as follows. We start by briefly reviewing the Langevin and the Fokker-Planck description of noisy overdamped systems with discrete time delay, as well as the two main earlier approximations for the one-time probability density $\rho_{1}$ (Sec. II). We then introduce the force linearization closure and provide the main steps to derive its key equation, a closed FPE for $\rho_{1}$ (Sec. II.3). In the second part of the paper, we apply the FLC to the two examples mentioned above and compare the results with quasi-exact numerical results obtained from Brownian dynamics (BD) simulations, i.e., direct simulations of the delayed Langevin equation. We demonstrate that the FLC renders a very good approximation of the one-time probability density, which is, moreover, more accurate than the predictions from the two main earlier approaches. Furthermore, we briefly discuss the possibility to estimate escape times (which are rigorously connected to two-time probability densities) within the novel approach (Sec. III.4). We summarize and conclude in Sec. V. The paper contains four appendices with additional theoretical results (e.g., the second member of the FPE hierarchy), and technical aspects of the calculations.

II Theoretical framework

II.1 Delayed Langevin and Fokker-Planck equation

We consider an overdamped Brownian particle at the time-dependent stochastic position $\chi(t)$ which represents the dynamical quantity of interest. The particle moves in a static external potential $V_{\mathrm{s}}(x)$, with spatial variable $x\!\in\Omega$, where $\Omega$ denotes the spatial domain of the system. The particle is further subject to a “delay force” $F_{\mathrm{d}}$, which depends on the instantaneous and on the delayed particles position, i.e., at a given time $t$ on $\chi(t)$ and $\chi(t-\tau)$, with $\tau\!>\!0$ being the delay time. The total deterministic force is hence

[TABLE]

with $x_{\tau}\!\in\Omega$ being a second spatial variable needed due to the two involved particle positions in $F_{\mathrm{d}}$, and $F_{\mathrm{s}}(x)\!=\!-\partial_{x}V_{\mathrm{s}}(x)$ being the negative of the spatial derivative of $V_{\mathrm{s}}$. The delayed Langevin equation (LE) reads

[TABLE]

where $\gamma$ is the friction coefficient of the surrounding medium at temperature ${\cal T}$ and $\Gamma(t)$ denotes the Langevin force. The latter introduces additive Gaussian white noise to the system, i.e., $\langle\Gamma(t)\rangle\!=\!0~{}\forall t$ and $\langle\Gamma(t)\Gamma(t^{\prime})\rangle\!=\!\delta(t-t^{\prime})~{}\forall t,t^{\prime}$, with $\langle...\rangle$ denoting the ensemble average and $\delta(t)$ the Delta distribution. Further, $D_{0}$ is the short-time diffusion coefficient satisfying $D_{0}\!=\!k_{\mathrm{B}}{\cal T}/\gamma$ (with $k_{\mathrm{B}}$ being the Boltzmann constant) Risken (1984). We measure the time in units of the “Brownian” time scale $\tau_{\mathrm{B}}\!=\!\sigma^{2}/D_{0}$, where $\sigma$ is the particle diameter. This is the time in which a free Brownian particle (no external forces) travels over a distance equal to its own size.

As usual for delay equations, the evolution of the dynamical variable (here $\chi$) depends on a given history function $\phi$, which serves as an initial condition $\chi(t)\!=\!\phi(t),~{}\forall t\in[-\tau,0]$.

We emphasize that the theoretical framework discussed in this paper does not just apply to Brownian particles, but to other natural and artificial systems as well. Having this in mind, the friction and diffusion coefficient, the Brownian time scale and the particle diameter introduced before just provide energy, time and spatial length scales, which one can adjust to the system under consideration.

As firstly shown in Guillouzic et al. (1999), it is possible to derive from the delayed LE a Fokker-Planck equation (FPE) for the (one-time) probability density $\rho_{1}(x,t)\!=\!\langle\delta(x-\chi(t))\rangle$ in a similar manner to the Markovian case. We particularly refer to Ref. Frank (2005a); *Frank2005a, which presents a derivation based upon Novikov’s theorem Novikov (1965). Novikov’s theorem generally links the variational derivative of an arbitrary functional of the Langevin force, $\Lambda[\Gamma]$, to the correlations between that functional and the Langevin force. For Gaussian white noise it reads Novikov (1965); Frank (2005a):

[TABLE]

A central step in the derivation in Frank (2005a) is the usage of Novikov’s theorem to express the correlations \big{\langle}\left\{\delta(x-\chi(t))\delta(x_{\tau}-\chi(t-\tau))\right\}\Gamma(t)\big{\rangle} by the related variational derivatives. For the case of additive noise, the delayed FPE for $\rho_{1}$ is given by Guillouzic et al. (1999); Frank and Beek (2001); Frank (2005a)

[TABLE]

For the sake of a shorter notation, we drop here and in the following the dependency on the history function, i.e., $\rho_{1}(x,t)\equiv\rho_{1}(x,t|\phi)$ and $\rho_{\mathrm{c}}(x_{\tau},t\!-\!\tau|x,t)\equiv\rho_{\mathrm{c}}(x_{\tau},t\!-\!\tau|x,t;\phi)$. Comparing Eq. (4) to an ordinary FPE for Markovian systems, one notes the appearance of a delay-averaged drift term $\widehat{F}$ instead of the usual drift term, which only depends on instantaneous quantities. Specifically, the delay-averaged drift involves the conditional probability $\rho_{\mathrm{c}}$, which is related to the two-time (joint) probability density via $\rho_{2}(x,t;x_{\tau},t\!-\!\tau)\!=\!\rho_{\mathrm{c}}(x_{\tau},t\!-\!\tau|x,t)\rho_{1}(x,t)$ [with $\int_{\Omega}\rho_{2}(x,t;x_{\tau},t-\!\tau)\mathrm{d}x_{\tau}\!=\!\rho_{\mathrm{1}}(x,t)$]. Equation (4) is hence not self-sufficient. The delayed FPE for $\rho_{2}$, on the other hand, involves $\rho_{3}$ (as explicitly shown in the Appendix A), and so forth. Thus, an infinite hierarchy of equations emerges, whose $n\!+\!1$st member depends on the $n$-time probability density $\rho_{n}(x,t;...;x_{n\tau},t-n\tau)$. By finding suitable approximations for $\rho_{\mathrm{c}}$, Eqs. (4) can be closed and the hierarchy truncated. This is our objective in the present work.

II.2 Earlier approaches to approximate the one-time probability density

One such approach, involving a first-order perturbation-theoretical (PT) ansatz for the density on level of the FPE, was previously introduced in Ref. Frank (2005a). Within PT, the delay force $F_{\mathrm{d}}$ is regarded as a small perturbation to the non-delayed dynamics, i.e., $|F_{\mathrm{d}}|\!\ll\!|F_{\mathrm{s}}|$. The resulting first-order equation has the same form as Eq. (4), but the $\widehat{F}$-term [Eq. (4b)] contains the conditional probability density with respect to the unperturbed, i.e., non-delayed, system, which we will refer to as $\rho^{F_{d}\equiv 0}_{c}$. The latter follows a closed (Markovian) FPE, and consequentially, the combination of the two FPEs is self-sufficient. From a practical perspective, the PT approach is particularly appealing for systems where an analytical expression for $\rho^{F_{d}\equiv 0}_{\mathrm{c}}$ 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 approach often requires additional approximations. In the present study we use two different approaches, which we will specify in Sec. III.1.

Yet another approach involves a small delay expansion on the level of the LE Guillouzic et al. (1999, 2000); Longtin (2010). More specifically, the authors of Ref. Guillouzic et al. (1999, 2000) suggest a Taylor expansion up to linear order of the total deterministic force (and of the noise intensity when delayed noise is considered) in powers of $\tau$. The Taylor expansion is somewhat problematic, since it involves the derivative of the position $\chi$, which is a stochastic variable and hence not continuously differentiable. Performing this Taylor expansion ad hoc hence implies a certain inaccuracy, which may become especially apparent in the range of large noise intensities (large $D_{0}$), see Frank (2005a) for a discussion.

The novel approach introduced in the present paper works on the basis of the FPE, like the perturbation theory. As we will demonstrate, our ansatz provides an improved approximation of the one-time probability density for the considered nonlinear static forces.

II.3 Force linearization closure

In the following, we focus on non-equilibrium steady states (NESS). Indeed, many Langevin systems described by Eq. (2) automatically reach a NESS, defined by ${\partial_{t}}\rho_{\mathrm{1,ss}}(x,t)\!=\!0$. Steady states generally allow for more analytical treatment than transients. Throughout this work we denote steady state quantities with the subscript “ss”. Since the steady state conditional probability only depends on the time difference, and not on the instances of time themselves, we will further use the shortened notation: $\rho_{\mathrm{c,ss}}(x_{\tau}|x;\tau)=\rho_{\mathrm{c,ss}}(x_{\tau},t\!-\!\tau|x,t)$.

The basis of our approach is that, within a NESS, the FPE hierarchy can be solved exactly when all deterministic forces are linear in $x$ and $x_{\tau}$ Küchler and Mensch (1992); Frank and Beek (2001). A further, yet less restrictive requirement is that the system obeys natural boundary conditions, i.e., $\lim_{x\to\pm\infty}\rho(x,t)\!=\!0$ and $\Omega\!=\!\mathbb{R}$. Given this background, the main idea of our approach is to achieve a closed approximate FPE for $\rho_{1}$ by linearizing the deterministic forces in all members of the infinite FPE hierarchy apart from the first one. Due to the involved linearization procedure, which we outline in detail below, we call our approach “force linearization closure” (FLC).

II.3.1 Linearization of the deterministic forces

We start by considering the (time-dependent) energy landscape resulting from the total force (1). As a simple estimate for the total steady state energy landscape, we assume at this step the system to be at rest, i.e., $\chi(t-\tau)\!\equiv\!\chi(t)$, such that the total (static) potential is given by

[TABLE]

In general, $V_{\mathrm{STAT}}$ may have multiple (local) minima (defined by $F\!=\!0$ and $F^{\prime}\!>\!0$). We number each minimum $x^{i}_{0}$ with an integer $i\in\!\mathcal{I}$. The index set $\mathcal{I}$ hence contains one element for each local minimum. To formulate the FLC, we split the spatial domain $\Omega$ into non-overlapping intervals $\Omega^{i}$, such that both boundaries of $\Omega$ and each local maximum of $V_{\mathrm{STAT}}$ represents a bound of an interval, and every $\Omega^{i}$ contains exactly one minimum $x^{i}_{0}$.

By this procedure, we obtain a splitting of $\Omega$ into subdomains $\Omega^{i}$, whose union is again the entire spatial domain: $\Omega\!=\!\bigcap_{{i\in\mathcal{I}}}\Omega^{i}$. For static potentials with a single minimum $x^{j}_{0}$ , this means $\Omega^{j}\!=\!\Omega$. For each $i$, we then perform a Taylor expansion of the entire deterministic force $F\!=\!F_{\mathrm{s}}+F_{\mathrm{d}}$ in both spatial variables $x$ and $x_{\tau}$. More specifically, we expand around the deviations with the enclosed minimum, i.e., $\Delta x^{i}\!=\!x-x^{i}_{0}$ and $\Delta x_{\tau}^{i}\!=\!x_{\tau}-x^{i}_{0}$. Neglecting all terms of quadratic orders $\mathcal{O}(\Delta{x^{i2}})$, $\mathcal{O}(\Delta x^{i}\Delta x^{i}_{\tau})$, $\mathcal{O}(\Delta{x_{\tau}^{i2}})$, or higher, we obtain

[TABLE]

where $-\alpha^{i}$ and $-\beta^{i}$ are the first-order derivatives of $F$ with respect to $x$ and $x_{\tau}$, respectively, evaluated at the minimum $x^{i}_{0}$. [Note that since we expand around the minima of $V_{\mathrm{STAT}}$ (5), we always have $F(x^{i}_{0},x^{i}_{0})\!=\!0$, such that the constant terms vanish.]

This linearization procedure yields an approximation of the steady state energy landscape composed of a sequence of quadratic polynomials with time-dependent centers, each subdomain $\Omega^{i}$ reaching from one local maximum to the following one (or to a bound of the entire spatial domain).

II.3.2 Analytical solution for linearized forces

Now, we will turn back to the full system which involves thermal noise. Without further reasoning, one would expect that the linearization of the deterministic forces renders a good approximation of the stochastic dynamics, whenever the probability density close to the minima of the (approximate) total potential in the steady state [Eq. (5)] is large. This is, for instance, the case, when the potential barriers are high compared to thermal fluctuations.

The following steps are performed separately for each subdomain $\Omega^{i}$. We first apply the linearization $F\approx F^{\mathrm{lin},i}$ [see Eq. (6)] to the force terms in all members of the infinite FPE hierarchy starting from the second, i.e., the equation for the two-time steady state density $\rho_{2,\mathrm{ss}}(x,t;x_{\tau},t-\tau)$ (see Appendix A for the general from of this FPE). Assuming natural boundary conditions at every subdomain bound, the FPEs for all $\rho_{\mathrm{n,ss}}$ can be solved by multivariate Gaussian distributions Küchler and Mensch (1992). The second member of the FPE hierarchy in the linear case and the corresponding analytical solutions $\rho^{\mathrm{lin}}_{\mathrm{2,ss}}$ and $\rho^{\mathrm{lin}}_{\mathrm{3,ss}}$ are given in the Appendices B and C. For each subdomain, we thus have access to the function $\rho_{\mathrm{2,ss}}\equiv\rho^{\mathrm{lin},i}_{\mathrm{2,ss}}$ (the explicit formula for $\rho^{\mathrm{lin},i}_{\mathrm{2,ss}}$ can also be found in Frank et al. (2003)). The corresponding conditional probability density reads

[TABLE]

with the coefficients

[TABLE]

We recall $\Delta x^{i}\!=\!x-x^{i}_{0}$, where $x^{i}_{0}$ is the enclosed minimum of $V_{\mathrm{STAT}}$ [Eq. (5)]. Note that $\omega^{i}$ becomes imaginary if $|\alpha^{i}|\!<\!\beta^{i}$, such that the hyperbolic functions in $K(\alpha^{i},\beta^{i},\tau)$ [see Eq. (8b)] convert to trigonometric ones. In this case, there exist critical $\tau$ values $\tau_{c}\omega=\arccos(-\alpha^{i}/\beta^{i})+2\pi\kappa,~{}\forall\kappa\in\mathbb{Z}$, for which $K^{i}$ diverges, and no stable NESS exists. There is also no stable state, when $-\beta^{i}\leq\alpha^{i}$, or $\tau\alpha^{i}\!\leq\!-1$. Only if $-\alpha^{i}\!<\!\beta^{i}\!<\!\alpha^{i}$, a steady state is approached for all $\tau$. A derivation and discussion of the steady state conditions is presented in Küchler and Mensch (1992).

With Eq. (II.3.2), the first member of the approximate FPE hierarchy (where the force linearization is applied to the second and all higher members) accordingly reads

[TABLE]

In Eq. (9a), $\rho^{\mathrm{FLC}}_{\mathrm{1,ss}}$ is the steady-state one-time probability density obtained within the FLC approach. Notice that the deterministic force $F$ appearing in the integral in Eqs. (9b) is not linearized.

Equations (II.3.2-9b) [together with the linearization rule (6)] form a closed set of equations. The here presented closure of the FPE hierarchy is the cornerstone of our approach. Since we have a general expression for $\rho_{\mathrm{c,ss}}^{\mathrm{lin},i}$ [Eq. (II.3.2)] with which $\widehat{F}^{\mathrm{FLC}}$ can be calculated separately, the FLC formally converts the delayed FPE for $\rho_{\mathrm{1,\mathrm{ss}}}$ itself into a closed, quasi-Markovian one [Eq. (9a)]. The fact that the delay is still present in our approximate FPE, becomes apparent by considering the special case of a linear force $F$. Then, Eqs. (9) coincide with the corresponding exact delayed FPE. Furthermore, one can easily see from Eqs. (9) that the usual FPE is recovered when the delay time vanishes, or when $F$ becomes independent of $x_{\tau}$: in both cases $\widehat{F}^{\mathrm{FLC}}\!=\!F$.

II.3.3 Vanishing steady state probability current

The key equations of the FLC [Eqs. (II.3.2-9b)] are valid for steady states of, in principle, arbitrary systems which can be meaningfully linearized according to the procedure described in Sec. II.3.1. A particularly simple situation arises, if additionally the steady-state probability current given by $J\!=\![\widehat{F}^{\mathrm{FLC}}(x)-\gamma D_{0}{\partial_{x}}]\rho^{\mathrm{FLC}}_{\mathrm{1,ss}}(x)$ vanishes. If $J=0$, the formal solution of Eq. (9a) takes the simple Boltzmann-like form

[TABLE]

with the effective static potential

[TABLE]

where $\widehat{F}^{\mathrm{FLC}}$ is given in Eq. (9b), and $\hat{x}\in\Omega$ is arbitrary but fixed. Here and in the following, we denote the normalization constant by $Z$.

We stress that the assumption $J\!=\!0$ simplifies the analysis, but it is not a necessary condition. Also for non-zero currents, Eq. (9a) can be treated using standard techniques for (Markovian) Fokker-Planck equations Risken (1984); Gardiner (2002).

III Applications

We now apply the FLC to two generic examples, involving a multistable and a bistable static potential combined with a linear delay force. In the subsequent section III.1, we first define this force and provide some results which apply to both model systems. In Secs. III.2 and III.3 we then present results from the FLC and compare its performance to the two main approaches known from the literature, i.e., the perturbation theory (PT) and the small delay expansion. As a test of the different approximations, we provide numerical results from Brownian dynamics (BD) simulations. Details about the numerical methods are given in the Appendix D.

III.1 Linear delay force, general results

We consider a linear delay force, with amplitude $k\geq 0$,

[TABLE]

which vanishes for $k\!\rightarrow\!0$ or $\tau\!\rightarrow\!0$. The particular force in Eq. (12) can be associated with the delayed confining potential $V_{\mathrm{d}}(x,x_{\tau})=(k/2)[x-x_{\tau}]^{2}$. Such quadratic feedback potentials are commonly used to model optical traps Kotar et al. (2010); Kim and Qian (2007); Wördemann (2012), which are implemented in many experimental setups to control colloidal particles Kotar et al. (2010); Kim and Qian (2007); Qian et al. (2013); *Balijepalli2012; Balijepalli et al. (2012). From a more general perspective, the delay force in Eq. (12) is of Pyragas type Pyragas (1992), and has been extensively studied in the context of chaos control Schöll and Schuster (2008). It is important to note that the FLC generally also applies to systems with nonlinear delay forces.

For both exemplary systems, we use natural boundary conditions. Therewith, the steady state probability current vanishes irrespective of the particular form of the static potential $V_{\mathrm{s}}$, and solution (10) readily applies. We can give a general expression for $\rho_{\mathrm{1,ss}}^{\mathrm{FLC}}$ as follows. The delay force is already linear, and yields $\beta^{i}\!=\!-k,\forall i\in\mathcal{I}$, where $\beta^{i}$ and $k$ are the coefficients appearing in Eqs. (6) and (12), respectively. After performing several Gaussian integrals, Eq. (9b) yields

[TABLE]

$\widehat{F}^{\mathrm{FLC}}$ can be further simplified by using the identity

[TABLE]

where the quantities $d^{i}_{1}$ and ${d^{i}_{2}}$ are given in Eqs. (8c) and (8d). This yields the piecewise defined effective potential

[TABLE]

for $x\in\Omega^{i}$. Inserting Eq. (15) into Eq. (10) one obtains $\rho^{\mathrm{FLC}}_{\mathrm{1,ss}}$. Furthermore, since $F_{\mathrm{d}}(x,x)\!=\!0\,\forall x\in\Omega$, our estimate for the steady state energy landscape coincides with the static potential $V_{\mathrm{STAT}}\!=\!V_{\mathrm{s}}$, and hence $x^{i}_{0}$ and the bounds of $\Omega^{i}$ are readily determined by the extrema of $V_{\mathrm{s}}$. Thus, for a given $V_{\mathrm{s}}$, one only needs to calculate $\alpha^{i}$ by linearizing $F_{\mathrm{s}}$, and therewith the coefficient $d^{i}_{2}$ [Eq. (8d)].

Before specifying $V_{\mathrm{s}}$, we review some results from the two main earlier approaches, for the case of the linear delay force [Eq. (12)] (and natural boundary conditions). Just like the FLC [Eq. (10)], the small delay expansion gives rise to a Boltzmann-distributed one-time probability density

[TABLE]

Here, the effective potential $V_{\mathrm{eff}}^{\mathrm{approx}}\equiv V^{\mathrm{small}\tau}_{\mathrm{eff}}$ reads Frank (2005a)

[TABLE]

The application of the perturbation-theoretical approach, on the other hand, is not straightforward. The closed approximate FPE from the PT involves the conditional probability of the corresponding unperturbed ($F_{\mathrm{d}}\!\equiv\!0$) system, $\rho^{F_{d}\equiv 0}_{\mathrm{c,ss}}(x_{\tau}|x;\tau)$. Since for most static potentials (including the ones we will consider in Secs. III.2 and III.3) no analytical expression for this quantity is available, further approximations become inevitable. As suggested in Ref. Frank (2005a), one can for this purpose utilize the “short time propagator” Risken (1984)

[TABLE]

which is derived and discussed in Ref. Risken (1984) on p. 73 and prior. The usage of the short time propagator implies a first-order approximation in $\tau$ Risken (1984); Frank (2005a). The resulting one-time probability density is also the Boltzmann distribution Eq. (16) with $V_{\mathrm{eff}}^{\mathrm{approx}}\equiv V^{\mathrm{PT}}_{\mathrm{eff}}$ given by

[TABLE]

Interestingly, when the short time propagator is used for $\rho^{F_{d}\equiv 0}_{\mathrm{c,ss}}$, the PT approach obviously yields qualitatively very similar results as the small $\tau$ expansion. In particular, both approaches render effective potentials that are proportional to the static potential [Eqs. (17) and (19)]. This means, also the approximate densities have the same functional form as in the case $F_{\mathrm{d}}\!\equiv\!0$. In fact, the only remaining effect of the delay force, according to these approximations, is a constant ($x$-independent) factor within the exponent of $\rho^{\mathrm{approx}}_{\mathrm{1,ss}}$. Moreover, the small $\tau$ expansion and the PT with short time propagator even become equivalent for small $\tau k$, since then these factors coincide $1/(1-k\tau)=(1+k\tau)+\mathcal{O}([\tau k]^{2})$.

Alternatively, we propose to approximate $\rho^{F_{d}\equiv 0}_{\mathrm{c,ss}}$ by the conditional probability from the corresponding unperturbed ($F_{d}\equiv 0$), linearized system, i.e., the corresponding Ornstein–Uhlenbeck process, given in Ref. Risken (1984) (p. 100). The resulting effective potential from the PT approach thus reads

[TABLE]

As opposed to the result obtained with the short-time propagator [Eq. (19)], the effective potential in Eq. (III.1) has a different functional form than $V_{\mathrm{s}}(x)$. More specifically, the delay force now effectively adds to the static potential a fixed quadratic potential around the minima of $V_{\mathrm{s}}$. Since the linear delay force [Eq. (12)] is indeed expected to trap the particle in a quadratic confining potential (with history-dependent center), this effective potential appears to be more realistic than Eq. (19). However, the history-dependency of the position of the “trap” imposed by $F_{\mathrm{d}}$ is, of course, also not captured by this approximations.

As a first consistency check, one readily sees from Eqs. (15-17),(19),(III.1) that all four approximations give rise to Boltzmann-distributed density profiles, which, in the limits of vanishing delay force ($\tau\rightarrow 0$ or $k\rightarrow 0$), all coincide with the exact result for the non-delayed (Markovian) FPE, i.e., $\rho_{\mathrm{1,ss}}\propto\exp[-{V_{\mathrm{s}}(x)/k_{\mathrm{B}}\mathcal{T}}]$. However, they yield different effective potentials for finite delay forces, which we will compare in the following.

Finally, we note that considering the force (12) rather than a force with two different amplitudes $F_{\mathrm{d}}(x,x_{\tau})\!=\!-k^{\prime}x-kx_{\tau},$ does not imply a loss of generality concerning the results in this section. A minor adjustment required if $k^{\prime}\!\neq\!k$, is that $x^{i}_{0}$ are then the minima of $V_{\mathrm{STAT}}\neq V_{\mathrm{s}}$.

III.2 Periodic static potential

We start by considering the periodic “washboard” potential

[TABLE]

with barriers of height $\Delta\!V_{\mathrm{}}$ and minima at $x^{i}_{0}\!=\!2\pi\sigma i$ and $i\!\in\!\mathcal{I}\!=\!\mathbb{Z}$. Particles in sinusoidal potentials represent a well-studied paradigm to model transport on rough surfaces Juniper et al. (2016). Linearizing the static potential yields $\alpha^{i}\!=\!\Delta\!V_{\mathrm{}}/(2\sigma^{2})+k$ (and $\beta^{i}\!=\!-k$) for all $i\in\mathcal{I}$. On all subdomains $\Omega^{i}\!=\![(2i-1)\pi\sigma,(2i+1)\pi\sigma]$, the density is hence given by the Boltzmann distribution (10) with effective potential (15), where $V_{\mathrm{s}}$ is given in Eq. (21) and

[TABLE]

with $\gamma\omega\sigma=\sqrt{({\Delta\!V_{\mathrm{}}}/{2\sigma})^{2}+k\Delta\!V_{\mathrm{}}}$ and $\Delta x^{i}=x-x^{i}_{0}$.

Figure 1 shows the density profile in one potential well generated with BD simulations of the delayed LE (2), on the one hand, and the corresponding FLC approximation (22), on the other hand, for an exemplary parameter choice. Here and throughout the entire paper, we set the barrier heights to $\Delta\!V_{\mathrm{}}\!=\!{8}\,k_{\mathrm{B}}{\cal T}$. One clearly sees that the FLC generates a very good approximation of the probability density distribution. Similar convincing results are found for all other considered values of $\tau$ and $k$ in the tested range $\tau/\tau_{\mathrm{B}}\in[10^{-2},10]$ and $k/(\sigma^{-2}k_{\mathrm{B}}\mathcal{T})\in[0,32]$. To compare approximation and simulation results more systematically, we calculate the moments $\mu_{n}$ of the distributions within one potential well (for $x$ values between two maxima of $V_{\mathrm{s}}$), e.g., of $\rho_{\mathrm{1,ss}}(-\pi\sigma\leq x\leq\pi\sigma)$. Due to the spatial inversion symmetry of the total potential around the enclosed minimum, all odd moments vanish. We hence consider the second central moment $\mu_{2}\!=\!\langle(\chi-\mu_{1})^{2}\rangle_{\mathrm{ss}}$ (with $\mu_{1}\!=\!\langle\chi\rangle_{\mathrm{ss}}$) of one peak of the distribution. Figure 2 shows $\mu_{2}$ vs. the delay time $\tau$ for two exemplary values of delay force amplitude. We find that in the presence of the delay force (12), the variance of the distribution is reduced. This is not surprising, since the delay force arises from a quadratic potential which “traps” the particle and hence enhances the confining effect of $V_{\mathrm{s}}$. Within the considered parameter regime, the variance is seen to decrease with increasing $\tau$, until a saturation value is approached at $\tau^{\mathrm{W}}_{\mathrm{sat}}\!\approx\!0.5\,\tau_{\mathrm{B}}$ (due to the strongly increasing simulation times, we did not simulate much larger $\tau$).

On the level of the probability densities and the FPE, the saturation of $\mu_{2}$ at large $\tau$ indicates that in the considered range of delay times, the delay-averaged drift $\widehat{F}$, and hence $\rho_{\mathrm{c,ss}}(x_{\tau}|x;\tau)$, are essentially constant for $\tau\!>\!\tau^{\mathrm{W}}_{\mathrm{sat}}$. In other words, the conditional probability for the time difference $\tau$ remains essentially unchanged, when $\tau$ is further increased. This suggests some kind of relaxation mechanism within the valley, where the history of the stochastic process becomes “washed out” on the level of ensemble averaged quantities, like the probability densities. In this context it is interesting to note that, for the case without delay force ($k\!=\!0$), the relaxation time within a potential well is of the order $\tau^{\mathrm{W}}_{\mathrm{ir}}\!\approx\!{\gamma}/{V_{\mathrm{s}}^{\prime\prime}(x_{\mathrm{min}})}=2\,\tau_{\mathrm{B}}k_{\mathrm{B}}\mathcal{T}/\Delta\!V_{\mathrm{}}$ (see Ref. Gardiner (2002), p. 348). In the present case, $\Delta\!V_{\mathrm{}}\!=\!{8}\,k_{\mathrm{B}}{\cal T}$, that is $\tau^{\mathrm{W}}_{\mathrm{ir}}\!\approx\!0.25\,\tau_{\mathrm{B}}$. The saturation of $\mu_{2}$ thus sets in after the density relaxation within a potential well, consistent with our expectation.

A different situation occurs for much larger $\tau$ values than the ones considered in Fig. 2, in particular, when $\tau$ gets into the range of the mean escape times (numerical results for the escape times are provided in Sec. III.4). Then we expect again a $\tau$-dependency of $\rho_{\mathrm{c,ss}}$ and therewith of $\mu_{2}$, since the transport between the potential valleys becomes important.

We now compare the density from the FLC approach with 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–Uhlenbeck approximation [see Eq. (III.1) and above]. As visible in Fig. 1 and 2, the PT generally overestimates the height of the density peak. Approximating the conditional probability with the short time propagator yields worse results than using $\rho_{\mathrm{c,ss}}$ of the corresponding Ornstein–Uhlenbeck process. This observation matches our expectations, see Sec. III.1. At least when the Ornstein–Uhlenbeck approximation is used, the PT reproduces the quantitative behavior of the function $\mu_{2}(\tau)$, see Fig. 2. On the contrary, the small delay expansion fails completely for larger delay times.

We conclude that the FLC generates the best approximation of the steady-state density distribution in the periodic potential, especially in the regime of large $\tau$ or large $k$. This is the regime where the “perturbation”, i.e. $|F_{\mathrm{d}}|$, is not small compared to the force applied by the static potential, such that the basic assumption of the perturbation theory is not fulfilled. The small delay expansion, on the other hand, involves a truncated Taylor expansion in $\tau$, making it plausible that also this approach fails for large delay times.

III.3 Bistable static potential

As a second example, we consider the double-well potential

[TABLE]

with a barrier of height $\Delta\!V_{\mathrm{}}$ (set to $\Delta\!V\!=\!8\,k_{\mathrm{B}}{\cal T}$) at $x/\sigma\!=\!0$ and two minima at $x^{i}_{0}=i\sigma$ with $i\in\mathcal{I}\!=\!\{-1,1\}$.

A Brownian particle in a double-well potential is a generic model for a bistable noisy system Risken (1984); Gardiner (2002), which in recent years has been considered also under the impact of a linear delay force Tsimring and Pikovsky (2001); Masoller (2003); Du and Mei (2015); Xiao (2016); Piwonski et al. (2005). As opposed to the corresponding expression given in Eq. (12) (which corresponds to a moving “optical tweezer”), these earlier theoretical studies consider a slightly different delay force involving solely the delayed particle position, i.e., $F_{\mathrm{d}}(x,x_{\tau})=-kx_{\tau}$. Our approach can, however, as well be applied to such a delay force.

The force linearization (according to the procedure described in Sec. II.3.1) yields $\alpha^{i}\!=\!8\,\Delta\!V_{\mathrm{}}/\sigma^{2}+k$ for both $i$ (and $\beta^{i}\!=\!-k$). For $x/\sigma\in\![-\infty,0]$ and $x^{i=-1}_{0}\!=\!-\sigma$, or $x/\sigma\in\![0,\infty]$ and $x^{i=1}_{0}\!=\!\sigma$, the effective potential $V^{\mathrm{FLC}}_{\mathrm{eff}}$ is given by Eq. (15) with

[TABLE]

with $\gamma\omega\sigma=\sqrt{({8\Delta\!V_{\mathrm{}}}/{\sigma})^{2}+k\Delta\!V_{\mathrm{}}}$.

Also for this model, the FLC clearly renders a very good 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 one peak of the bimodal distribution, i.e., $\rho_{\mathrm{1,ss}}(x\!>\!0)$, as a function of the delay time. More specifically, we use the skew $\mu^{x>0}_{3}\!=\!\int_{x>0}[(x-\mu_{1})/\sqrt{\mu_{2}}]^{3}\rho_{1,\mathrm{ss}}(x)\mathrm{d}x$, i.e., the third central moment normalized with the standard deviation. The comparison with the numerical data in Fig. 4 reveals that the FLC approach renders reasonably good predictions of $\mu^{x>0}_{3}$. We moreover find convincing results for all other tested parameters within the considered ranges $\tau/\tau_{\mathrm{B}}\in[10^{-2},5]$ and $k/(\sigma^{-2}k_{\mathrm{B}}\mathcal{T})\in[0,18]$. Please note that this range spans from rather small delay forces that merely affect the shape of the energy landscape, to $k$ values so high that there is no second minimum of the total potential, when $\chi(t-\tau)$ rests in the first minimum. For all these qualitatively different cases, the agreement between FLC and numerical results is good. However, it is not as accurate as that for the second moment in the case of the periodic potential (see Fig. 2). One reason for that might be the fact that the approximate potential is, by construction, symmetric with respect to the enclosed minimum (within each subdomain). In this sense, the double-well potential, which is asymmetric around each minimum, is not as well approximated as the symmetric “washboard” potential. This manifests in a larger magnitude of the higher order Taylor terms that are neglected within the linearization procedure (6). We also note that the third moment is per se very sensitive against inaccuracies, since it involves cubic terms in the relative position. The deviations between the FLC and the exact result are hence expected to be larger.

Regarding the impact of increasing $k$ and $\tau$, we see from Fig. 4 that the skew $\mu^{x>0}_{3}$ becomes smaller, which means that the probability distribution becomes more symmetric. This is not surprising since the pure delay potential is quadratic (i.e., symmetric) in the system state variable. Moreover, similar to the second moment in the case of the “washboard” potential (see Fig. 2), $\mu^{x>0}_{3}$ approaches a saturation value (within the considered parameter range). The corresponding delay time is about one order of magnitude smaller than in the “washboard” case: $\tau^{\mathrm{D}}_{\mathrm{sat}}\!\approx\!0.1\tau_{\mathrm{B}}$ for $\Delta\!V_{\mathrm{}}\!=\!8\,k_{\mathrm{B}}\mathcal{T}$. This value of $\tau^{\mathrm{D}}_{\mathrm{sat}}$ is, as in the case of the “washboard”, significantly larger than that for the relaxation time $\tau^{\mathrm{D}}_{\mathrm{ir}}$ within a well. For the bistable system at $k\!=\!0$ one finds Gardiner (2002) (see p. 348) $\tau^{\mathrm{D}}_{\mathrm{ir}}/\tau_{\mathrm{B}}\!\approx\!1/8\,(k_{\mathrm{B}}\mathcal{T}/\Delta\!V_{\mathrm{}})=1/64$ at $\Delta\!V_{\mathrm{}}=8\,k_{\mathrm{B}}\mathcal{T}$. Thus, the saturation of $\mu_{3}$ can be explained by the same arguments as in the case of the “washboard” potential.

Figure 3 and 4 also show the results according to the PT [Eqs. (19) and (III.1)], and from small delay expansion, see Eq. (17). Very similar to the case of the “washboard”, only the PT with Ornstein–Uhlenbeck approximation is capable of reproducing the quantitative behavior (especially for large $\tau$), and the FLC clearly provides the best approximation of the one-time probability.

III.4 Escape times

III.4.1 The Kramers-FLC estimate

The escape time $\tau_{\mathrm{K}}$ denotes the average time $\chi(t)$ spends in the vicinity of a potential minimum, until it leaves the valley by jumping to an adjacent one. By construction, mean escape times involve two locations and two instances in time, and hence are connected to two-time probabilities. Within our FLC approach, the higher order probabilities $\rho^{\mathrm{FLC}}_{\mathrm{(n>2),ss}}$ are multivariate Gaussian densities. These belong to a quadratic static potential with only one valley and no escape process. Therefore, there is no direct route to calculate the mean escape time on the basis of the FLC approach.

However, in the limit of small noise intensities, we can find an approximation for $\tau_{\mathrm{K}}$ by using Kramers theory for (Markovian) systems characterized by a static potential landscape $U$. When the potential barriers $\Delta U$ are large compared to the thermal energy, i.e., when $\gamma D_{0}\!\ll\!\Delta U$, the Kramers theory provides an Arrhenius formula Gardiner (2002); Gernert et al. (2014); Hänggi et al. (1990) for the escape rate, $r_{\mathrm{K}}\!=\!{\sqrt{U^{\prime\prime}(x_{\mathrm{min}})|U^{\prime\prime}(x_{\mathrm{max}})|}}/({2\pi\gamma})\exp\left(-{\Delta U}/{\gamma D_{0}}\right)$, (with $U^{\prime\prime}(x_{\mathrm{ex}})\!=\!\partial_{xx}U(x)|_{x=x_{\mathrm{ex}}}$, with $x_{\mathrm{ex}}\in\{x_{\mathrm{min}},x_{\mathrm{max}}\}$). This estimate for $r_{\mathrm{K}}$ relies on an equilibrium approximation for $\rho_{1}$ and does not involve $\rho_{2}$. Due to the two symmetric ways to leave a valley of the “washboard” potential, the resulting mean escape time can be approximated by $\tau_{\mathrm{K}}\!=\!{1}/({2r_{\mathrm{K}}})$ Hänggi et al. (1990) and for the double-well potential with the unique direction to exit each valley, by $\tau_{\mathrm{K}}={1}/{r_{\mathrm{K}}}$.

To use the Kramers theory in the context of the FLC approach, we note that the first member of the FPE hierarchy [Eq. (9a)] is formally identical to a Markovian (non-delayed) FPE for $\rho_{\mathrm{1,ss}}$ with a static effective potential $V^{\mathrm{FLC}}_{\mathrm{eff}}$ given by Eq. (11). Thus, we can directly apply the Arrhenius formula with $U\equiv V^{\mathrm{FLC}}_{\mathrm{eff}}$. For the linear delay force in our examples, one just needs to substitute $\Delta U\!=\!\Delta\!V_{\mathrm{}}+(k/2)(x_{\mathrm{max}}-x_{\mathrm{min}})$ and $U^{\prime\prime}(x_{\mathrm{ex}})=V^{\prime\prime}_{\mathrm{s}}(x_{\mathrm{ex}})+k(1-|d^{i}_{2}|)$. Taken altogether, our estimate of $\tau_{\mathrm{K}}$ involves two separate approximations: the FLC to obtain $V^{\mathrm{FLC}}_{\mathrm{eff}}$ and the application of the Arrhenius formula to this non-Markovian system (and implicitly, all simplifications made within the Kramers theory).

III.4.2 Escape times in the delayed “washboard” potential

We now apply the Kramers-FLC approximation to calculate the escape times for the delayed “washboard” potential. Exemplary results for $\tau_{\mathrm{K}}$ are plotted in Fig. 5 together with numerical data. For $F_{\mathrm{d}}\!\equiv\!0$, the BD simulation results and the Kramers-FLC estimate coincide. In the regarded parameter range, $\tau_{\mathrm{K}}$ generally increases with $k$ and $\tau$, and the FLC results are roughly in agreement with the numerical results. This is rather remarkable, due to the crudeness of applying the Kramers theory to the non-Markovian system. At large values of $\tau$, one observes a qualitatively different behavior: while the escape times resulting from the Kramers-FLC estimate saturate, the corresponding BD data continue to increase with $\tau$. This difference becomes particularly prominent for large $k$. We propose that the discrepancy can be explained as follows. As already discussed in Sec. III.2, the conditional probability saturates at a finite value $\tau_{\mathrm{sat}}$, which is related to the intra-well relaxation time $\tau_{\mathrm{ir}}$. Accordingly, also $V^{\mathrm{FLC}}_{\mathrm{eff}}$, and therewith $\tau_{\mathrm{K}}$ from the Kramers-FLC approximation, saturate at $\tau_{\mathrm{sat}}$. In the true non-Markovian system, the escape times not only depend on $\rho_{\mathrm{2,ss}}$, but also on higher $n$-time probabilities. Therefore they don’t necessarily saturate together with $\rho_{\mathrm{2,ss}}$. Indeed, our numerical investigations reveal that the saturation of $\tau_{\mathrm{K}}$ sets in at much higher $\tau$ values, see Fig. 5. For larger $\tau$, another, not yet discussed time scale becomes increasingly important, that is the jump duration time. In the present system, this time is roughly $2\,\tau_{\mathrm{B}}$ (according to BD). For $\tau$ in the range of the jump duration times, the temporal changes of the total potential and the significant changes of the system state variable $\chi$ occur on similar time scales. The interplay between both motions then leads to new (inter-well) dynamical behavior not captured by our approach. It hence becomes less justified to treat the delayed system as a quasi-Markovian one, whose stochastic (non-delayed) dynamics evolve in a static (effective) potential.

One example for such new dynamical behavior are quasi-regular oscillations of $\chi$ between two valleys. In fact, we have observed analogous dynamics in the delayed bistable system. For both static potentials, the delay-induced oscillations have a mean period of about $\tau$ and start at random times. We further observe that they occasionally pause and randomly set in again. Figure 6 shows BD results for the distribution of waiting times between sequential jumps in the delayed double-well potential for an exemplary parameter setting. One sees that most waiting times lie within an interval $\Delta t_{\mathrm{jump}}\in[0,\tau]$ around a single, yet broad peak at about $\tau/2$. This illustrates the stochastic character of the oscillations with mean period of about $\tau$. Similar quasi-regular oscillations can been found in the double-well potential with linear delay force $F_{\mathrm{d}}(x,x_{\tau})\propto x_{\tau}$ as reported in Tsimring and Pikovsky (2001) and further discussed in Masoller (2003). Interestingly, this version of delayed bistable system has moreover been shown to exhibit coherence resonance, i.e., that the regularity of the oscillations is maximal at a certain finite noise intensity. Preliminary numerical studies suggest that our model systems also display coherence resonance. Further investigations on the delay-induced oscillatory behavior in our systems are in progress.

IV Spatial autocorrelation function

In this final section, we aim to discuss a fundamental difference between the FLC and the PT, which has not been in the focus so far. This point also provides an explanation why the FLC performs better on the exemplary systems considered in the present work.

We recall that both approaches yield closed FPEs for $\rho_{1}$ via an approximation of the conditional probability $\rho_{\mathrm{c,ss}}(x_{\tau},t\!-\!\tau|x,t)$, and therewith of the two-time probability density $\rho_{\mathrm{2,ss}}(x,t;x_{\tau},t\!-\!\tau)$ of the delayed system. The essential difference is that in the framework of the FLC, $\rho_{\mathrm{c,ss}}$ is approximated by the corresponding function of the linearized delayed system, while in the PT approach, $\rho_{\mathrm{c,ss}}$ stems from the corresponding non-delayed system. The FLC hence involves a $\rho_{2}$ related to a non-Markovian system, while the PT utilizes a “Markovian $\rho_{2}$”. On the level of the two-time probability density, the difference between a Markovian and a non-Markovian system manifests itself, e.g., via the spatial autocorrelation $C(\tau)=\langle\chi(t)\chi(t-\tau)\rangle_{\mathrm{ss}}-\langle\chi(t)\rangle^{2}_{\mathrm{ss}}$ obtained by integrating $\rho_{2}$ times $xx_{\tau}$ over both spatial variables and subtracting the squared first moment, see Eq. (33). While in the Markovian case $C(\tau)$ always vanishes for large $\tau$, reflecting the loss of memory, it stays finite for typical non-Markovian systems, revealing that the system states at the times $t$ and $t-\tau$ are correlated for large $\tau$. This is illustrated in Fig. 7, where we have plotted $C(\tau)$ for two of the considered non-Markovian systems (the “washboard” potential and the quadratic potential received by the force linearization, both with linear delay force), and for the corresponding Markovian systems (obtained by setting the delay force zero). The fact that only the FLC involves a non-Markovian quantity to approximate $\rho_{\mathrm{c,ss}}$, which accounts for these finite correlations, gives an explanation why it yields better results for large delay times than the perturbation theory.

V Conclusion

In this work, we introduce the FLC as a novel approach to close the FPE for the steady state probability density of delayed systems. The new closure is achieved by linearizing the total deterministic force in all members of the FPE hierarchy starting from the second. We have specified on a linear delay force, which can experimentally be implemented, for instance, by an optical laser tweezer acting on a polarizable colloidal particle Kotar et al. (2010). However, the novel approach generally also applies to nonlinear delay forces.

By considering two exemplary model systems, we present numerical evidence that the novel approach generates indeed a very good approximation of the one-time probability density. In fact, the FLC performs much better here compared to the small delay (Taylor) expansion introduced in Ref. Guillouzic et al. (1999), and to the perturbation theory on the level of the FPE from Ref. Frank (2005a). We also note that the FLC is somewhat more general than the PT in so far as it does not rely (contrary to the PT) on an analytical expression for the conditional probability for the corresponding system without delay force. And even* if* this quantity is available, the respective quantity of the corresponding linearized, delayed system $\rho^{\mathrm{lin}}_{\mathrm{c,ss}}$ (which is the key ingredient within the FLC), often renders a better approximation, especially for large delay times. This is because $\rho^{\mathrm{lin}}_{\mathrm{c,ss}}$ is a non-Markovian quantity which can account for the finite correlations between the system states at the times $t$ and $t\!-\!\tau$. These correlations are not represented in the PT. Moreover, the FLC is not restricted to small noise levels, as the small delay expansion, or to small feedback forces, as the PT approach.

However, the FLC also has its limitations. It only works in the steady state, in contrast to the earlier approaches. Additionally, also the corresponding force linearized system must have steady state. Another issue occurs when our estimate for the total potential in the steady state $V_{\mathrm{STAT}}$ [Eq. (5)] has no minima giving rise to an ambiguity regarding the points around which the linearization should be carried out. However, one can easily adjust the ansatz of the FLC, by simply using other (e.g., random) centers of the linearization, preferring those where the probability density is high (in order to minimize the error in the final approximation). For instance, the minima of the static potential can be reasonable choices.

Having in mind these restrictions, the FLC can be applied to a much wider class of delayed systems than the ones considered in the present work. One important extension are systems with inertia, for which the Fokker-Planck description has been discussed extensively in Ref. Rosinberg et al. (2015). Such systems have been shown to exhibit complex dynamical behavior including chaos; an example is a delayed double-well potential with inertia and periodic forcing Nbendjo et al. (2003). Furthermore, it is in principle feasible to generalize the FLC towards systems with multiple discrete delays, which, for instance, can play a crucial role in the stochastic process of gene regulation Parmar et al. (2015); *Bratsun2005 or in the collective noisy motion of animals McKetterick and Giuggioli (2014). The force-linearized case of such systems is still solvable Giuggioli et al. (2016); McKetterick and Giuggioli (2014), with the $n$-time probabilities being again given by multivariate Gaussian distributions. Also for the more general class of noisy systems with distributed delays, the force linearized case can be handled analytically for arbitrary memory kernels Giuggioli et al. (2016); Hänggi and Talkner (1978); *Hanggi1982; *Hernandez1983; Budini and Cáceres (2004); McKetterick and Giuggioli (2014). This generally paves the path to adapt the basic idea of the FLC. However, to the best of our knowledge, the derivation of the FPE with multiple or distributed delays and nonlinear forces has not been carried out to date, at least not in the general case. A lack of this equation inherently limits the applicability of the FLC. Another important generalization would be the case of multiplicative noise (and delay), for which so far not many analytical results are known. It is clear that the situation will substantially depend on the specific functional form of the noise term. While at least the delayed Fokker-Planck equation has been worked out for the general case (independent of the specific multiplicative noise term) Guillouzic et al. (1999); Frank (2004), there is to the best of our knowledge no general solution for the linear force case at hand. For some specific multiplicative noise terms, however, the solution for the linear force model is indeed known Frank and Beek (2001); Frank (2004); Budini and Cáceres (2004), making the FLC in fact feasible (with non-Gaussian conditional probability distributions). Generally speaking, for all types of delayed stochastic systems, the applicability of the FLC depends on the availability of both, the general FPE and the solution of the force-linearized case.

A further aspect touched in the present work concerns the calculation of transport-related quantities such as escape times or waiting time distributions. In general, they require the two-time probability $\rho_{2}$ or even higher $n$-time densities. In the present form, the FLC does not provide access to $\rho_{n\geq 1}$. Still, we have demonstrated that a reasonable estimate of the escape times is possible by combining the FLC and the Kramers theory. This estimate, however, breaks down when the delay times are in the range of (or large compared to) the jump duration times. Under such conditions, the interplay between the dynamics of the system state variable $\chi$ and the time-dependent energy landscape causes oscillatory motion, whose description is clearly beyond the rather crude FLC-Kramers approximation. For this reason, an extension of the theory towards higher members of the FPE hierarchy would be worthwhile. The analytical results concerning the second member of the FPE hierarchy and $\rho_{3}$ for the linear case, given in the Appendices A-C, are essential steps in this direction.

A different way to access transport-related quantities was presented in Tsimring and Pikovsky (2001), where the autocorrelation function in a very similar model (a double-well potential supplemented with linear delay force) was successfully approximated by that of an appropriate discrete model (where an entire potential valley is regarded as one discrete state). In this case the Master equation can analytically be solved. We aim to stress that, in its present form, the FLC does not compete with this approach, but yields a description on a different time and length scale. While Ref. Tsimring and Pikovsky (2001) yields a convincing approximation of the inter-well dynamics, our approach rather targets the intra-well dynamics. Hence, to some extent, both approaches complement each other.

Finally, we want to briefly comment on the possible relevance of the FLC for stochastic thermodynamics. Indeed the thermodynamics of delayed systems is a subject of strong current interest Rosinberg et al. (2015); Jiang et al. (2011); *Munakata2009; *Munakata2014. Of particular interest are the information and entropy change rates, which rely on two-time probabilities. Time-delayed feedback forces have been shown to yield nontrivial contributions to these quantities Rosinberg et al. (2015); Xiao (2016). In its present form, the framework of the FLC does not provide access to the involved probabilities. This is another motivation to extend our approach towards higher $n$-time probability densities. One could then use the FLC predictions in order to estimate these intriguing thermodynamic quantities far from equilibrium.

Acknowledgements.

This work was supported by the Deutsche Forschungsgemeinschaft through SFB 910 (project B2).

Appendix A SECOND MEMBER OF THE DELAYED FPE, GENERAL CASE

In the following we present a derivation of the delayed FPE for $\rho_{2}$, which, to the best of our knowledge, has not been reported elsewhere so far. We extend the derivation presented in Ref. Frank (2005a) [on which we comment before Eqs. (4)], towards the second member of the infinite FPE hierarchy. Accordingly, we start by considering the formal time derivative of \rho_{2}(x,t;x_{\tau},t\!-\!\tau)\!=\!\big{\langle}\delta[x-\chi(t)]\delta[x_{\tau}-\chi(t\!-\!\tau)]\big{\rangle}, and use the delayed LE (2) to substitute the derivatives ${\partial\chi(t)}/{\partial t}$ and ${\partial\chi(t\!-\!\tau)}/{\partial t}$. We then introduce the functional $\Lambda[\Gamma]\equiv\delta(x-\chi(t))\delta(x_{\tau}-\chi(t-\tau))$, just like in Ref. Frank (2005a). Using Novikov’s theorem [see Eq. (3)], we can express the emerging correlations \big{\langle}\Lambda[\Gamma]\Gamma(t)\big{\rangle} by functional derivatives. We now use basic variational derivative rules:

[TABLE]

and analogously for ${\delta\Lambda[\Gamma]}/{\delta\Gamma(t-\tau)}$. The last term in Eq. (25) vanishes, since it is noncausal: The particle position at time $t-\tau$ cannot depend on the noise at the later time $t$, since $\Gamma$ is a random force with vanishing temporal correlation Frank et al. (2003). In the case of additive noise, the variational derivatives ${\delta\chi(t)}/{\delta\Gamma(t)}$ and ${\delta\chi(t)}/{\delta\Gamma(t-\tau)}$ do not explicitly depend on $\Gamma$, so we can evaluate the ensemble averages and obtain the FPE for $\rho_{2}$

[TABLE]

This is the second member of the Fokker-Planck hierarchy for the general case of a system with delay force $F(x,x_{\tau})$. As expected, it contains the next higher-order probability distribution function, that is $\rho_{3}$. We also note that, in contrast to the first member of the FPE hierarchy [see Eq. (4)], Eq. (A) contains the variational derivatives ${\delta\chi(t)}/{\delta\Gamma(t)}$ and ${\delta\chi(t)}/{\delta\Gamma(t-\tau)}$. These need to be specified for the specific $F$ under consideration.

Appendix B SECOND MEMBER OF FPE FOR LINEAR FORCES

Now we specialize the delayed FPE for $\rho_{2}$ [see Eq. (A)] for forces that are linear in $x$ and $x_{\tau}$, i.e., forces of the form $F^{\mathrm{L}}(x,x_{\tau})=-{\alpha x}-\beta x_{\tau}$. The corresponding delayed LE can be solved iteratively by using the method of steps Frank et al. (2003): to solve the equation in a time interval $t\in[n\tau,(n\!+\!1)\tau],\forall n\in\mathbb{N}$, one substitutes the solution of the preceding interval $[(n\!-\!1)\tau,n\tau]$ to get rid of the delay term (for $n\!=\!0$ the history function is used). One finds that on each interval, the formal solution has the same structure

[TABLE]

with $\chi(t)=\mathcal{Y}_{n+1}(t)$, for $t\in[n\tau,(n+1)\tau]$. From Eq. (B) we find the functional derivatives

[TABLE]

Inserting Eqs. (28) and (29) into Eq. (A) we obtain for the case of linear forces $F^{\mathrm{L}}(x,x_{\tau})=-{\alpha x}-\beta x_{\tau}$:

[TABLE]

In the steady state, Eqs. (B) and the FPE for $\rho_{1}$ [Eq. (4)] are solved by the respective (multivariate) Gaussian distributions Küchler and Mensch (1992)

[TABLE]

with $\bm{x}\!=\!(x,x_{\tau},..,x_{(n-1)\tau})^{\mathrm{T}}$, $\mathbf{D_{1}}\!=\!C(0)$, and

[TABLE]

The covariance matrices $\mathbf{D_{n}}$ involve the steady state spatial autocorrelation function between the systems state at time $t$ and time $t+z$:

[TABLE]

Explicit expressions for $C(0)$, $C(\tau)$, and $C(2\tau)$ (which occur in $\mathbf{D_{k}}$, k$=1,2,3$) are given in the next section.

Appendix C SPATIAL AUTOCORRELATION FUNCTION

Due to the symmetry property $C(z)=C(-z)$, we only consider non-negative arguments $z$ in the following. As shown in Refs. Guillouzic et al. (1999); Frank et al. (2003), the spatial autocorrelation function on the interval $z\!\in\![0,\tau]$ is given by

[TABLE]

with $\omega\!=\!\sqrt{\!\alpha^{2}-\beta^{2}}/\gamma\in\!\mathbb{C}$. In order to obtain $C(z)$ for $z>0$, one can directly deduce from the delayed LE Frank et al. (2003):

[TABLE]

For $z>0$, the last term vanishes, since a non-zero correlation between the system state $\chi$ and the future noise would be noncausal. The resulting differential equation for $z>0$

[TABLE]

can be solved iteratively using the method of steps (see text at the beginning of the Appendix B). Using the solution for the interval $z\in[0,\tau]$ in Eqs. (34) and (35), we obtain the correlation function on the interval $z\in[\tau,2\tau]$,

[TABLE]

where

[TABLE]

Inspection of Eq. (38) shows that there exist critical $\tau$-values for which $C(z)$ diverges, see also text after Eq. (8d) (with $\alpha\!=\!\alpha_{i}$ and $\beta\!=\!\beta_{i}$ in the discussion there).

We would like to note that in Ref. Frank et al. (2003), the solution (35) is erroneously stated to be valid for all $z\!\in\!\mathbb{R}$. Along their derivation, Eq. (18) of Ref. Frank et al. (2003) is actually not valid for $|z|\!>\!\tau$, although they claim otherwise. A correlation function (35) $\forall z\!\in\!\mathbb{R}$ would in fact be unphysical. Consider for example $\beta\!<\!\alpha$, where the system is spatially confined. Then, it is not reasonable that the autocorrelation function grows without bounds for large time differences, as Eq. (35) would predict. The iterative solution obtained from Eq. (37) by using the method of steps, is indeed bounded, consisting with the physical intuition.

Appendix D BROWNIAN DYNAMICS SIMULATIONS

In order to test our theoretical results against quasi-exact data, we perform Brownian dynamics (BD) simulations of the delayed LE (2). We use the Euler-Maruyama integration scheme Buckwar (2000); Kloeden and Platen (1992), with a varying temporal discretization of $\Delta t\!\in\![10^{-3},10^{-6}]\,\tau_{\mathrm{B}}$, such that $1000\,\Delta t\leq\tau$. We perform each simulation multiple times (with different random number seed), and take the average over at least $10^{5}$ realizations. The (pseudo) random numbers are generated with the algorithm “Mersenne Twister” Mer (2016) and the Box-Müller method Kloeden and Platen (1992).

As initial condition, we use an equilibrium configuration in the respective static potential without delay force. Before we measure the steady state properties, we let $100\,\tau_{\mathrm{B}}$ times pass in the presence of the delay force to let the system reach a steady state. The simulation time thereafter is more than $200\,\tau_{\mathrm{B}}$.

The density profiles from the BD simulations are obtained from histograms with a spatial resolution of $\Delta x\!=\!0.005\,\sigma$. These histograms are also used to calculate the moments $\mu_{n}$. With the same procedure the moments of the analytical distributions resulting from the FLC, PT and small delay expansion are calculated (we use the same bin sizes and positions as in the simulations).