Fokker-Planck equations for time-delayed systems via Markovian Embedding

TL;DR

This paper introduces a novel Markovian embedding method to derive a complete Fokker-Planck equation for time-delayed stochastic systems, providing new insights and improved approximations for the two-time probability density function.

Contribution

It presents a new approach to obtain a closed Fokker-Planck equation for delayed systems using Markovian embedding, simplifying previous methods and offering better approximations for the two-time PDF.

Findings

Derived a closed, infinite-dimensional FPE for delayed systems.

Proposed an effective approximation for the two-time PDF.

Outperformed earlier methods in large delay regimes.

Abstract

For stochastic systems with discrete time delay, the Fokker-Planck equation (FPE) of the one-time probability density function (PDF) does not provide a complete, self-contained probabilistic description. It explicitly involves the two-time PDF, and represents, in fact, only the first member of an infinite hierarchy. We here introduce a new approach to obtain a Fokker-Planck description by using a Markovian embedding technique and a subsequent limiting procedure. On this way, we find a closed, complete FPE in an infinite-dimensional space, from which one can derive a hierarchy of FPEs. While the first member is the well-known FPE for the one-time PDF, we obtain, as second member, a new representation of the equation for the two-time PDF. From a conceptual point of view, our approach is simpler than earlier derivations and it yields interesting insight into both, the physical meaning, and the mathematical structure of delayed processes. We further propose an approximation for the two-time PDF, which is a central quantity in the description of these non-Markovian systems as it directly gives the correlation between the present and the delayed state. Application to a prototypical bistable system reveals that this approximation captures the non-trivial effects induced by the delay remarkably well, despite its surprisingly simple form. Moreover, it outperforms earlier approaches for the one-time PDF in the regime of large delays.

Equations

$$\mathrm{d}{X}(t)={F}[X(t),{X(t-\tau)}]\mathrm{d}t+{\color[rgb]{0,0,0}\sqrt{2D_{0}}\,\mathrm{d}W(t)},$$

$${\partial }_t{\rho }_1(x,t)=-{\partial }_x{\int }_{\Omega }\left[F(x,x_{\tau }){\rho }_2(x,t;x_{\tau },t-\tau )\right]\mathrm{d}{x}_{\tau }+D_0{\partial }_x^2{\rho }_1(x,t),$$

$${\partial }_t{\rho }_2(x,t;x_{\tau },t-\tau )=$$

$$\displaystyle+D_{0}\,\Bigg{(}{\partial_{x}^{2}}+{\partial_{x_{\tau}}^{2}}+{\partial_{x}\partial_{x_{\tau}}}\frac{\delta X(t)}{\delta\xi(t\!-\!\tau)}\Big{|}_{\!\!\!\!\!\!\!\!\!\!X(t)=x\atop\!\!X(t-\tau)=x_{\tau}}\Bigg{)}\rho_{2}(x,t;x_{\tau},t\!-\!\tau)$$

$$\mathrm{d}{X}_{j\tau }(z)=F[X_{j\tau }(z),X_{(j-1)\tau }(z)]\mathrm{d}z+\sqrt{2D_0}\mathrm{d}{W}_{j\tau }(z).$$

$${\partial }_z{\rho }_{N+1}=$$

$${\partial }_t{\rho }_2=-{\partial }_x\left[F(x,x_{\tau }){\rho }_2\right]-{\partial }_{x_{\tau }}{\int }_{\Omega }\left[F(x,x_{\tau }){\rho }_3\right]\mathrm{d}{x}_{2\tau }+D_0\left[{\partial }_x^2+{\partial }_{x_{\tau }}^2\right]{\rho }_2,$$

$${\dot{X}}_0(t)=$$

$${\dot{X}}_j(t)=$$

$$\dot{X}_{0}(t)=F\left[X_{0}(t),\int\limits_{-\tau}^{t}\underbrace{\left(\frac{n}{\tau}\right)^{n}\frac{(t-s)^{n-1}}{(n-1)!}\,e^{-n(t-s)/\tau}}_{=K_{n}}X_{0}(s)\,\mathrm{d}s\right]+{\color[rgb]{0,0,0}\sqrt{2D_{0}}\,\,\xi(t)},$$

$$X_0(t)=X(t),X_n(\tau )=X(t-\tau ).$$

$$X_j(t)=X(t-j\tau /n),\forall j\in \{0,1,2,..,n\}$$

$$\displaystyle\underbrace{\partial_{t}\rho_{n+1}{\color[rgb]{1,1,1}\Bigg{|}}}_{\mathrm{(I)}}=\underbrace{-\partial_{x_{0}}\big{[}F(x_{0},x_{n})\rho_{n+1}\big{]}{\color[rgb]{1,1,1}\Bigg{|}}}_{\mathrm{(II)}}\underbrace{-\frac{n}{\tau}\sum_{j=1}^{n}\partial_{x_{j}}\left[(x_{j-1}-x_{j})\rho_{n+1}\right]}_{\mathrm{(IV)}}\underbrace{+D_{0}\,\partial_{x_{0}}^{2}\rho_{n+1}{\color[rgb]{1,1,1}\Bigg{|}}}_{\mathrm{(III)}}$$

$$\underset{x_j\to {\omega }_{1,2}}{\lim }{\rho }_m(x_j,...)=\underset{x_j\to {\omega }_{1,2}}{\lim }{\partial }_{x_j}{\rho }_m(x_j,...)=0,\forall j,m,$$

$$(\mathrm{I})$$

$$(\mathrm{IV})\to$$

$$\displaystyle~{}-\int_{\Omega}x_{1}\,\partial_{x_{1}}\rho_{2}(x_{0};x_{1})\mathrm{d}x_{1}+x_{0}\underbrace{\left[\rho_{2}(x_{0};x_{1})\right]_{\omega_{1}}^{\omega_{2}}}_{=0}\Bigg{\}}$$

$$\displaystyle=\frac{n}{\tau}\Bigg{\{}n\,\rho_{1}(x_{0})+\sum_{j=2}^{n}\underbrace{\left[x_{j}\rho_{3}(x_{0};x_{j-1};x_{j})\right]_{\omega_{1}}^{\omega_{2}}}_{=0}-\sum_{j=2}^{n}\int_{\Omega}\rho_{2}(x_{0};x_{j})\mathrm{d}x_{j}-\underbrace{\left[x_{1}\rho_{2}(x_{0};x_{1})\right]_{\omega_{1}}^{\omega_{2}}}_{=0}-\rho_{1}(x_{0})\Bigg{\}}=0.$$

$${\partial }_t{\rho }_1(x_0,t)=-{\partial }_{x_0}{\int }_{\Omega }F(x_0,x_n){\rho }_2(x_0;x_n)\mathrm{d}{x}_n+D_0{\partial }_{x_0}^2{\rho }_1(x_0,t).$$

$${\partial }_t{\rho }_2(x_0,t;x_n,t)=$$

$$+{\partial }_{x_n}{\int }_{\Omega }\frac{n}{\tau }(x_n-x_{n-1}){\rho }_3(x_0,t;x_{n-1},t;x_n,t)\mathrm{d}{x}_{n-1},$$

$${\partial }_t{\rho }_m(x_0;x_{n-m+2},..,x_n)=$$

$$+\sum _{l=0}^{m-3}\frac{n}{\tau }{\partial }_{x_{n-l}}\left[(x_{n-l}-x_{n-l-1}){\rho }_m(x_0;x_{n-m+2},..,x_n)\right]$$

$$+\frac{n}{\tau }{\partial }_{x_{n-m+2}}{\int }_{\Omega }(x_{n-m+2}-x_{n-m+1}){\rho }_{m+1}(x_0;x_{n-m+1},..,x_n)\mathrm{d}{x}_{n-m+1}$$

$$+D_0{\partial }_{x_0}^2{\rho }_m(x_0,t;x_{n-m+2},t,..,x_n,t),$$

$$\frac{n}{\tau }{\partial }_{x_n}{\int }_{\Omega }(x_n-x_{n-1})$$

$${\iiint }_{\Omega }{\rho }_3\left(x_0,t;x_{n-1},t;x_n,t\right)\mathrm{d}{x}_0\mathrm{d}{x}_{n-1}\mathrm{d}{x}_n=\frac{\tau }{n}{\int }_{\Omega }{\int }_{\tilde{\Omega }}{\int }_{\Omega }{\rho }_3\left(x_0,t;\tilde{x}+x_n,t;x_n,t\right)\mathrm{d}{x}_0\mathrm{d}\tilde{x}\mathrm{d}{x}_n$$

$${\rho }_3\left(x_0,t;\frac{\tau }{n}\tilde{x}+x_n,t;x_n,t\right)\stackrel{!}{=}\frac{\tau }{n}\widetilde{{\rho }_3}\left(x_0,t;\tilde{x},t;x_n,t\right).$$

$$\frac{n}{\tau }{\partial }_{x_n}{\int }_{\Omega }(x_n-x_{n-1}){\rho }_3(x_0;x_{n-1};x_n)\mathrm{d}{x}_{n-1}⟶$$

$$=$$

$$\{x_0,t\}$$

$$\{x_n,t\}$$

$$\{\tilde{x},t\}$$

$${\partial }_{x_n}{\int }_{\tilde{\Omega }}\tilde{x}\widetilde{{\rho }_3}\left(x_0,t;\tilde{x},t;x_n,t\right)\mathrm{d}\tilde{x}$$

$${\partial }_t{\rho }_2(x,t;x_{\tau },t-\tau )=$$

$$+{\partial }_{x_{\tau }}\left[⟨\dot{X}(t-\tau )|x=X(t),x_{\tau }=X(t-\tau )⟩{\rho }_2(x,t;x_{\tau },t-\tau )\right]$$

$$-⟨\dot{X}⟩{\partial }_{x_{\tau }}{\rho }_2(x,t;x_{\tau },t-\tau )=$$

$${\rho }_2(x,t;x_{\tau },t-\tau )=Z^{-1}h(x_{\tau })e^{\int F(x,x_{\tau })\mathrm{d}x},$$

$${\rho }_2^{\mathrm{approx}}(x,t;x_{\tau },t-\tau )=Z^{-1}{e}^{\int F_{\mathrm{nd}}(x)\mathrm{d}x+kx{x}_{\tau }+\int F_{\mathrm{nd}}(x_{\tau })\mathrm{d}{x}_{\tau }},$$

Full text

1 Introduction

Predicting stochastic dynamics far from thermal equilibrium is a major goal of statistical physics and a key step towards an understanding of realistic, fluctuation-dominated systems from biology, chemistry, or socio-economics, just to name a few of the many fields of application. An important mechanism that drives system out of equilibrium and, at the same time, renders the dynamics non-Markovian is time delay. Delays can have various origins, such as unavoidable signal transmission time lags (e. g., in laser systems [1, 2]), information processing times (as in neural systems [3, 4, 5, 6]), or decision making times (e. g., in the context of financial markets [7, 8]). Moreover, they are omnipresent in systems subject to feedback control [9, 10, 11, 12, 13, 14, 15, 16, 17, 18]. In many cases, it is sufficient and suitable to model such processes by including a single, discrete delay and a white noise. This is the situation considered in this paper. Discrete delays are known to induce complex behavior (e. g., delay-induced oscillations, chaos), and, on the other hand, to stabilize unstable orbits [19, 20, 21]. Therefore, they are often implemented intentionally, for example in the framework of Pyragas control [22]. It is important to note that the time-delayed systems considered here are different from those modeled by generalized Langevin equations [23]. The latter involve a memory kernel in the deterministic (friction) force resulting from conservative interactions between a colloidal particle (or, more generally, a “slow” variable) and the bath particles (as well as among the bath particles). This delayed friction is related to the noise via a fluctuation-dissipation relation [24] (see, e. g., [25, 26, 27, 28]). In contrast, the systems considered here involve a discrete delay within a deterministic force, and there is no relation to noise correlations.

While they have numerous applications, the mathematical description of time-delayed stochastic systems is still far from being complete. The Langevin equation (LE) is used to model such processes on a stochastic level, but it is infinite-dimensional due to the discrete delay. Exact analytical results based on the LE have, so far, been obtained only for linear systems [29]. On the other hand, the Fokker-Planck (FP) description for the probability density function (PDF), which we focus on in this paper, is given by an infinite hierarchy of coupled equations of increasing complexity [30, 31, 32, 33]. As shown by Frank [30] (who only considered the first member), this hierarchy can be derived using Novikov’s theorem [34]. The first member of the FP hierarchy (originally derived by Guillouzic et al. [35]) is the FP equation (FPE) for the one-time PDF. It explicitly involves the two-time PDF for the spatio-temporal correlations between times $t$ and $t-\tau$, with $\tau$ being the delay time. Although not being self-sufficient, the equation for the one-time PDF has proven itself to be quite useful for the analytical treatment. On the one hand, it is an important tool in the search for exact results [36, 37]. Moreover, on the base of this FPE different approximation schemes have been developed in the past years (see [31] for an overview). Examples are the perturbation theory [38] and the “force-linearization closure” [31, 37], both of which yield reasonable approximations for the one-time PDFs in quite large parameter regimes. (As opposed to this, the so called small $\tau$ expansion [35], which is based on the LE, is only appropriate for very small values of the delay time, where the system behaves more or less like a Markovian one.)

While these analytical results have indeed contributed to an understanding of delayed processes, major challenges remain. In particular, the insights provided by the one-time PDF alone are admittedly limited as we deal with non-Markovian dynamics, which is crucially determined by memory, and, hence, temporal correlations. Having this in mind, it is not surprising that many non-Markovian effects are only captured by two-time PDFs (or higher PDFs). In fact, only by studying them, non-Markovian steady-states can be qualitatively distinguished from thermal equilibria at all111As an illustrative example, consider the following. Let $\rho^{\mathrm{nM}}_{1}(x)$ be the non-Markovian one-time PDF of an arbitrary nonequilibrium steady state. This PDF can in general not be distinguished from a PDF $\rho^{\mathrm{M}}_{1}(x)$ of a Markovian (equilibrium) process in a (fictive) potential defined by $U_{\mathrm{ficitve}}(x)\equiv k_{\mathrm{B}}\mathcal{T}\ln{\rho^{\mathrm{nM}}_{1}(x)}$, which would equally yield a distribution $\rho^{\mathrm{M}}_{1}=Z^{-1}e^{-U_{\mathrm{fictive}}/k_{\mathrm{B}}\mathcal{T}}=\rho^{\mathrm{nM}}_{1}$.. Moreover, the two-time PDF is essential to obtain stochastic thermodynamical quantities like the fluctuating heat, work, or entropy [39]. As an example, the $m^{\mathrm{th}}$ moment of the heat rate [40, 37] of a system described by the LE (1) is given by $\langle\dot{q}(t)^{m}\rangle=\langle F[X(t),X(t-\tau)]^{2m}\rangle+const.$, which means that all moments can be evaluated via the two-time PDF. Furthermore, calculations of transport properties such as, e. g., the Kramers rates, that are solely based on the one-time PDF have been shown to be limited to parameter regimes where non-Markovian effects are rather weak [31]. However, despite the clear benefits of a probabilistic treatment and the apparent limits of a study focusing on the first member of the FP hierarchy alone, no general approximation schemes for the two-time PDF or other temporal correlation functions are known (one exception being the two-state reduction employed in [41, 37]). In fact, the higher FP hierarchy members have rarely been discussed at all. This hinders the development of further approximation schemes and solution methods. In the present work we concentrate on these higher members. By introducing a different FP approach which yields an alternative representation of the hierarchy, we provide a step towards a full understanding of the probabilistic description.

Contrary to earlier approaches, we here derive a FP hierarchy by using a Markovian embedding technique [42, 43, 44, 45, 26, 46], which is inspired by the linear chain trick [47, 48, 49]. In particular, we introduce $n\in\mathbb{N}$ auxiliary variables whose dynamics is given by a set of $n$ coupled, linear, stochastic differential equations without delay. This set generates the same stochastic process as the delayed equation in the limit $n\to\infty$. The resulting Markovian system is described by a closed, self-sufficient FPE on an $n$-dimensional space. By marginalization, i. e., projection onto a one-dimensional subspace, and subsequent limiting procedure, we obtain the well-known delayed FPE for the one-time PDF. Furthermore, by projection onto higher-dimensional subspaces, we derive a FP hierarchy which is different from the one obtained by Novikov’s theorem [34]. Interestingly, the limiting procedure contains one step where position-like variables are converted into velocity-like variables. Explicitly considering the second member, we show how the two hierarchies can be converted into each other.

Furthermore, we propose an approximation for the steady-state probability density in the regime of large delay times, based on the new equation for the two-time PDF. More specifically, by neglecting correlations between the system state and the displacement, we obtain a closed equation that can be solved analytically. To test the approximation, we consider a bistable system involving a doublewell potential and a linear delay force. This is an established prototypical model to study the combined impact of noise and time delay on nonlinear dynamics [41, 37, 31]. In particular, the model exhibits delay-induced particle oscillations between the two minima. As these oscillations manifest in the two-time PDF, this system is an ideal candidate to test approximations, and to study non-Markovian effects on the level of the two-time PDF. By comparing with numerical simulations, we demonstrate that the proposed approximation yields surprisingly good results. In particular, the approximation captures the main characteristics of the two-time PDF, both, when the delay force is directed towards, and when it points away from the delayed position. In this way, we present to our knowledge the first approximation for the two-time probability density of stochastic systems with time delay. In the regime of large delay times, the new approximation even outperforms earlier approaches known from the literature, demonstrating that a framework on the level of the second FP member, is somewhat superior compared to approaches on the first member.

This paper is organized as follows. First, we review earlier derivations of the FPE known from the literature. Then, we discuss the Markovian embedding and show how the delayed FPE can be derived from it by marginalization. We further discuss the higher members stemming from this approach, particularly focusing on the second member, which are different from those obtained by earlier approaches. At the end of the main text, we introduce the aforementioned approximation scheme and discuss the application to the bistable, delayed system. In the appendix, we show how the second member of our hierarchy can be transformed into the one obtained from Novikov’s theorem, and present a short derivation of Novikov’s theorem via path integrals. In this work, we restrict ourselves to systems in one spatial dimension with a single time delay. However, our results can be readily generalized to multidimensional systems, and multiple delays can, as well, be treated analogously.

2 Langevin equation and Fokker-Planck equation for delayed systems

We consider a system described by a stochastic delay differential equation, in particular, by the overdamped LE

[TABLE]

with $X\in\Omega=[\omega_{1},\omega_{2}]\subseteq\mathbb{R}$ and time $t>0$. $F$ denotes the deterministic force222The force is scaled with the friction coefficient $\gamma$, i. e., $F\to F/\gamma$, which is related with the bath’s temperature $\mathcal{T}$ via $\gamma D_{0}=k_{\mathrm{B}}\mathcal{T}$, where $k_{\mathrm{B}}$ is the Boltzmann constant. given by some (generally nonlinear) function depending on the instantaneous and on the delayed system state, $X(t)$ and $X(t-\tau)$, with $\tau\geq 0$ being the single discrete time delay. $W$ is a Wiener process with independent increments, $\langle W\rangle=0$, and $\langle W(t)W(t^{\prime})\rangle=\min(t,t^{\prime})$, generated by the additive Gaussian white noise $\xi$, with $\langle\xi\rangle=0$, $\langle\xi(t)\xi(t^{\prime})\rangle=\delta(t-t^{\prime})$. Thereby, $\langle...\rangle$ denotes the ensemble average (i. e., average over all noise realizations), w. r. t. a given initial condition. The latter is specified by a history function $\phi(\tau\leq t\leq 0)$, i. e., $X(t)=\phi(t),t\in[-\tau,0]$, which can be fixed for the entire ensemble or be drawn from a distribution $P(\phi)$. $D_{0}$ is the (constant) strength of the noise. We further consider natural boundary conditions (details are given below).

The FP description provides a complementary way to model the process on the probabilistic level. It is well-known that the FPE for the one-time PDF related to (1) reads [50, 35, 30, 31]

[TABLE]

where $\rho_{N}(x,t;x_{\tau},t-\tau;...;x_{N\tau},t-N\tau)=\langle\delta(x=X(t))\delta(x_{\tau}=X(t-\tau))...\delta(x_{N\tau}=X(t-N\tau))\rangle$ denotes the $N$-time PDF (with normalization $\iint_{\Omega}...\int_{\Omega}\rho_{N}\mathrm{d}x\mathrm{d}x_{\tau}...\mathrm{d}x_{N\tau}=1$). As opposed to Markovian dynamics, the FPE for $\rho_{1}$ is not closed but involves the two-time (joint) PDF, $\rho_{2}$. The FPE (2) can be derived by different strategies [35, 30, 38, 50, 51]. However, only the approach which bases on Novikov’s theorem, results in a complete FP description (in the form of an infinite hierarchy, see below), while all the others only yield the first, not self-sufficient member. In the following, we briefly review earlier work, in particular Refs. [30, 50, 51], in order to better understand the technical and conceptual difference of our approach.

2.1 Second member of the Fokker-Planck hierarchy from Novikov’s theorem

In this section, we briefly review the approach [30, 31] based on Novikov’s theorem (28) (a short derivation of the latter is given in the Appendix C). Besides leading to (2), this approach can also be used to obtain the FPE for $\rho_{2}$ and, in general, $\rho_{N},~{}N\in\mathbb{N}$. It renders an infinite hierarchy of coupled equations of increasing complexity, whose $N$th member involves both, the $N$-time (joint) PDF $\rho_{N}(x,t;x_{\tau},t-\tau;x_{2\tau},t-2\tau;..;x_{N},t-N\tau)$ and the $(N+1)$-time PDF $\rho_{N+1}$. It should be noted, however, that the higher members (starting from the second one) contain functional derivatives w. r. t the noise, which need to be obtained separately from the LE. For example, the second member reads

[TABLE]

(see [31] for a derivation). This equation resembles a Markovian FPE of a two-dimensional (or two-particle) process in the variables $x$ and $x_{\tau}$, with the exception of the term which involves a derivative w. r. t. both variables and the functional derivative w. r. t. the noise. For the special case of linear forces $F$, the latter can be calculated using the method of steps [52], see Appendix B. However, to the best of our knowledge, this functional derivative is unknown for general nonlinear forces. We note that this term stems from the fact that if $X(t)$ and $X(t-\tau)$ are regarded as “two” random processes, then these processes are not only connected by the force $F$ (via the LE), but additionally by the fact that the corresponding noise processes $\xi(t)$ and $\xi(t-\tau)$ are identical with time-shift $\tau$. To elucidate this last point, we briefly review an alternative derivation of the FPE known from the literature, which does not yield this term, but instead imposes additional constraints on the solutions.

The mentioned approach is used in [50, 51] to derive a FPE for $\rho_{1}$. Here we extended it towards the higher members of the hierarchy. The approach uses a description with two time arguments $j$, $z$. In particular, an integer $j=0,1,2,..$ counts the number of intervals of length $\tau$ that have passed since the beginning $t=0$, and a continuous variable $z\in[0,\tau]$ measures the time within the current interval. Then, auxiliary phase-space variables $X_{j\tau}(z):=X(z+j\tau)$ are introduced, which, by construction, all follow a LE of the identical form

[TABLE]

(The steady-state dynamics can be accessed if $j\to\infty$.) A sketch of the construction of the auxiliary variables is given in Fig. 1. In the probability space belonging to the phase space $\{X_{0\tau}(z),...,X_{N\tau}(z)\}$, the process is Markovian, and the corresponding FPE is given by

[TABLE]

with $\rho_{N+1}(x_{0\tau},x_{1\tau},...,x_{N\tau},z)=\langle\delta[x_{0\tau}-X(z)]\delta[x_{1\tau}-X(z+\tau)]...\delta[x_{N\tau}-X(z+N\tau)]\rangle$. By integrating the equation over the entire domain of all variables but the first (and the last) one, this many-variable FPE can be used to derive the FPE (2) for $\rho_{1}$, and a FPE for $\rho_{2}$, reading

[TABLE]

with $\rho_{2}(x,t;x_{\tau},t-\tau)$ and $\rho_{3}(x,t;x_{\tau},t-\tau;x_{2\tau},t-2\tau)$. As opposed to the corresponding equation from Novikov’s theorem, Eq. (2.1), this equation lacks of the term with the unknown functional derivative. This is because the stochastic set of Eqs. (4) is, in fact, only equivalent to the delayed process, if one further imposes additional constraints $X_{j}(\tau){=}X_{j+1}(0),~{}\forall j$, ensuring that the end point of the stochastic trajectories on $[j\tau,(j+1)\tau]$ matches the starting point on $[(j+1)\tau,(j+2)\tau]$. Translating these constraints to the level of PDFs is a non-trivial task on its own, thus, the usefulness of these equations as a starting point for approximation schemes is very limited. The comparison between Eq. (2.1) from Novikov’s theorem and (6) suggests that the terms with the unknown functional derivatives in (2.1) are indeed such terms ensuring these additional constraints on the trajectories.

3 Markovian embedding

As an alternative approach to formulate a complete FP description and to derive the FPE (2), we here employ the Markovian embedding technique [42, 43, 44, 45, 26, 46], as also suggested in earlier literature [49]. The key idea is to introduce $n\in\mathbb{N}$ auxiliary variables $X_{j}$ (with one of them replacing the delayed system state in the original equation), giving rise to a $(n+1)$-dimensional Markovian system that generates the same dynamics as the original delay equation (which is hence “embedded” in that extended system). Specifically, we consider the set of dynamical equations given by [47]

[TABLE]

with $j\in\{1,2,...,n\}$, $t\geq-\tau$, $X_{j}\in\mathbb{R}$ and the initial conditions $X_{j\in\{1,..,n\}}(-\tau)\equiv 0$ (see Appendix A for a discussion of the initial conditions). Projection of (7b) onto the variable $X_{0}$ (see Appendix A for details) yields

[TABLE]

with the Gamma-distributed memory kernel $K_{n}$. In the limit $n\rightarrow\infty$, the kernel collapses onto a delta peak at $\tau$, i. e., $K_{n}(t-s)\to\delta(t-s-\tau)$. Then, the projected Eq. (8) becomes identical to the original, delayed LE (1), with

[TABLE]

In addition to $X_{0}$ and $X_{n}$, also the other auxiliary variables can be linked to the variable $X$ of the delayed LE (1). In particular, using the difference quotient to approximate $\dot{X}_{j}$ in the dynamical Eqs. (7b), that is, $[X_{j-1}(t)-X_{j}(t)](n/\tau)=\dot{X}_{j}(t)\approx[X_{j}(t+\Delta t)-X_{j}(t)]/\Delta t$, gives rise to the interpretations $\Delta t=\tau/n$, and $X_{j}(t+\Delta t)\approx X_{j-1}(t)$. Iteration implies $X_{j}(t+j\tau/n)\approx X_{j-i}(t)$. In combination with Eq. (9), this reveals

[TABLE]

which becomes accurate in the limit $n\to\infty$ (at $\Delta t\to 0$ the difference quotient equals the differential quotient).

The here presented scheme of rewriting the delayed differential equation as a set of non-delayed ones is called the linear chain trick [47, 48, 49]. A schematic visualization of the system described by (7) is given in Fig. 2. Importantly, the $n$ auxiliary variables generally follow linear Eqs. (7b), even if the delayed process itself involves nonlinear forces, $F$. Furthermore, their dynamical equations are deterministic, as opposed to $X_{0}$ which is subject to the white noise term $\xi$ [see (7a)]. We stress that the set of auxiliary variables (10) is different from the one (4) used in the earlier approach, which we have reviewed in the preceding section. A direct comparison is provided in Fig. 1.

Next, we will use the Markovian embedding to find a probabilistic description of the delayed system.

3.1 Fokker-Planck equation of the Markovian $(n+1)$-dimensional system

The Markovian $(n+1)$-dimensional system given by the set of coupled LEs (7), can also be described by the FPE

[TABLE]

for the $(n+1)$-point (joint) PDF $\rho_{n+1}(x_{0},t;x_{1},t;...;x_{n},t)=\langle\delta[x_{0}-X_{0}(t)]\delta[x_{1}-X_{1}(t)]..\delta[x_{n}-X_{n}(t)]\rangle$. Equation (11) is a self-sufficient FPE, which describes the dynamics of the delayed process in a closed form in the limit $n\to\infty$.

Importantly, $n$ will always be considered a finite value here (until after marginalization). Finding a mathematical framework in which the limit $n\to\infty$ of Eq. (11) is meaningful and well-defined, is subject of ongoing work. Indeed, it is not clear how to formulate a Fokker-Planck equation on an infinite-dimensional space [49]. Moreover, one drift term is in fact proportional to $n$ itself, making this limit highly non-trivial. Here, we can safely ignore this problem, since we focus on the equations for the marginalized (low-dimensional) PDFs. As shown in Secs. 4 and 5, they can be derived by keeping $n$ finite and performing the limit only after the marginalization.

Equation (11) resembles the FPE of a many-particle system. However, in sharp contrast to a system of interacting (colloidal) particles, the “interactions” are here unidirectional (see Fig. 2), which means that the “particle” $x_{j}$ only feels the “particle” $x_{j-1}$ and *not *vice-versa. Hence, the forces between them are non-Newtonian (violating Newtons 3rd law) and cannot be described by an interaction potential.

4 Derivation of first member of Fokker-Planck hierarchy via Markovian embedding

4.1 Marginalization

In the following, we show that the delayed FPE for $\rho_{1}$ can be derived by a “coarse-graining” of (11). To be more specific, we marginalize the FPE (11) for the joint PDF, that is we integrate it over the whole domain of all $n$ auxiliary variables, i. e., $\iiint_{\Omega}\mathrm{d}x_{1}\mathrm{d}x_{2}..\mathrm{d}x_{n}$, and, therewith obtain an equation for the marginal PDFs333We use the general relation between joint, $\rho_{2}$, and marginal distributions $\rho_{1}(z)=\int_{a}\rho_{2}(y,z)\mathrm{d}y$, with the normalization $\int_{b}\rho_{1}(z)\mathrm{d}z=\int_{b}\int_{a}\rho_{2}(y,z)\mathrm{d}y\mathrm{d}z=1$, with $y\in a$, $z\in b$. . Importantly, we take the limit $n\rightarrow\infty$ only after the marginalization. We further use natural boundary conditions at the endpoints of the domain $\Omega=[\omega_{1},\omega_{2}]$ specified by

[TABLE]

which are justified in most physical scenarios. For sake of a shorter notation, we will omit the time arguments of the PDFs in lengthy expressions. The terms (I – III) in Eq. (11) can be readily marginalized to

[TABLE]

[TABLE]

So far, the marginalized FPE (13) still pertains to a finite (but arbitrary) value of $n\in\mathbb{N}$. As a last step, we take the limit. To this end, we recall the meaning of the respective stochastic variables in this limit [Eq. (9)]. Based on that we perform the phase space transformation $\{x_{0},t\}\rightarrow\{x,t\}$ and $\{x_{n},t\}\rightarrow\{x_{\tau},t-\tau\}$. Therewith, we finally obtain the delayed FPE for the one-time PDF (2).

Next, we demonstrate that our procedure can also be applied to find higher members of the FP hierarchy.

5 Derivation of higher members via Markovian embedding

In the following, we derive a FPE for the two-time PDF $\rho_{2}(x,t;x_{\tau},t-\tau)$. This is a particularly important quantity in the description of delayed processes, as it explicitly describes the correlation between the delayed and present state. This time, we marginalize the complete Fokker-Planck equation w. r. t. all auxiliary variables but the last one (again keeping $n$ finite until after the marginalization). We do not integrate over $x_{n}$ as this variable represents $X(t-\tau)$ in the limit $n\to\infty$. In a second step, we again perform a coordinate transformation to take the limit.

5.1 Marginalization

We again start from the FPE (11) of the $n$-variable system and integrate the equation over $x_{1},..,x_{n}$. Using similar steps as presented in Sec. 4, e. g., partial integration and subsequent application of the boundary conditions, we obtain the marginalized FPE

[TABLE]

which involves a three-point PDF. Please note that integrating (5.1) over $x_{n}\in\Omega$ yields the first member of the hierarchy (13), which shows the consistency of our calculation. Higher-order equations can be obtained in an analogous manner, yielding the infinite hierarchy of coupled equations

[TABLE]

for $m\geq 3$. Note that in contrast to the FPE (13) for $\rho_{1}$, all higher members of the FP hierarchy (5.1,5.1) explicitly contain the linear forces stemming from the auxiliary variable Eqs. (7b).

We now inspect in more detail the second member (5.1) in the limit $n\to\infty$, to recover the system with discrete delay. We again need to transform the probability space, hence, the phase space and the PDFs defined on it.

5.2 Transformation of phase space, limiting procedure

In order to take the limit $n\to\infty$ of Eq. (5.1), we will transform the phase space and the PDFs defined on it. We perform the transformation in two steps. First we replace the phase space variable $x_{n-1}$ by $\widetilde{x}=(n/{\tau})(x_{n-1}-x_{n})$, which will later in the limit $n\to\infty$ receive the meaning of a velocity. In a second step, we transform back to the original variables of the delayed LE (1), whereby mainly the time arguments will be affected.

As first step, we change the phase space from $\{x_{0},t;x_{n-1},t;x_{n},t\}$ to $\{x_{0},t;\widetilde{x},t;x_{n},t\}$. To this end, we start with performing the substitution $x_{n-1}~{}\longrightarrow~{}\widetilde{x}=(n/{\tau})(x_{n-1}-x_{n})$ to rewrite the integral $\int_{\Omega}\mathrm{d}x_{n-1}$ as $({\tau}/{n})\int_{\widetilde{\Omega}}\mathrm{d}\widetilde{x}$.

[TABLE]

with $\dot{x}_{\tau}\in\widetilde{\Omega}=\mathbb{R}$. We note that the appearance of a velocity-like variable is unusual, as we deal with overdamped dynamics. Since $\dot{X}$ has no well-defined corresonding PDF for overdamped dynamics [53], we omit the notation used on the left side of (16i). With these steps, we finally obtain from (5.1), which is the FPE for the two-point PDF of the Markovian system (7), the second member of the FP hierarchy for the delayed system

[TABLE]

Comparison with the corresponding FPE (2.1) obtained from Novikov’s theorem reveals that Eq. (17) is in fact quite different. Our approach thus yields an alternative representation of the second member. However, we show in the Appendix B, that one can indeed transform (17) into (2.1). As we will argue below, this alternative representation is a new starting point to approximate the two-time PDF, $\rho_{2}(x,t;x_{\tau},t-\tau)$, for which, to our knowledge, no approaches have been reported. Moreover, also the Eqs. (5.1,5.1) might potentially serve for this purpose.

The main difference of (17) as compared to (2.1) (from Novikov’s theorem), is that it involves cross-correlations of the present and delayed system state with the delayed velocity, instead of the correlations among the system states at three different times [via the three-time PDF $\rho_{3}(x,t;x_{\tau},t-\tau,{x_{2\tau},t-2\tau})$]. This indicates that, instead of going back even further into the past with every member [by first taking into account $X(t-\tau)$, then $X(t-2\tau)$, then $X(t-3\tau)$ …], the here derived hierarchy rather collects more and more information about the dynamics at time $t-\tau$ [by taking into account $X(t-\tau)$, then $\dot{X}(t-\tau)$].

Both representations are valid and describe the same process. This shows that to (probabilistically) predict the future after a time $t$, either knowledge of the system states at all times $t,t-\tau,t-2\tau,...$ would be needed, or, complete knowledge about the preceding interval $[t-\tau,t]$ would suffice. On the stochastic level of description, the process which underlies the non-Markovian one, can either be considered as a process of the state vector $(X(t),X(t-\tau),X(t-2\tau),X(t-3\tau),..)$, or, equivalently (and actually more natural), of the state vector of all $X(s)$ $\in\mathbb{R}^{\infty}$, with $s\in[t-\tau,t)$. Moreover, the here presented results show that the “discretization of the past” (i. e., the discretization of $[t-\tau,t)$ into $n\in\mathbb{N}$ slices), as implied by the employed embedding technique, is a valid approach to time-continuous dynamics with delay.

A technical advantage of the representation (17) is that it does not involve unknown functional derivatives w. r. t. the noise, as opposed to (2.1), and only contains one diffusion term w. r. t. $x$ (and no additional diffusion term in $x_{\tau}$). In the next section, we will discuss an approximation scheme based on (17).

6 Approximation for the two-time PDF

In this section, we discuss the possibility of finding an approximate two-time PDF based on (17). We are particularly interested in an approach that is justified when the delay time $\tau$ is comparable with other timescales of interest, because in this regime, the non-Markovian effects are most prominent. We further focus on non-equilibrium steady states, which are typically approached after the transient behavior due to the initial condition has decayed. By making an ansatz for the higher order correlations, we find from (17) a closed equation that can be solved analytically.

As a first step, we assume that, if the delay time is not very small, the correlations between the system state $X(t)$ and the delayed displacement $\mathrm{d}X(t-\tau)$ can be neglected, as well as the correlations between system state and instantaneous displacement. This amounts to the ansatz $\langle\dot{X}(t-\tau)|X(t),X({t-\tau})\rangle\approx\langle\dot{X}(t-\tau)\rangle$. Inserting this into (17), one obtains

[TABLE]

involving the particle current $\langle\dot{X}\rangle$. In order to find an ansatz for the latter, one can proceed by utilizing the framework [54] (developed for Markovian systems), which gives a connection between $\langle\dot{X}\rangle$ and the probability current $J$ via $\langle\dot{X}\rangle=\int_{\Omega}J(x,t)\mathrm{d}x$, where $J$ is found from the FPE $\partial_{t}\rho_{1}=\partial_{x}J$ (written in the form of a continuity equation). For the delayed systems considered here $J(x,t)=\int_{\Omega}F(x,x_{\tau})\rho_{2}(x,x_{\tau})\mathrm{d}x_{\tau}-D_{0}\partial_{x}\rho_{1}(x)$. These assumptions give rise to the approximate, closed FPE (18) with $\langle\dot{X}\rangle=\iint_{\Omega}F(x,x_{\tau})\rho_{2}(x,x_{\tau})\mathrm{d}x_{\tau}\mathrm{d}x$. We expect the approximation to break down when $\tau$ is very small, i. e., when the correlations between $X(t)$ and $\mathrm{d}X(t-\tau)$ are not negligible.

For a common type of systems, this approximate equation becomes even simpler without requiring further assumptions or simplifications, that is, when $\langle\dot{X}\rangle=0$ in the steady-state. For example, this is always the case for systems that obey a symmetry w. r. t. coordinate inversions $x\to-x$. Then, (18) has the formal solution

[TABLE]

where $Z$ is a normalization constant. Further, $h(x_{\tau})$ is a function that is not determined by (18) alone, but is, however, determined by the steady-state constraint $\partial_{t}\rho_{1}=0$, which implies $\rho_{1}(x)=\rho_{1}(x_{\tau})\Rightarrow\int_{\Omega}\rho_{2}(x,x_{\tau})\mathrm{d}x=\int_{\Omega}\rho_{2}(x,x_{\tau})\mathrm{d}x_{\tau}$.

As an important example, we consider systems with *linear *delay force, and (inversion-symmetric) nonlinear, nondelayed force $F_{\mathrm{nd}}$, i. e., $F(x,x_{\tau})=F_{\mathrm{nd}}(x)+kx_{\tau}$. Then $h=e^{\int F_{\mathrm{nd}}(x_{\tau})\mathrm{d}x_{\tau}}$ yields

[TABLE]

which fulfills (18) and the steady-state constraint.

A peculiar feature of (20) is that it does not explicitly depend on the precise value of $\tau$ (unless $F_{\mathrm{nd}}$ does). However, it has been observed in earlier studies for both, nonlinear [31] and linear systems [29] (where the PDFs can be computed analytically), that the moments and correlation functions show a saturation in $\tau$. In other words, $\rho_{1}$ and $\rho_{2}$ become in fact independent of $\tau$ in the regime of large delay times. Since we expect our approximation to be justified only in the large $\tau$ regime, we suspect that we can understand the independence of (20) on $\tau$ in this sense. We recall that the approximation is expected to break down for small values of $\tau$, because then the correlations between $\dot{X}(t-\tau)$ and $X(t)$ are non-negligible.

6.1 Application to a bistable delayed system

To test the approximation (20), we consider a concrete example, namely a bistable system [41, 37, 31] given by the rescaled LE444The dimensionless LE is obtained by rescaling position $x/\sigma\to x$ and time $(D_{0}/\sigma^{2})t\to t$, $(D_{0}/\sigma^{2})\tau\to\tau$, with the the minima $x\pm\sigma$ of the doublewell potential $V=V_{0}[(x/\sigma)^{4}-2(x/\sigma)^{2}]$, and by rescaling the parameters $V_{0}/(\gamma D_{0})\to V_{0}$ and $k/(\gamma D_{0})\to k$ with the bath’s thermal energy $\gamma D_{0}$.

[TABLE]

Equation (21) describes the dynamics of a stochastic system in a static doublewell potential $V=V_{0}(x^{4}-2x^{2})$ with minima at $x\pm 1$, subject to a linear delay force stemming from an “optical tweezers” potential $(k/2)(x-x_{\tau})^{2}$ [55]. Using Eq. (20), we can readily write down the approximate steady-state two-time PDF

[TABLE]

It is well-known that the bistable system (21) exhibits delay-induced oscillations between the two potential minima, when the delay time is larger than the the typical time between noise-induced jumps, and when $|k|$ has the same order of magnitude as $V_{0}$. Because these jumps manifest in the two-time PDF, as will be explained below, we consider this parameter regime (we fix $V_{0}\!=\!3.5$ and $\tau\!=\!10$) and only vary $k$, particularly focusing on the question whether the approximation can capture the delay-induced jumps. To test the results from a quantitative point of view, we have additionally performed direct Brownian dynamics simulations of Eq. (21) with $>2000$ realizations, temporal discretization $10^{-4}$, and cutoff $>10^{6}$ time steps to reach steady-state conditions.

First, we consider the system in the absence of the delay force ($k=0$), where the dynamics is dominated by the noise, and the probability to find $X$ is always highest around the potential minima $x=\pm 1$. As thermal jumps are very likely on the timescale of $\tau$, the chance of finding the system in one well is independent of the question whether it was on the same or on the other side $\tau$ ago, meaning that the probabilities to find $X(t)$ and $X(t-\tau)$ are independent from each other. We hence expect the joint two-time PDF555Due to the symmetry of the external potentials, we further always expect that $\rho_{2}(x,x_{\tau})$ is symmetric w. r. t.  $\{x,x_{\tau}\}\to-\{x,x_{\tau}\}$. $\rho_{2}(x,x_{\tau})$ to have pronounced peaks of identical heights around $\{x,x_{\tau}\}\approx\{\pm 1,\pm 1\}$ and $\{\pm 1,\mp 1\}$. This expectation is confirmed by numerical results, and is correctly described by our approximation (22), as plotted in Fig. 3 (left plots, upper and lower panel, respectively). (We note that for $k=0$ the dynamics is Markovian and the system asymptotically approaches thermal equilibrium.) We now consider the impact of $k\neq 0$ (where a non-equilibrium steady-state is approached).

For negative $k<0$, i. e, when the linear feedback force is directed towards the delayed system state, the induced oscillations are known to have a mean period of roughly $2\tau$ [31, 37, 41]. This means that the jumps reduce the joint probability to find the system in the same well at times $t-\tau$ and $t$. We hence expect the delay-induced oscillations to lower the two-time PDF around $\{\pm 1,\pm 1\}$, and increase it around $\{\pm 1,\mp 1\}$. This effect is indeed seen in both, the numerically obtained PDF and the approximate results from (22), see Fig. 3. Interestingly, by changing the sign of $k$, the mean period of the oscillations can be switched to about $\tau$ [41]. Thus, we expect the delay-induced oscillations to have, in fact, the opposite effect for positive $k$, i. e., enhance the PDF at $\{\pm 1,\pm 1\}$, and reduce it around $\{\pm 1,\mp 1\}$. Again, this expected behavior is confirmed by simulation results, and as well captured nicely by the proposed approximation, see Fig. 4.

The comparison between numerical results and the approximate $\rho_{2}$ from (22) in Figs. 3, 4 includes a broad range of $k$ values. Remarkably, we find that the proposed approximation yields very good results in all cases considered. In particular, it captures the main characteristics of the two-time PDF described above, for negative and for positive $k$. Even quantitatively, the agreement is very reasonable. We further considered higher barrier heights ($V_{0}=5$), smaller delays ($\tau=5$), as well as higher feedback strength (up to $|k|=20$) and found again very good agreement (not shown here). On the other hand, the approximation breaks down if the delay time is very short (we have tested, e. g., $\tau=0.1$). We suspect that in this regime, the correlations between $\dot{X}(t-\tau)$ and $X(t)$ are not negligible (hence, the assumption $\langle\dot{X}(t-\tau)|X(t),X({t-\tau})\rangle\approx\langle\dot{X}(t-\tau)\rangle$ is not justified).

The accuracy of the approximate $\rho_{2}$ is indeed a remarkable result, in particular, because the underlying assumption appears to be rather crude. This demonstrates that an approximation scheme on the level of the second member of the FP hierarchy is somewhat superior as compared to an approach on the level of the first member. To further support this argument, we compare the prediction of approximation (22) to earlier approaches for the one-time PDF known from the literature (to the best of our knowledge, approaches for the two-time PDF have not been presented elsewhere). To this end, we marginalize $\rho_{2}$ from Eq. (22), i. e., $\rho_{1}(x)=\int_{\Omega}\rho_{2}(x,x_{\tau})\mathrm{d}x_{\tau}$ and compare the result with the small $\tau$ expansion from [35], and the force-linearization closure (FLC) [31]. As can be seen in Fig. 5, the novel approximation yields the best results in the regime of large delay times. Admittedly, these other approaches are, in fact, known to be not justified for large $\tau$ values, or when the delay-induced jumps are very likely, respectively [31]. We nevertheless use them to compare with, due to the lack of other approximation schemes for this nonlinear system. In Appendix D, a short delay times is considered, where the other approximations are more reasonable.

We note that besides the features discussed above, the numerically obtained two-time PDFs in Figs. 3 and 4 show very interesting further details. An example is the broken symmetry $\rho_{2}(a,t;a^{\prime},t-\tau)\neq\rho_{2}(a^{\prime},t,a,t-\tau)$, which is not captured by the approximation (22). A study of these features is beyond the scope of this work, but it would be interesting to see how these details are connected with dynamical and thermodynamical properties of the non-equilibrium steady-state.

7 Conclusions

In this work, we have discussed the probabilistic description of delayed stochastic systems. As a starting point, we have reviewed an earlier approach based on Novikov’s theorem, from which a Fokker-Planck description can be derived in the form of an infinite hierarchy [30, 31, 32]. The first member is the well-known FPE for the one-time PDF, which contains the two-time PDF and is thus not closed. Still, this equation has in the past proven itself to be an important tool in the search for exact results [36, 37] and a valuable base for approximations [31, 38]. The main purpose of our work was to shed light onto the higher members of this hierarchy, which have rarely been discussed in earlier literature.

To this end, we introduced an alternative FP description basing on a Markovian embedding technique, where the delayed process is represented as the dynamics of one variable of an $(n+1)$-variable Markovian system. By projecting the corresponding closed, Markovian FPE onto lower-dimensional subspaces, and taking the limit $n\to\infty$ in the projected equations, we derived a hierarchy of FPEs. While the first member is identical, the representation of the higher members derived in this way differs from the one obtained by Novikov’s theorem. We demonstrated that both hierarchies can, however, be converted into each other. This shows that a “discretization of the past” (i. e., the discretization of $[t-\tau,t)$ into $n\in\mathbb{N}$ equidistant slices), which is implied by the employed embedding technique, is a valid approach to time-continuous dynamics with delay. Furthermore, the higher members derived here [Eqs. (5.1,5.1) or (17)], might serve as a new starting point for approximations.

For the second member, which we focused on in this paper, we have proposed one possible approximation scheme, and applied it to a bistable model with linear delay force. We found remarkably good agreement with simulation results, despite the simple form of the obtained approximate two-time PDF. With this, we present for the first time approximate results for this quantity which is indeed central in the probabilistic description of stochastic systems with time delay. The proposed scheme renders very good results in the regime of large delay times, for both cases, feedback forces that are directed towards, or away from the delayed position. In particular, it captures the main features of the two-time probability density which has very distinct characteristics depending on the direction of the feedback, related to delay-induced oscillations. Moreover, the corresponding approximate one-time PDF (which can be obtained by marginalization) outperforms earlier approaches known from the literature.

The benefit of an approximation schemes based on Eqs. (17) or (2.1), as compared to an approach operating on the first member (like the ones known from the literature), is twofold. Besides expectantly yielding an improved approximation of the one-time PDF, an approximate two-time PDF itself could be used directly to estimate various quantities of interest. An example is the distribution of the fluctuating heat in nonlinear, delayed systems [40, 37]. In future studies, we want to explore these possibilities with the new approximation discussed here.

The development of further approximation schemes for the PDFs is beyond the scope of the present paper. However, besides the approach discussed in Sec. 6, another, quite different type of approximation is in fact readily given by the Markovian embedding approach, namely, by directly considering the extended Markovian system with a *finite *number of auxiliary variables $n$ [49]. Then, Eq. (11) is a well-defined, closed FPE, which can be used to calculate the $n$-point PDF by means of standard techniques for Markovian FPEs (at least numerically). By marginalization of the resulting PDFs, one would recover the corresponding (approximate) PDFs of the delayed process. On the level of the stochastic delay differential equation (the LE), this is in fact a common method [56, 57, 58, 6]. It amounts to approximating the discrete delay by a Gamma-distributed memory kernel (8) with finite width. However, a systematic investigation of this approximation (including an error estimate), has, to the best of our knowledge, not been performed yet. A similar finite-dimensional approximation for linear systems is discussed in [59]. A more detailed study of these approximations is subject of ongoing work.

Finally, we would like to reconsider the Markovian embedding from a conceptual point of view. Non-Markovianity, in general, and specifically the discrete time delay considered here, is a consequence of an incomplete description of a system, i. e., missing degrees of freedom. The delayed force might either stem from a specifically designed external agent, a feedback controller, or, more generally, from any part of a more complex system (the “super-system”), which is not explicitly modeled within the considered description, but interacts with the system at hand in form of delayed forces. The Markovian auxiliary system, in which we here embed the delayed process, allows us to mimic the external agent (or other parts of the “super-system”) without explicitly modeling it in detail. In this sense, the auxiliary variables can be interpreted as a “substitute” of the hidden degrees of freedom in the system. On the other hand, it is well-known that hidden degrees of freedom have a significant impact on the entropy production and other thermodynamical notions, see for example [60, 61]. Therefore, a study of the consequences of our treatment of delayed systems from a stochastic thermodynamical (and information-theoretical) perspective would be very interesting. This will be a focus of future work.

Acknowledgements

This work was funded by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) - Projektnummer 163436311 - SFB 910. We further thank Christian Kuehn for fruitful discussions.

Appendix A APPENDIX: The discrete delayed case as limit of the coarse-grained dynamics

In this appendix, we show how the dynamical equation of the delayed system (1) can be recovered by taking the limit $n\to\infty$ of the corresponding equation for the embedded system (7). To this end, we project onto the the dynamics of the variable $X_{0}$. Until completion of the projection, we will keep $n$ finite. Formally integrating Eq. (7b) (by first solving the homogeneous problem and then finding a specific solution by variation of constants) yields

[TABLE]

In principle, this solution can be simplified using an appropriate initial condition for the auxiliary variables, for example $X_{j}(0)=0,~{}\forall j>0$. However, we here proceed with the general case, i. e., without specifying the initial condition. As a next step (assuming $n>1$), we iteratively plug Eq. (23) into the corresponding solution for $j+1$ [given by Eq. (23) with $j+1$ instead of $j$]. For $1\leq j<n$, we find

[TABLE]

with the Gamma-distributed memory kernel defined in Eq. (8). The first term in (24d) vanishes trivially, if the auxiliary variables satisfy the initial condition $X_{j\in\{1,2,..,n\}}(-\tau)=0$. Interestingly, this terms also vanishes if the asymptotic behavior at $t\to\infty$ is considered, due to the exponential damping with time. Plugging (24d) into (7a) [neglecting the initial condition term] then yields the projected Eq. (8).

Now we consider the limit $n\to\infty$ of the memory kernel (8). First, we find that the mean value of the kernel and the integral over it are, independently of the value of $n$, given by $\mu=\int_{0}^{\infty}K_{n}(T)T\mathrm{d}T=\tau$ and $\int_{0}^{\infty}K_{n}(T)\mathrm{d}T=1$, respectively. Furthermore, for large $n$, the value of the kernel at $\tau$ can be estimated by using the Stirling formula $n!\approx\sqrt{2\pi n}(n/e)^{n}$, yielding $K_{n}(\tau)=\frac{{n}^{n}\,e^{-n}}{(n-1)!\tau}\approx\frac{n!}{(n-1)!\sqrt{2\pi n}\tau}=\sqrt{\frac{n}{2\pi}}\frac{1}{\tau}\to\infty.$ On the other hand, the variance of the kernel vanishes for $n\to\infty$, since $\int_{0}^{\infty}K_{n}(T)T^{2}\mathrm{d}T-\mu^{2}=\tau^{2}/{n}\to 0$. This implies

[TABLE]

Thus, the projected equation (8) is indeed equal to Eq. (1) in the limit $n\to\infty$. We note that even the transient dynamics (i. e., finite $t$) is only equivalent for the initial conditions $X_{j\in\{1,2,..,n\}}(-\tau)\equiv 0$, as mentioned above.

Appendix B APPENDIX: Connection to Fokker-Planck hierarchy from Novikov’s theorem

While the FPE for $\rho_{1}$ obtained by our approach is equivalent to the first member of the FP hierarchy from Novikov’s theorem (2), the second members differ. However, the apparent disagreement can be resolved. Indeed, as we show below, our result (17) can be transformed to the representation obtained from Novikov’s theorem [given in (2.1)]. We focus on the term which differs, i. e., the drift term $\partial_{x_{\tau}}[\rho_{2}(x,t;x_{\tau},t-\tau)$ $\langle\dot{X}(t-\tau)|x=X(t),x_{\tau}=X({t-\tau})\rangle]$. To this end, we first use the definition of $\rho_{2}$ via delta-distributions, plug in the LE (1), which yields

[TABLE]

In combination with (26a,26), we obtain the identity

[TABLE]

proving that the second member of the here presented FP hierarchy, Eq. (17), is identical to the corresponding one (2.1) obtained from Novikov’s theorem.

Appendix C APPENDIX: Novikov’s theorem

In this section, we establish the relation [34]

[TABLE]

for a functional $\Lambda$ of a Gaussian white noise $\xi$. It links the functional derivative w. r. t. the noise to the cross-correlation between functional and noise. First, we consider the case where the ensemble is w. r. t. a *fixed *initial condition $\phi=X(-\tau\leq t\leq t)$, yielding $\Lambda(0)={\Lambda_{0}}$. In a second step we will generalize towards initial conditions drawn from an arbitrary distribution.

[TABLE]

This result can readily be generalized to ensembles where the initial conditions are instead drawn from an arbitrary, normalized distribution $P(\Lambda_{0})$. In particular, Eq. (29f) implies

[TABLE]

This is the relation (28), which is often referred to as Novikov’s theorem.

Appendix D APPENDIX: The approximation in the regime of short delay times

The approximate two-time PDF (22) is, in general, not reliable, if the delay time is very short, as the underlying assumption is violated. For the moderate $|k|/V_{0}$ and small $\tau$ values considered in Fig. 6 (a), the corresponding one-time PDF (obtained by marginalization) still gives a reasonable approximation. Here also the other approaches from the literature are justified. For the large (and negative) $|k|/V_{0}$ value in Fig. 6 (b), only the small $\tau$ expansion yields a reasonable approximation.

• [1]

D.M. Kane, K.A. Shore (eds.), Unlocking dynamical diversity: optical feedback effects on semiconductor lasers (John Wiley & Sons, 2005)

• [2]

Y. Liu, H. Chen, J. Liu, P. Davis, T. Aida, Phys. Rev. A 63(3), 031802 (2001)

• [3]

A. Waibel, T. Hanazawa, G. Hinton, K. Shikano, K.J. Lang, Backpropagation: Theory, Architectures and Applications (Lawrence Erlbaum Assoc. Mahwah, NJ, 1995)

• [4]

A. Longtin, J.G. Milton, J.E. Bos, M.C. Mackey, Phys. Rev. A 41(12), 6992 (1990)

• [5]

J.R. Cabral, H. Luckhoo, M. Woolrich, M. Joensson, H. Mohseni, A. Baker, M.L. Kringelbach, G. Deco, Neuroimage 90, 423 (2014)

• [6]

V. Gupta, K. Kadambari, et al., Biological cybernetics 104(6), 369 (2011)

• [7]

A. Tambue, E.K. Brown, S. Mohammed, Commun. Nonlinear Sci. Numer. Simul. 20(1), 281 (2015)

• [8]

J.L. Callen, M. Khan, H. Lu, Contemp. Account. Res. 30(1), 269 (2013)

• [9]

E. Schöll, S.H.L. Klapp, P. Hövel (eds.), Control of self-organizing nonlinear systems (Springer, 2016)

• [10]

A. Carmele, J. Kabuss, F. Schulze, S. Reitzenstein, A. Knorr, Phys. Rev. Lett. 110(1), 013601 (2013)

• [11]

N. Német, S. Parkins, Phys. Rev. A 94(2), 023809 (2016)

• [12]

M. Durve, A. Saha, A. Sayeed, Eur. Phys. J. E 41(4), 49 (2018)

• [13]

U. Khadka, V. Holubec, H. Yang, F. Cichos, Nat. Commun. 9, 3864 (2018)

• [14]

S.A.M. Loos, R. Gernert, S.H.L. Klapp, Phys. Rev. E 89(5), 052136 (2014)

• [15]

M. Mijalkov, A. McDaniel, J. Wehr, G. Volpe, Phys. Rev. X 6(1), 011008 (2016)

• [16]

N. Bruot, L. Damet, J. Kotar, P. Cicuta, M.C. Lagomarsino, Phys. Rev. Lett. 107(9), 094101 (2011)

• [17]

C. Masoller, Phys. Rev. Lett. 88, 034102 (2002)

• [18]

R. Gernert, S.A.M. Loos, K. Lichtner, S.H.L. Klapp, in Control of Self-Organizing Nonlinear Systems (Springer, 2016), pp. 375–392

• [19]

I. Schneider, Philos. Trans. Royal Soc. A 371(1999), 20120472 (2013)

• [20]

A. Zakharova, S.A.M. Loos, J. Siebert, A. Gjurchinovski, J.C. Claussen, E. Schöll, in Control of self-organizing nonlinear systems (Springer, 2016), pp. 3–23

• [21]

I. Schneider, M. Bosewitz, Disc. Cont. Dyn. Syst. A 36(1), 451 (2016)

• [22]

E. Schöll, H.G. Schuster (eds.), Handbook of chaos control (John Wiley & Sons, 2008)

• [23]

R. Zwanzig, J. Stat. Phys. 9(3), 215 (1973)

• [24]

R. Kubo, Reports on progress in physics 29(1), 255 (1966)

• [25]

M. Krüeger, C. Maes, J. Phys.: Condens. Matter 29(6), 064004 (2016)

• [26]

J.D. Bao, P. Hänggi, Y.Z. Zhuo, Phys. Rev. E 72(6), 061107 (2005)

• [27]

C. Maes, S.R. Thomas, Phys. Rev. E 87(2), 022145 (2013)

• [28]

C. Maes, J. Stat. Phys. 154(3), 705 (2014)

• [29]

U. Küchler, B. Mensch, Stoch. Stoch. Rep. 40(1-2), 23 (1992)

• [30]

T.D. Frank, Phys. Rev. E 72(1), 011112 (2005)

• [31]

S.A.M. Loos, S.H.L. Klapp, Phys. Rev. E 96(13), 012106 (2017)

• [32]

M.L. Rosinberg, T. Munakata, G. Tarjus, Phys. Rev. E 91, 042114 (2015)

• [33]

L. Giuggioli, T.J. McKetterick, V.M. Kenkre, M. Chase, J. Phys. A 49(38), 384002 (2016)

• [34]

E.A. Novikov, Sov. Phys. JETP 20(5), 1290 (1965)

• [35]

S. Guillouzic, I. L’Heureux, A. Longtin, Phys. Rev. E 59(4), 3970 (1999)

• [36]

T.D. Frank, P.J. Beek, Phys. Rev. E 64(2), 021917 (2001)

• [37]

S.A.M. Loos, S.H.L. Klapp, Sci. Rep. 9, 2491 (2019)

• [38]

T.D. Frank, Phys. Rev. E 71(3), 031106 (2005)

• [39]

U. Seifert, Rep. Prog. Phys. 75(12), 126001 (2012)

• [40]

K. Sekimoto, Stochastic energetics, vol. 799 (Springer, 2010)

• [41]

L.S. Tsimring, A. Pikovsky, Phys. Rev. Lett. 87(25), 250602 (2001)

• [42]

A. Puglisi, D. Villamaina, EPL 88(3), 30004 (2009)

• [43]

A. Crisanti, A. Puglisi, D. Villamaina, Phys. Rev. E 85(6), 061127 (2012)

• [44]

P. Siegle, I. Goychuk, P. Talkner, P. Hänggi, Phys. Rev. E 81(1), 011136 (2010)

• [45]

P. Siegle, I. Goychuk, P. Hänggi, EPL 93(2), 20002 (2011)

• [46]

D. Villamaina, A. Baldassarri, A. Puglisi, A. Vulpiani, J. Stat. Mech. Theory Exp 2009(07), P07024 (2009)

• [47]

A. Longtin, Complex time-delay systems: theory and applications (Springer-Verlag Berlin Heidelberg, 2010)

• [48]

H.L. Smith, An introduction to delay differential equations with applications to the life sciences, vol. 57 (Springer New York, 2011)

• [49]

S.I. Niculescu, K. Gu, Advances in time-delay systems, vol. 38 (Springer Science & Business Media, 2012)

• [50]

T.D. Frank, Phys. Rev. E 69, 061104 (2004)

• [51]

Y. Zheng, X. Sun, Discrete Continuous Dyn. Syst. Ser. B 22(9) (2017)

• [52]

T.D. Frank, P.J. Beek, R. Friedrich, Phys. Rev. E 68(2), 021912 (2003)

• [53]

C.W. Gardiner, Handbook of Stochastic Methods, 2nd edn. (Springer, 2002)

• [54]

P. Reimann, Phys. Rep. 361 (2002)

• [55]

J. Kotar, M. Leoni, B. Bassetti, M.C. Lagomarsino, P. Cicuta, Proc. Natl. Acad. Sci. U.S.A. 107(17), 7669 (2010)

• [56]

M. Rosinberg, G. Tarjus, T. Munakata, Phys. Rev. E 98(3), 032130 (2018)

• [57]

B. De Vries, J.C. Príncipe, in Advances in neural information processing systems (1991), pp. 162–168

• [58]

K. Rateitschak, O. Wolkenhauer, Mathematical biosciences 205(2), 163 (2007)

• [59]

A. René, A. Longtin, Chaos 27(11), 114322 (2017)

• [60]

J. Mehl, B. Lander, C. Bechinger, V. Blickle, U. Seifert, Phys. Rev. Lett. 108(22), 220601 (2012)

• [61]

K. Kawaguchi, Y. Nakayama, Phys. Rev. E 88(2), 022147 (2013)

• [62]

M.J. Hall, M. Reginatto, Fundam. Theor. Phys. 184 (2016)