Brownian molecules formed by delayed harmonic interactions

1.1. Feedback systems

1.2. Stochastic delay differential equations

1.3. Outline

We consider a two-dimensional system of N overdamped Brownian particles coupled via time-delayed harmonic pair interactions given by the potential

$$\begin{eqnarray}&&{V}_{2}[{r}_{{ij}}(t-\tau )]=\displaystyle \frac{k}{2}{\left[{r}_{{ij}}(t-\tau )-R\right]}^{2},\end{eqnarray} \tag{ 1 }$$

depicted in figure 1 by springs connecting the individual particles. In equation (1), R > 0 denotes the equilibrium spring length, k their stiffness, and ${r}_{{ij}}(t-\tau )=\left|{{\bf{r}}}_{i}(t-\tau )-{{\bf{r}}}_{j}(t-\tau )\right|$ is the distance between the particles i and j located at positions r_i and r_j at an earlier time t − τ. Clearly, the picture of linear springs can properly represent the time-delayed interactions only for a vanishing time delay τ.

Zoom In Zoom Out Reset image size

Altogether, the particles diffuse in the composity potential

$$\begin{eqnarray}&&V=\displaystyle \frac{1}{2}\displaystyle \sum _{(i,j)}{V}_{2}({r}_{{ij}}),\end{eqnarray} \tag{ 2 }$$

where the summation runs over all pairs (i, j), so that the ith particle is driven by the force ${{\bf{F}}}_{i}=-{{\rm{\nabla }}}_{i}V\,=({\partial }_{{x}_{i}}V,{\partial }_{{y}_{i}}V)$. Because at time t the particle feels the value of the potential corresponding to its position at time t − τ, here x_i and y_i denote the Cartesian coordinates of the position vector r_i in the past. In effect, the N-particle system therefore obeys the set of nonlinear delayed Langevin equations

$$\begin{eqnarray}&&{\dot{{\bf{r}}}}_{i}(t)=-\displaystyle \frac{k}{\gamma }\displaystyle \sum _{j\ne i}\left[{r}_{{ij}}(t-\tau )-R\right]{{\bf{e}}}_{{ij}}(t-\tau )+\sqrt{2{D}_{0}}{{\boldsymbol{\eta }}}_{i}(t),\ \ i=1,\,\ldots ,\,N.\end{eqnarray} \tag{ 3 }$$

The unit vector e_ij = r_ij/r_ij points from particle j to particle i and the diffusion coefficient D₀ = (βγ)⁻¹ is related to the inverse temperature β = 1/k_BT and the friction coefficient γ via the Einstein relation (k_B denotes the Boltzmann constant). The vectors ${{\boldsymbol{\eta }}}_{i}$ comprise independent Gaussian white noises satisfying the relations

$$\begin{eqnarray}&&\langle {{\boldsymbol{\eta }}}_{i}(t)\rangle =0,\qquad \left\langle {\eta }_{i}^{\alpha }({t}_{1}){\eta }_{j}^{\beta }({t}_{2})\right\rangle ={\delta }_{{ij}}{\delta }_{\alpha \beta }{\delta }_{{\rm{D}}}({t}_{1}-{t}_{2}).\end{eqnarray} \tag{ 4 }$$

The numbers α and β label the components of the vector ${{\boldsymbol{\eta }}}_{i}$, while i refers to the specific particle. Note that the noise and the friction in equation (3) are related by the fluctuation-dissipation theorem [59] for a vanishing delay τ = 0 only, and that the system is always out of thermodynamic equilibrium for τ > 0 [56]. In order to obtain a model that would obey the fluctuation-dissipation theorem, one should consider a different noise correlation function (4).

For small enough values of the time delay, the particles form, after an initial transient period, highly symmetric molecular-like structures, some of which are displayed in figure 4(a) in section 3. For N = 2 (dimer) and N = 3 (trimer) the steady-state structures occupy the global minimum of the potential V. For N > 3 the global minimum becomes inaccessible due the chosen two-dimensional geometry and the resulting structures are thus frustrated in the sense that some of the springs do not reach their equilibrium length in the steady-state. The structures are dynamical, due to the Brownian motion of the particles, which persistently kicks the system out of the minimum of the potential energy (2). The effect of the delay is that the system may exhibit exponentially decaying oscillations on its return to the energy minimum. The decay rate of these oscillations decreases with increasing delay, and, for delays larger than a certain threshold, their amplitude exponentially increases. This is because large delays induce in the system a 'swing effect', when the repulsive force from one side of the potential propels the particle to a 'higher' position at its opposite side, and so on.

Within the equilibrium model that obeys the fluctuation-dissipation theorem, the stationary probability density function (PDF) for positions of the individual particles would simply be the Boltzmann distribution $\exp (-\beta V)/Z$ with potential V, inverse temperature β, and normalization Z. However, the physical situation at hand, where the delay is interpreted as a result of a feedback control mechanism and thus is independent of the noise, requires the more involved description with equation (4) that leads to non-trivial non-equilibrium steady states. Consequently, the Boltzmann distribution can no longer be assumed. For the simplest case of a dimer with short delay time, we will find an approximately Gaussian distribution, corresponding to a ('deformed') Boltzmann factor at an effective temperature. For larger particle numbers and longer delay times, the situation becomes more complicated.

A similar system with a quasi-constant force between the particles (constant upto a change of sign at distance R), i.e. obeying the set of Langevion equations

$$\begin{eqnarray}&&{\dot{{\bf{r}}}}_{i}(t)=-\displaystyle \frac{k}{\gamma }\displaystyle \frac{{\displaystyle \sum }_{j\ne i}\mathrm{sign}\left[{r}_{{ij}}(t-\tau )-R\right]{{\bf{e}}}_{{ij}}(t-\tau )}{\left|{\displaystyle \sum }_{j\ne i}\mathrm{sign}\left[{r}_{{ij}}(t-\tau )-R\right]{{\bf{e}}}_{{ij}}(t-\tau )\right|}+\sqrt{2{D}_{0}}{{\boldsymbol{\eta }}}_{i}(t),\end{eqnarray} \tag{ 5 }$$

i = 1, ..., N, with sign(x) denoting the signum function, was discussed earlier in [11]. The main difference from our setting (3) is that, in equation (5), the absolute value of the force does not depend on the interparticle distances and the particle number N. The main benefit of assuming the harmonic potential in equation (3) is that it allows much more complete analytical treatment. To allow for an easy comparison of the two models, we illustrate our results using parameters inspired by [11]. In the following, we first review some analytical results for stochastic dynamics of dimers and trimers. On this basis, we will return to the discussion of the emerging structures in section 3.

2.1. Center of mass

Similarly as for the dynamics considered in [11], the center of mass coordinate ${{\bf{r}}}_{{\rm{c}}}\equiv \left({\sum }_{i=1}^{n}{{\bf{r}}}_{i}\right)/N$ of the system obeys the Langevin equation

$$\begin{eqnarray}&&{\dot{{\bf{r}}}}_{{\rm{c}}}(t)=\sqrt{2{D}_{c}}{{\boldsymbol{\eta }}}_{{\rm{c}}}(t),\end{eqnarray} \tag{ 6 }$$

where ${{\boldsymbol{\eta }}}_{{\rm{c}}}\equiv {\sum }_{i=1}^{N}{{\boldsymbol{\eta }}}_{i}/\sqrt{2}$ denotes Gaussian white noise satisfying equation (4) (with the labels i, j replaced by c) and the diffusion coefficient D_c = D/N. Regardless of the interactions, the center of mass performs ordinary Brownian motion and, assuming the center of mass is at time t = 0 located at the point r₀, the PDF for r_c reads

$$\begin{eqnarray}&&{P}_{N}^{{\rm{c}}}({\bf{r}},t)=\sqrt{\displaystyle \frac{N}{4\pi {Dt}}}\exp \left[-N\displaystyle \frac{{\left({\bf{r}}-{{\bf{r}}}_{0}\right)}^{2}}{4{Dt}}\right].\end{eqnarray} \tag{ 7 }$$

2.2. Dimer

Let us now consider the simplest case of two interacting particles. For N = 2, equation (3) yields the system of four coupled equations of motion:

$$\begin{eqnarray}&&{\dot{{\bf{r}}}}_{1}(t)=-\displaystyle \frac{k}{\gamma }\left[r(t-\tau )-R\right]{\bf{e}}(t-\tau )+\sqrt{2D}{{\boldsymbol{\eta }}}_{1}(t),\end{eqnarray} \tag{ 8 }$$

$$\begin{eqnarray}&&{\dot{{\bf{r}}}}_{2}(t)=+\displaystyle \frac{k}{\gamma }\left[r(t-\tau )-R\right]{\bf{e}}(t-\tau )+\sqrt{2D}{{\boldsymbol{\eta }}}_{2}(t),\end{eqnarray} \tag{ 9 }$$

where we have used the abbreviations e(t) ≡ e₁₂(t) and r(t) ≡ r₁₂(t). In the previous section, we have already resolved the dynamics of the center of mass coordinate for arbitrary N. Now, we consider only the dynamics of the relative coordinate r = r₁₂ = r₁ − r₂ which obeys the equation of motion

$$\begin{eqnarray}&&\dot{{\bf{r}}}(t)=-\omega \left[r(t-\tau )-R\right]{\bf{e}}(t-\tau )+\sqrt{4D}{{\boldsymbol{\eta }}}_{{\rm{r}}}(t)\end{eqnarray} \tag{ 10 }$$

with frequency

$$\begin{eqnarray}&&\omega \equiv 2k/\gamma \end{eqnarray} \tag{ 11 }$$

and ${{\boldsymbol{\eta }}}_{{\rm{r}}}\equiv \left({{\boldsymbol{\eta }}}_{1}-{{\boldsymbol{\eta }}}_{2}\right)/\sqrt{2}$ describing Gaussian white noise satisfying equation (4) with the vector components i, j replaced by the label r.

Projecting equation (10) on the direction of the bond at time t (by multiplication with ${\bf{e}}(t)=\left(\cos \varphi (t),\sin \varphi (t)\right))$ and on the direction perpendicular to the bond (by multiplying it with ${{\bf{e}}}_{\varphi }=\left(-\sin \varphi (t),\cos \varphi (t)\right)$), we obtain the equations

$$\begin{eqnarray}&&\dot{r}(t)=-\omega \left[r(t-\tau )-R\right]\cos \left[\varphi (t,t-\tau )\right]+\sqrt{4D}{\eta }_{{\rm{r}}}^{r}(t),\end{eqnarray} \tag{ 12 }$$

$$\begin{eqnarray}&&\dot{\varphi }(t)=\omega \displaystyle \frac{r(t-\tau )-R}{r(t)}\sin \left[\varphi (t,t-\tau )\right]+\sqrt{\displaystyle \frac{4D}{{r}^{2}(t)}}{\eta }_{{\rm{r}}}^{\varphi }(t),\end{eqnarray} \tag{ 13 }$$

where φ(t, t − τ) = φ(t) − φ(t − τ) denotes the change of orientation of the vector e(t − τ) during time τ. Above, we used the formulas r ≡ re, $\dot{{\bf{r}}}\equiv \dot{r}{\bf{e}}+r\dot{\varphi }{{\bf{e}}}_{\varphi }$ and ${{\boldsymbol{\eta }}}_{{\rm{r}}}\equiv {\eta }_{{\rm{r}}}^{r}{\bf{e}}+{\eta }_{{\rm{r}}}^{\varphi }{{\bf{e}}}_{\varphi }$.

From symmetry considerations, it follows that the stationary PDF for the orientation must be constant in φ. To gain analytical insight into the dynamics and PDF of the bond-length r, we linearize the coupled Langevin equations (12) and (13). If the angle dependent stiffness $2k/\gamma \cos \left[\varphi (t,t-\tau )\right]$ in equation (12) is strong enough such that the terms proportional to [r(t − τ) − R]/R can safely be neglected independently of the noise strength, the formula (13) for the angle assumes the form⁵

$$\begin{eqnarray}&&\dot{\varphi }(t)=\sqrt{\displaystyle \frac{4D}{{R}^{2}}}{\eta }_{{\rm{r}}}^{\varphi }(t).\end{eqnarray} \tag{ 14 }$$

The corresponding transition PDF (Green's function) for change of the orientation by φ(t, t − τ) = φ(t) − φ(t − τ) during the time interval τ reads [60–62] $p\left[\varphi (t,t-\tau ),\tau \right]={\left(2\pi \right)}^{-1}+{\pi }^{-1}{\sum }_{m=1}^{\infty }\cos \left[m\varphi (t,t-\tau )\right]\exp \left[-2{m}^{2}\tau D/{R}^{2}\right]$. Using this function, we average equation (12) over φ(t, t − τ) obtaining

$$\begin{eqnarray}&&\dot{r}(t)=-{\omega }_{\tau }\left[r(t-\tau )-R\right]+\sqrt{4D}{\eta }_{{\rm{r}}}^{r}(t),\end{eqnarray} \tag{ 15 }$$

where ${\omega }_{\tau }=\omega \exp (-2D\tau /{R}^{2})$ is the natural relaxation rate. Note that the same formula with ω_τ substituted by ω is obtained by simply assuming that the change of the orientation φ(t, t − τ) of the bond per delay time τ is small, i.e. for $2\tau D/{R}^{2}\ll 1$. The main difference between the two approximations is that D/R² does not necessarily have to be small in the first case. We will discuss the regime of validity of the equation (15) in more detail around equation (25) below.

Equation (15) is a linear SDDE which can be solved analytically for $r\in (-\infty ,\infty )$. In our setting, $r\geqslant 0$ and thus we should solve equation (15) with a reflecting boundary at the origin. However, since we have assumed that $| r(t-\tau )-R| /r(t)\ll 1$, we already work in the regime where r only seldom significantly deviates from R and thus the solution of equation (15) on the full real axis should approximate well the desired solution on the positive half-line. The solution of equation (15) for $r\in (-\infty ,\infty )$ and $t\geqslant 0$ in terms of deviations of the bond length from its equilibrium length (which we call as shifted bond length),

$$\begin{eqnarray}&&x(t)=r(t)-R,\end{eqnarray} \tag{ 16 }$$

can be derived by several methods. We review two of them (time-convolutionless transform and Gaussian ansatz) for a general multidimensional linear SDDE in appendices A–B. Here, we present just the main formulas. Assuming that the system was initially in state x(0) = x₀ and that x(t) = 0 for t < 0, the formal solution of equation (15) for $r\in (-\infty ,\infty )$ and $t\geqslant 0$ reads

$$\begin{eqnarray}&&x(t)={x}_{0}\lambda (t)+\sqrt{4D}{\int }_{0}^{t}{\rm{d}}s\ \lambda (t-s){\eta }_{{\rm{r}}}^{r}(s),\end{eqnarray} \tag{ 17 }$$

where the dimensionless Green's function

$$\begin{eqnarray}&&\lambda (t)=\displaystyle \sum _{k=0}^{\infty }\displaystyle \frac{{\left(-{\omega }_{\tau }\right)}^{k}}{k!}{\left(t-k\tau \right)}^{k}\theta (t-k\tau )\end{eqnarray} \tag{ 18 }$$

solves equation (15) with $D=0$, $\lambda (t)=0$ for t < 0 and $\lambda (0)=1$. The symbol θ(x) in equation (18) stands for the Heaviside step function. For an arbitrary initial condition x(t), t < 0 and x(0) = x₀, the expression x₀λ(t) in equation (17) must be substituted by ${x}_{0}\lambda (t)-{\omega }_{\tau }{\int }_{-\tau }^{0}{\rm{d}}s\ \lambda (t-\tau -s)x(s)$. Based on the value of the reduced delay ω_ττ which is a dimensionless measure for the relevance of the delay relative to the natural relaxation time, the Green's function λ(t) in equation (17) exhibits three different dynamical regimes discussed in detail in appendix A, in figure A1 and also below: (i) monotonic exponential decay to zero for short delays $0\leqslant {\omega }_{\tau }\tau \leqslant 1/e$, (ii) oscillatory exponential decay to zero for intermediate ('resonant') delays $1/e\leqslant {\omega }_{\tau }\tau \leqslant \pi /2$, and (iii) oscillatory exponential divergence for long delays ω_ττ > π/2.

The stochastic process x(t) in equation (17) arises as a linear combination of white noises and thus the corresponding PDFs must be Gaussian. Indeed, we find that one- and two-time conditional PDFs for x(t) with the initial condition δ(x) for t < 0 and δ(x − x₀) at t = 0 read

$$\begin{eqnarray}&&{P}_{1}(x,t| {x}_{0},0)=\displaystyle \frac{1}{\sqrt{2\pi \nu (t)}}\exp \left\{-\displaystyle \frac{1}{2}{\left(\displaystyle \frac{x-\mu (t)}{\sqrt{\nu (t)}}\right)}^{2}\right\},\end{eqnarray} \tag{ 19 }$$

$$\begin{eqnarray}\begin{array}{rcl}{P}_{2}(x,t;x^{\prime} ,t^{\prime} | {x}_{0},0) & = & \displaystyle \frac{1}{\sqrt{4{\pi }^{2}\nu (t)\nu (t^{\prime} )[1-{w}^{2}(t^{\prime} ,t)]}}\\ & & \times \exp \left\{\displaystyle \frac{1}{2\left[1-{w}^{2}(t^{\prime} ,t)\right]}\left[{\left(\displaystyle \frac{x-\mu (t)}{\sqrt{\nu (t)}}\right)}^{2}+{\left(\displaystyle \frac{x^{\prime} -\mu (t^{\prime} )}{\sqrt{\nu (t^{\prime} )}}\right)}^{2}\right.\right.\\ & & -\,\left.\left.2w(t,t^{\prime} )\left(\displaystyle \frac{x-\mu (t)}{\sqrt{\nu (t)}}\right)\left(\displaystyle \frac{x^{\prime} -\mu (t^{\prime} )}{\sqrt{\nu (t^{\prime} )}}\right)\right]\right\},\end{array}\end{eqnarray} \tag{ 20 }$$

where $t\geqslant t^{\prime} \geqslant 0$ and

$$\begin{eqnarray}&&\mu (t)\equiv \left\langle x(t)\right\rangle ={x}_{0}\lambda (t),\end{eqnarray} \tag{ 21 }$$

$$\begin{eqnarray}&&\nu (t)\equiv \left\langle {x}^{2}(t)\right\rangle -{\mu }^{2}(t)=4D{\int }_{0}^{t}{\rm{d}}s\ {\lambda }^{2}(s),\end{eqnarray} \tag{ 22 }$$

$$\begin{eqnarray}&&w(t,t^{\prime} )\equiv \displaystyle \frac{\left\langle x(t)x(t^{\prime} )\right\rangle -\mu (t)\mu (t^{\prime} )}{\sqrt{\nu (t)\nu (t^{\prime} )}}=\displaystyle \frac{4D}{\sqrt{\nu (t)\nu (t^{\prime} )}}{\int }_{0}^{t^{\prime} }{\rm{d}}s\ \lambda (t-s)\lambda (t^{\prime} -s)\end{eqnarray} \tag{ 23 }$$

denote the mean of the shifted bond-length (16), its variance, and normalized time correlation, respectively.

The function ${P}_{1}(x,t| {x}_{0},0){\rm{d}}x$ stands for the probability that the system which departs with certainty from state x₀ at time 0 is found at time t somewhere in the interval (x, x + dx). Similarly, ${P}_{2}(x,t;x^{\prime} ,t^{\prime} | {x}_{0},0){\rm{d}}x{\rm{d}}x^{\prime}$ denotes the probability that the (shifted) bond length is in the interval $(x^{\prime} ,x^{\prime} +{\rm{d}}x^{\prime} )$ at time $t^{\prime}$ and at a later time t in (x, x + dx) under the condition that it has started at time t = 0 at x₀. The one-time PDF P₁ possesses the same structure as the corresponding PDF for τ = 0 [63]. The non-Markov character of the process (15) with nonzero delay manifests itself in the fact that the two-time PDF P₂ cannot be constructed from the one-time PDF P₁, while this is always possible for a Markov process.

The FPEs corresponding to the PDFs (19) and (20) are given by equations (B.6) and (B.7) in appendix B, respectively. Interestingly enough, the diffusion and drift terms in the FPEs are given by the natural scales 2D and ω_τ only in the limit τ → 0. Moreover, both coefficients acquire a time dependence, determined by the function λ(t). Specifically, the diffusion and drift coefficients in equation (B.6) for P₁ read ${D}_{\tau }(t)\,=D{\lambda }^{2}(t){\rm{d}}\left[{\int }_{0}^{t}{\rm{d}}s{\lambda }^{2}(s){/\lambda }^{2}(t)\right]/{\rm{d}}t$ and ${\omega }_{\tau }(t)=-\dot{\lambda }(t)/\lambda (t)$, respectively⁶ . While the drift coefficient in equation (B.7) for P₂ is also given by ω_τ(t), the diffusion coefficient reads $2{D}_{\tau }(t)+4D\lambda (t){\int }_{0}^{t^{\prime} }{\rm{d}}s\,{\rm{d}}\left[\lambda (t-s)\lambda (t^{\prime} -s)/\lambda (t)\right]/{\rm{d}}t$. This difference in diffusion coefficients is the reason why the PDF P₂ can be constructed from P₁ in the standard way for Markov processes only for τ = 0, where both diffusion coefficients coincide.

For τ = 0 the PDF P₁ always eventually relaxes to a time independent stationary state which does not depend on the initial condition and which is described by the equilibrium Gibbs formula ${P}_{1}(x,\infty ;{x}_{0},0)\propto \exp \left(-\omega {x}^{2}/2D\right)$. For a nonzero delay in the regimes (i) and (ii), i.e. when the noiseless solution (18) and thus the mean value μ(t) = x₀λ(t) of x converges to 0 for $t\to \infty$, the system relaxes into a time-independent non-equilibrium steady state ${P}_{1}(x,\infty ;{x}_{0},0)\propto \exp \left[-{\omega }_{\tau }(\infty ){x}^{2}/{D}_{\tau }(\infty )\right]$ with ${\omega }_{\tau }(\infty )\ne \omega$ and ${D}_{\tau }(\infty )\ne 2D$, see figure 8 in section 5.1. In these cases, our approximate model thus predicts that, in the long run, the delay merely 'deforms' the (approximate) Gaussian equilibrium distribution through a parameter renormalization. This state is characterized by a nonzero entropy production rate [56]. For long time delays, no stationary state exists. In comparison to the analogous setting with a piece-wise constant force discussed previously [11], this destabilization for long delay times τ is a new feature, due to increasingly high systematic forces that may occur for long delays.

In the regimes (i) and (ii), the variance $\nu (t)=\sqrt{\left\langle {x}^{2}(t)\right\rangle -{\left\langle x(t)\right\rangle }^{2}}$ converges to the stationary value [30, 31, 34]

$$\begin{eqnarray}&&{\nu }_{\mathrm{ss}}=\mathop{\mathrm{lim}}\limits_{t\to \infty }\nu (t)=\displaystyle \frac{2D}{{\omega }_{\tau }}\displaystyle \frac{1+\sin ({\omega }_{\tau }\tau )}{\cos ({\omega }_{\tau }\tau )}.\end{eqnarray} \tag{ 24 }$$

The derivation of this formula is given in appendix C, where we also derive an analytical expression for the stationary time correlation function $C(t)={\mathrm{lim}}_{s\to \infty }\left\langle x(s)x(s+t)\right\rangle$. Note that the variance (24) diverges upon entering the unstable regime (iii) for ${\omega }_{\tau }\tau \to \pi /2$.

The formula (24) finally allows us to specify the regime ${\nu }_{\mathrm{ss}}\ll {R}^{2}$ where the approximation $\left[r(t-\tau )-R\right]/r(t)\approx 0$ used in the derivation of equation (15) from equations (12) and (13) is reasonable because the PDF for r is relatively sharply peaked around the mean bond length R. As we already noted, equation (15) is also valid in the case $2\tau D/{R}^{2}\ll 1$ when the bond rotates only slightly in each delay interval and thus the angle φ(t, t − τ) in equations (12) and (13) is small. However, also in this case we need to additionally assume that ${\nu }_{\mathrm{ss}}\ll {R}^{2}$ in order to ensure that the error caused by considering the wrong boundary condition at r = −R is negligible. Altogether, the used approximation is expected to give good results if the condition

$$\begin{eqnarray}&&{\nu }_{\mathrm{ss}}\ll {R}^{2}\end{eqnarray} \tag{ 25 }$$

is fulfilled.

An example of the stochastic evolution of the dimer obtained from BD simulations of the exact system (12) and (13) is depicted in figure 2(a). In figure 2(b), we compare the results obtained from BD simulations with the time evolution of the average shifted bond length (21) for different values of the equilibrium length R. As expected, the approximate analytical formula (21) describes well the exact result for large enough R satisfying the inequality (25). For larger values of ν_ss/R², the analytical result underestimates the correct value. This is because the bond length in the BD simulation is obtained from equation (12) with the reflecting boundary at the origin, while we allow negative values of r(t) in the approximate analytical description. Similarly as for the mean value, the analytical formula (22) for the bond length variance ν(t) approximates very well the value obtained from BD simulations for large enough R, as shown in figure 2(c). In figure 2(d), we depict the monotonous rapid divergence of the stationary variance (24) as the time delay τ reaches the unstable regime (iii). This means that the delay tends to delocalize structures. On the other hand, the variance can be optimized as a function of the frequency ω_τ. The best localized structure is obtained for the frequency fulfilling the formula $\cos ({\omega }_{\tau }\tau )={\omega }_{\tau }\tau$, i.e. ω_τ ≈ 0.74/τ and thus $2k\tau /\gamma \approx 0.74\exp (2D\tau /{R}^{2})$. For the corresponding value 0.74 of the reduced delay ω_ττ the system is in the dynamical regime (ii) with exponentially decaying oscillations.

Zoom In Zoom Out Reset image size

2.3. Trimer

Let us now consider the system composed of three particles. Then, equation (3) gives the system of six coupled equations of motion:

$$\begin{eqnarray*}\begin{array}{rcl}{\dot{{\bf{r}}}}_{1}(t) & = & \displaystyle \frac{k}{\gamma }\left\{\left[{r}_{12}(t-\tau )-R\right]{{\bf{e}}}_{21}(t-\tau )+\left[{r}_{13}(t-\tau )-R\right]{{\bf{e}}}_{31}(t-\tau )\right\}+\sqrt{2D}{{\boldsymbol{\eta }}}_{1}(t),\\ {\dot{{\bf{r}}}}_{2}(t) & = & \displaystyle \frac{k}{\gamma }\left\{\left[{r}_{12}(t-\tau )-R\right]{{\bf{e}}}_{12}(t-\tau )+\left[{r}_{32}(t-\tau )-R\right]{{\bf{e}}}_{32}(t-\tau )\right\}+\sqrt{2D}{{\boldsymbol{\eta }}}_{2}(t),\\ {\dot{{\bf{r}}}}_{3}(t) & = & \displaystyle \frac{k}{\gamma }\left\{\left[{r}_{32}(t-\tau )-R\right]{{\bf{e}}}_{23}(t-\tau )+\left[{r}_{31}(t-\tau )-R\right]{{\bf{e}}}_{13}(t-\tau )\right\}+\sqrt{2D}{{\boldsymbol{\eta }}}_{3}(t).\end{array}\end{eqnarray*}$$

For the relative coordinates r₁₂(t) = r₁(t) − r₂(t) we obtain

$$\begin{eqnarray}\begin{array}{rcl}{\dot{{\bf{r}}}}_{12}(t) & = & -\omega \left\{\left[{r}_{12}(t-\tau )-R\right]{{\bf{e}}}_{12}(t-\tau )+\displaystyle \frac{1}{2}\left[{r}_{13}(t-\tau )-R\right]{{\bf{e}}}_{13}(t-\tau )\right.\\ & & \left.+\displaystyle \frac{1}{2}\left[{r}_{32}(t-\tau )-R\right]{{\bf{e}}}_{32}(t-\tau )\right\}+\sqrt{2D}\left[{{\boldsymbol{\eta }}}_{1}(t)-{{\boldsymbol{\eta }}}_{2}(t)\right],\end{array}\end{eqnarray} \tag{ 26 }$$

where $\omega =2k/\gamma$ and similarly for ${\dot{{\bf{r}}}}_{13}(t)$ and ${\dot{{\bf{r}}}}_{32}(t)$. To get analytical results for bond lengths ${r}_{{ij}}(t)=| {{\bf{r}}}_{{ij}}(t)|$, we multiply the formulas for ${\dot{{\bf{r}}}}_{{ij}}(t)$ by the corresponding unit vectors ${{\bf{e}}}_{{ij}}(t)={{\bf{r}}}_{{ij}}(t)/| {{\bf{r}}}_{{ij}}(t)|$ and linearize the resulting equations. To this end, we need to deal with the expressions e_α(t − τ) · e_β(t), α, β = I, ..., III, where we introduced roman numbers as a shorthand indexing I ≡ 12, II ≡ 13, III ≡ 32. For a vanishing delay τ = 0, e_α · e_α = 1 and the scalar products e_α · e_β describe the angles of the triangle formed by the three particles (see figure 1 in section 2). We find that up to the leading order in the equilibrium bond length R the triangle is equilateral and thus the internal angles are π/3, leading to the relations e_I · e_II = e_I · e_III = −e_II · e_III ≈ 1/2 + O[(r_α − r_β)/R] with a correction that is on the order of (r_I − r_II)/R for e_I · e_II and similarly for the other scalar products. The linearized equation for the relative coordinate x_α ≡ r_α − R thus reads

$$\begin{eqnarray}\begin{array}{rcl}{\dot{x}}_{\alpha }(t) & = & -\omega \left[{x}_{\alpha }(t-\tau )+\displaystyle \frac{1}{4}{x}_{\alpha +{\rm{I}}}(t-\tau )+\displaystyle \frac{1}{4}{x}_{\alpha +\mathrm{II}}(t-\tau )\right]+\sqrt{2D}\left(\displaystyle \sum _{i=1}^{3}{{ \mathcal A }}_{\alpha i}{{\boldsymbol{\eta }}}_{i}(t)\right)\cdot {{\bf{e}}}_{\alpha }(t),\\ { \mathcal A } & = & \left(\begin{array}{ccc}1 & -1 & 0\\ 1 & 0 & -1\\ 0 & -1 & 1\end{array}\right),\end{array}\end{eqnarray} \tag{ 27 }$$

where the lower index $\alpha \equiv \alpha \,\mathrm{mod}\,\mathrm{III}$ is considered as periodic with the period 3, i.e. x_IV ≡ x_I and x_V ≡ x_II. Similarly as in the case of the dimer, equation (27) describes the dynamics of x_α(t) well for large equilibrium bond lengths R and for time delays small compared to reorientation times of the unit vectors e_α.

For an analytical treatment, it is advantageous to rewrite the system (27) in the matrix form

$$\begin{eqnarray}&&\dot{{\bf{x}}}(t)=-\omega M{\bf{x}}(t-\tau )+\sqrt{2D}{\boldsymbol{\xi }}(t),\end{eqnarray} \tag{ 28 }$$

for the three-dimensional column vector ${\bf{x}}(t)={\left[{x}_{{\rm{I}}}(t),{x}_{\mathrm{II}}(t),{x}_{\mathrm{III}}(t)\right]}^{\top }$. In equation (28), the noise vector ${\boldsymbol{\xi }}(t)$ is given by ${\boldsymbol{\xi }}(t)\equiv {A}^{1}{{\boldsymbol{\eta }}}^{1}(t)+{A}^{2}{{\boldsymbol{\eta }}}^{2}(t)$ with the auxiliary noise vectors ${{\boldsymbol{\eta }}}^{j}(t)\equiv {[{\eta }_{1}^{j},{\eta }_{2}^{j},{\eta }_{3}^{j}]}^{\top }$, j = 1, 2, containing the jth components of the original noises ${{\boldsymbol{\eta }}}_{i}(t)$, i = 1, ..., 3. From the system (27) follows that the matrices M, A¹ and A² read

$$\begin{eqnarray}M=\displaystyle \frac{1}{4}\left(\begin{array}{ccc}4 & 1 & 1\\ 1 & 4 & 1\\ 1 & 1 & 4\end{array}\right),{A}^{1}=\left(\begin{array}{ccc}{e}_{{\rm{I}}}^{1} & -{e}_{{\rm{I}}}^{1} & 0\\ {e}_{\mathrm{II}}^{1} & 0 & -{e}_{\mathrm{II}}^{1}\\ 0 & -{e}_{\mathrm{III}}^{1} & {e}_{\mathrm{III}}^{1}\end{array}\right),{A}^{2}=\left(\begin{array}{ccc}{e}_{{\rm{I}}}^{2} & -{e}_{{\rm{I}}}^{2} & 0\\ {e}_{\mathrm{II}}^{2} & 0 & -{e}_{\mathrm{II}}^{2}\\ 0 & -{e}_{\mathrm{III}}^{2} & {e}_{\mathrm{III}}^{2}\end{array}\right),\end{eqnarray} \tag{ 29 }$$

where ${e}_{\alpha }^{j}$ denote jth component of the two-dimensional unit vector e_α. The time-correlations between the three components of the noise vector ${\boldsymbol{\xi }}(t)$ are not mutually independent and read

$$\begin{eqnarray}&&\left\langle {\xi }_{\alpha }(t){\xi }_{\beta }(t^{\prime} )\right\rangle ={\left({A}^{1}{A}^{1\top }+{A}^{2}{A}^{2\top }\right)}_{\alpha \beta }{\delta }_{{\rm{D}}}(t-t^{\prime} )=2{M}_{\alpha \beta }{\delta }_{{\rm{D}}}(t-t^{\prime} ).\end{eqnarray} \tag{ 30 }$$

Due to the linearity of equation (28) and Gaussianity of the noise, the Green's function for the one-time PDF ${P}_{1}({\bf{r}},t| {{\bf{r}}}_{0},0)$ is Gaussian [37], and determined by the mean value ${\boldsymbol{\mu }}(t)=\left\langle {\bf{x}}(t)\right\rangle$ and the covariance matrix ${ \mathcal K }(t)=\left\langle {\bf{x}}(t){{\bf{x}}}^{\top }(t)\right\rangle -{\boldsymbol{\mu }}(t){\left[{\boldsymbol{\mu }}(t)\right]}^{\top }$. We review in detail the derivation of these functions in appendices A and B.2.

For the initial condition x(t) = 0, t < 0 and x(0) = x₀ we get

$$\begin{eqnarray}&&{\boldsymbol{\mu }}(t)=\lambda (t){{\bf{x}}}_{0},\end{eqnarray} \tag{ 31 }$$

$$\begin{eqnarray}&&{ \mathcal K }(t)=4{DM}{\int }_{0}^{t}{\rm{d}}{s}{\lambda }^{2}(s),\end{eqnarray} \tag{ 32 }$$

where λ(t) denotes the Green's function for equation (28) given by

$$\begin{eqnarray}&&\lambda (t)=\displaystyle \sum _{k=0}^{\infty }\displaystyle \frac{{\left(-\omega M\right)}^{k}}{k!}{\left(t-k\tau \right)}^{k}\theta (t-k\tau ).\end{eqnarray} \tag{ 33 }$$

Multiplying λ(t) by the vector ${[1,1,1]}^{\top }$, using the formula $M{[1,1,1]}^{\top }=3/2{[1,1,1]}^{\top }$ and comparing the result to the one-dimensional Green's function (18), we find that λ(t) undergoes with increasing t (i) monotonous exponential decay to 0 for $0\lt 3\omega \tau /2\leqslant 1/e$, (ii) oscillatory exponential decay to 0 for 1/e < 3ωτ/2 < π/2 and (iii) oscillatory exponential divergence for π/2 < 3ωτ/2. In the regimes (i) and (ii) the stationary value of the covariance matrix reads

$$\begin{eqnarray}&&{{ \mathcal K }}_{{ss}}=\mathop{\mathrm{lim}}\limits_{t\to \infty }{ \mathcal K }(t)=\displaystyle \frac{2D}{\omega }\displaystyle \frac{{ \mathcal I }+\sin \left(\omega \tau M\right)}{\cos \left(\omega \tau M\right)},\end{eqnarray} \tag{ 34 }$$

where ${ \mathcal I }$ denotes the identity matrix. This formula follows from the results of appendix C after substituting the matrices ω and $\sigma {\sigma }^{\top }$ from the formula (C.7) in the appendix by ωM and 4DM. In the appendix, we also derive an analytical expression for the stationary time correlation matrix $C(t)={\mathrm{lim}}_{s\to \infty }\left\langle {\bf{x}}(s){{\bf{x}}}^{\top }(s+t)\right\rangle$. The regime of stability 3ωτ/2 < π/2 of the trimer can also be determined from the condition that the matrix $\cos \left(\omega \tau M\right)$ is not singular, i.e. its determinant ${\cos }^{2}\left(3\omega \tau /2\right)\cos \left(3\omega \tau \right)$ is nonzero.

Due to the symmetry of the problem, all diagonal elements of the matrix ${{ \mathcal K }}_{{ss}}$ are identical and the same holds also for all its off-diagonal elements. The approximate analytical, time-dependent solution (32) for the covariance matrix is compared to the exact covariance matrix obtained by BD simulations of the complete model in figure 3. Given the approximations made, we find very good agreement. The analytical results only slightly underestimate the diagonal elements (probably for the same reason as for the dimer) and overestimate the off-diagonal elements. The behavior of the covariance matrix as a function of the frequency ω and delay τ is similar to the behavior of the variance (24) for the dimer. Specifically, the diagonal elements of ${{ \mathcal K }}_{{ss}}$ monotonicly increase (the PDF for the bond lengths become broader) with τ and exhibit a minimum as functions of ω, opening a possibility to optimize the width of the bond length PDF. The off-diagonal elements of ${{ \mathcal K }}_{{ss}}$ monotonously increase (the individual bonds of the trimer become more correlated) both with the delay and with the natural relaxation frequency.

Zoom In Zoom Out Reset image size

The approximate analytical study of the dimer and trimer revealed that both systems obey three dynamical regimes: (i) and (ii) a monotonous and an oscillatory exponential relaxation towards a steady state with the average bond length $\mu (\infty )=R$, respectively, and (iii) an exponential divergence. The performed BD simulations confirmed that for large R, when the model is well described by the approximate analytical formulas, these regimes can indeed be observed also in the complete model (3). Furthermore, the analytical study predicted that the dimer is in the unstable regime (iii) for ωτ > π/2 and the trimer for ωτ > π/3. Let us now discuss how general the presented findings are.

The stationary average bond length μ can be determined by minimizing the potential energy $V=\tfrac{1}{2}{\sum }_{(i,j)}{V}_{2}({r}_{{ij}})$ with the two-particle potential V₂ given by equation (1). Minimizing the potential in our two-dimensional geometry yields the highly symmetric molecule-like structures shown in figure 4(a). Due to the confinement to 2d, the global minimum of the potential corresponding to μ = R is accessible only for the dimer (N = 2) and the trimer (N = 3). For larger molecules, the average bond length decreases as a result of the infinite range of the potential. The system asymptotically relaxes to the depicted structures if the noise D vanishes and the reduced delay time ωτ is small enough such that the system is in the dynamical regime (i) or (ii). Nonzero D leads to fluctuations around the asymptotic structures and large ωτ causes exponentially diverging oscillations.

Zoom In Zoom Out Reset image size

We have solved the complete model using BD simulations and depict the behavior of the average bond length r(t) for several values of N in the dynamical regimes (i) and (ii)–(iii) in figures 4(b) and (c), respectively. In the regime (i), we observe that larger systems relax faster than those with smaller N. Furthermore, in figure 4(c), we see that larger systems oscillate with larger amplitudes and that the threshold between the regimes (ii) and (iii) is reached at smaller values of ωτ. More precisely, all the curves in figure 4(c) are plotted using the same parameters except for N and, while the curves for N ≤ 3 are in the regime (ii), the curves for N > 3 correspond to the regime (iii). These observations are in accord with our analytical findings for the dimer and trimer.

By analyzing the mean bond length at late times, we have evaluated the critical reduced delay ${\left(\omega \tau \right)}^{\star }$ determining the threshold between the regimes (ii) and (iii) for N = 2, ..., 10. In figure 5, we show its rescaled value

$$\begin{eqnarray}&&{c}_{\mathrm{crit}}\equiv \displaystyle \frac{2}{\pi }{\left(\omega \tau \right)}^{\star },\end{eqnarray} \tag{ 35 }$$

where the coefficient 2/π is introduced for the comparison to the approximate result for the dimer. We find that the stability factor is well described by C/N as suggested by the approximate analytical results for the dimer $\left[{\left(\omega \tau \right)}^{\star }=\pi /2\right]$ and trimer [${\left(\omega \tau \right)}^{\star }=\pi /3$]. However, the analytical results would imply that the constant C equals to 2, which is smaller than the value C ≈ 3 obtained from the complete model. The actual dimer and trimer are thus more stable than their linearized versions considered in our analytical study. The found scaling c_crit ≈ C/N implies that the stability of a system with rescaled potential stiffness k → k/N (to render the energy extensive [81]) would be (almost) independent of the particle number.

Zoom In Zoom Out Reset image size

To better understand this behavior, let us first consider the approximate analytical model described in section 2.2. Imagine a particle in a harmonic potential centered around x = 0, which is initially located at x(0) = x₀ > 0. Assuming that x(t) = 0 for t < 0 and neglecting the noise, the particle does not feel any force in the time interval t ∈ [0, τ] such that x(t) = x₀ for all t ∈ [0, τ]. In the subsequent time interval t ∈ [τ, 2τ], the particle experiences the force F = −ω_τx₀ pushing it towards the opposite wall of the trap. For times t > 2τ the force starts changing dynamically according to the earlier position at time t − τ. The particle keeps its direction of motion until it reaches the position x₁ where the force changes its sign. For large delays, the particle may stop significantly later than crossing the minimum, so that x₁ < 0 and $| {x}_{1}| \gt | {x}_{0}|$. A similar process then repeats when the particle returns back, with the difference that now it reaches a maximum position x₂ > 0, $| {x}_{2}| \gt | {x}_{1}| \gt | {x}_{0}|$, etc. The amplitude thus increases after each half-period of oscillation causing a diverging behavior.

In order to understand the difference between the approximate analytical and the complete (numerical) solutions of the model, it is helpful to project the latter to one dimension, where a particle moves in the double-well potential depicted figure 7 in section 5.1. We assume that the particle starts in the right well and oscillates with increasing amplitude as discussed above. After some time the amplitude becomes large enough that the particle crosses the barrier to the left well. Due to the presence of the additional well, the potential now contains much wider low-energy region compared to the purely harmonic case. The particle needs longer time to travel from one (unbounded) side of the potential to its other side, and hence also the (reduced) delay ωτ required for inducing diverging oscillations is larger than in the harmonic case. As a consequence, the complete model is seen to be more stable than foreseen by our analytical considerations. Moreover, our discussion reveals the existence of a fourth dynamical regime, preceding the unstable regime (iii), where the particle hops between the individual wells of the potential and the amplitude of the oscillations remains finite.

Compared to the quasi-constant velocity model investigated in [11], our analysis thus reveals two qualitative differences. First, the structures formed in the quasi-constant velocity model oscillate for arbitrary nonzero delay τ, while in the harmonic model these oscillations appear only if ωτ is large enough. Second, the amplitude of the oscillations in the quasi-constant velocity model is always constant in time, while the osculations in the harmonic model either vanish with time, if the system is in the regime (ii), or explode with time, if the system is in the regime (iii). The behavior observed in the harmonic model can be traced back to the increase of the force with the particle distance and thus we expect an analogous behavior also for other systems with time-delayed forces increasing with distance.

The investigated model, much as the model discussed in [11], is inspired by self-organized systems, where a feedback based on the information about the state of the system at a previous time leads to structure formation. Interpreting the delayed interactions in our model as a result of such feedback control, we can investigate the entropy flow out of the system caused by the feedback. Due to the non-analyticity of the model with quasi-constant forces considered in [11], the analysis of entropy flows in the supplementary material therein was performed for vanishing delay only. Using the approximate Gaussian PDFs found in sections 2.2 and 2.3 for the dimer and the trimer, respectively, we can repeat this analysis with nonzero delay.

Without feedback, i.e. without the time-delayed harmonic interactions (2), the particles would spread diffusively and the system entropy would increase accordingly. The feedback control utilizes the information about the particle positions to drive the system into a non-equilibrium steady state with a time independent PDF P(x) and thus with a time constant configurational entropy ${ \mathcal S }=-{k}_{{\rm{B}}}\int {\rm{d}}{\bf{x}}P({\bf{x}})\mathrm{log}P({\bf{x}})$. The smaller the entropy of the non-equilibrium steady state, the more localized the steady-state structure and thus the better the result of the feedback control. Another measure of the performance of the feedback is the rate $-{\dot{{ \mathcal S }}}_{-}$ of entropy taken from the system per unit time that can also be interpreted as the amount of information pumped into the system per unit time by the feedback device. This entropy flow balances the diffusive spreading in the steady state and is thus moreover a measure of the useful 'work' (in units of J K⁻¹) performed by the feedback device against thermal dispersal. Evaluating the stationary entropy production ${\dot{{ \mathcal S }}}_{{\rm{F}}}$ due to the feedback control mechanism and the stationary entropy production ${\dot{{ \mathcal S }}}_{{\rm{D}}}$ due to the breaking of the fluctuation-dissipation theorem in equations (3) and (4) for τ > 0, one can define the feedback efficiency as the ratio ${\eta }_{{\rm{F}}}=-{\dot{{ \mathcal S }}}_{-}/({\dot{{ \mathcal S }}}_{{\rm{F}}}+{\dot{{ \mathcal S }}}_{{\rm{D}}})$. The entropy production ${\dot{{ \mathcal S }}}_{{\rm{D}}}$ can be calculated along the lines of [56]. The entropy production ${\dot{{ \mathcal S }}}_{{\rm{F}}}={q}_{{\rm{H}}}/T$ is the housekeeping heat flux q_H flowing to the bath at temperature T, due to the overall operation of the feedback device, divided by the bath temperature. It clearly depends on the specific technical realization of the feedback. In all known relevant realizations of the feedback in microscopic systems [11, 13, 16, 18, 19], the housekeeping heat flux is very large compared to the 'functional' energy fluxes in the controlled system, resulting in a large ${\dot{{ \mathcal S }}}_{{\rm{F}}}$ compared to $-{\dot{{ \mathcal S }}}_{-}$, so that the efficiency η_F of such devices is usually negligibly small.

To evaluate the entropy flow due to the feedback (the time-delayed harmonic interaction) in the present setup, we proceed along the similar lines as in [11, 64]. The center of mass coordinate of the system is not affected by the feedback and diffuses freely (see section 2.1). The structure formation due to the feedback thus occurs only on the level of the bonds. Let us now consider the time-dependent PDF P(x, t) for the bonds that converges to a time-independent non-equilibrium steady state due to the competition between feedback and diffusion. The rate of change of its Shannon entropy ${ \mathcal S }(t)=-{k}_{{\rm{B}}}\int {\rm{d}}{\bf{x}}P({\bf{x}},t)\mathrm{log}P({\bf{x}},t)$ can formally be written as

$$\begin{eqnarray}&&\dot{{ \mathcal S }}(t)={\dot{{ \mathcal S }}}_{+}(t)+{\dot{{ \mathcal S }}}_{-}(t),\end{eqnarray} \tag{ 36 }$$

where ${\dot{{ \mathcal S }}}_{+}(t)$ stands for the positive entropy flowing into the system due to the diffusive spreading of the particles and ${\dot{{ \mathcal S }}}_{-}(t)$ corresponds to the outflow of entropy due to the feedback.

Assuming that the stochastic dynamics of the column vector x(t) describing the bonds obeys the generalized Langevin equation (GLE) $\dot{{\bf{x}}}(t)={\bf{F}}[{\bf{x}}(t),{\bf{x}}(t-\tau )]+\sigma \eta (t)$, where η denotes a zero mean Gaussian white noise with the covariance matrix $\left\langle {\eta }_{i}(t){\eta }_{j}(t^{\prime} )\right\rangle ={\delta }_{{ij}}\delta (t-t^{\prime} )$, the dynamical equation for P(x, t) can be written in the form [30, 65]

$$\begin{eqnarray}&&\displaystyle \frac{\partial }{\partial t}P({\bf{x}},t)=\displaystyle \frac{1}{2}\displaystyle \sum _{{ij}}{\left(\sigma {\sigma }^{\top }\right)}_{{ij}}\displaystyle \frac{{\partial }^{2}}{\partial {x}_{i}\partial {x}_{j}}P({\bf{x}},t)+{ \mathcal L }[{\bf{x}},t].\end{eqnarray} \tag{ 37 }$$

In this equation, the fist term on the right stands for the diffusive spreading of the PDF. The term ${ \mathcal L }[{\bf{x}},t]$ corresponds to the time-delayed force F[x(t − τ)] in the Langevin equation and thus it describes the effect of the feedback. Its concrete form is not relevant for the discussion below and thus we refer to the works [30, 65] for more details about its structure.

Inserting equation (37) into the formal time derivative $\dot{{ \mathcal S }}(t)=-{k}_{{\rm{B}}}\int {\rm{d}}{\bf{x}}[{\partial }_{t}P({\bf{x}},t)]\mathrm{log}P({\bf{x}},t)$ of the system entropy ${ \mathcal S }(t)$, we find that

$$\begin{eqnarray}&&{\dot{{ \mathcal S }}}_{+}(t)=-\displaystyle \frac{{k}_{{\rm{B}}}}{2}\int {\rm{d}}{\bf{x}}\displaystyle \sum _{{ij}}{\left(\sigma {\sigma }^{\top }\right)}_{{ij}}\left(\displaystyle \frac{{\partial }^{2}}{\partial {x}_{i}\partial {x}_{j}}P\right)\mathrm{log}P=-{k}_{{\rm{B}}}\int {\rm{d}}{\bf{x}}{\left({\rm{\nabla }}P\right)}^{\top }\displaystyle \frac{\sigma {\sigma }^{\top }}{P}{\rm{\nabla }}P,\end{eqnarray} \tag{ 38 }$$

$$\begin{eqnarray}&&{\dot{{ \mathcal S }}}_{-}(t)=-{k}_{{\rm{B}}}\int {\rm{d}}{\bf{x}}{ \mathcal L }[{\bf{x}},t]\mathrm{log}P({\bf{x}},t)=\dot{{ \mathcal S }}(t)-{\dot{{ \mathcal S }}}_{+}(t).\end{eqnarray} \tag{ 39 }$$

The last equation allows us to calculate the amount of entropy taken out of the system due to the feedback per unit time, ${\dot{{ \mathcal S }}}_{-}(t)$, from the PDF P(x, t) without knowing the explicit form of the operator ${ \mathcal L }$. It is interesting to adopt the Seifert's idea of trajectory-dependent entropy [57, 66] and use equation (39) to define the stochastic (position dependent) entropy flux

$$\begin{eqnarray}&&{\dot{s}}_{-}({\bf{x}},t)={k}_{{\rm{B}}}\left[{\left({\rm{\nabla }}P\right)}^{\top }\displaystyle \frac{\sigma {\sigma }^{\top }}{{P}^{2}}{\rm{\nabla }}P-\displaystyle \frac{1}{2}{\partial }_{t}{\left(\mathrm{log}P\right)}^{2}\right].\end{eqnarray} \tag{ 40 }$$

The average flux (39) then follows as the average ${\dot{{ \mathcal S }}}_{-}(t)=\left\langle {\dot{s}}_{-}({\bf{x}},t)\right\rangle$ either over the PDF P(x, t) or over the individual stochastic trajectories generated in a BD simulation. To the best of our knowledge, the statistics of the entropy flux (40) has not been investigated yet and thus it is not known whether its PDF fulfills some fluctuation symmetries. Such an investigation would clearly be beyond the scope of the present paper and we leave it for a future work.

Let us now evaluate the three entropy fluxes (36), (38) and (39) for a general d-dimensional Gaussian PDF

$$\begin{eqnarray}&&P({\bf{x}},t)=\displaystyle \frac{1}{\sqrt{{\left(2\pi \right)}^{d}| { \mathcal K }(t)| }}\exp \left[-\displaystyle \frac{1}{2}{\left({\bf{x}}-\mu (t)\right)}^{\top }{{ \mathcal K }}^{-1}(t)\left({\bf{x}}-\mu (t)\right)\right],\end{eqnarray} \tag{ 41 }$$

where $| { \mathcal K }(t)|$ denotes determinant of the covariance matrix ${ \mathcal K }(t)$. The corresponding entropy ${ \mathcal S }(t)$ reads

$$\begin{eqnarray}&&\displaystyle \frac{{ \mathcal S }(t)}{{k}_{{\rm{B}}}}=\displaystyle \frac{d}{2}+\displaystyle \frac{1}{2}\mathrm{log}\left[{\left(2\pi \right)}^{d}| { \mathcal K }(t)| \right]\end{eqnarray} \tag{ 42 }$$

leading to the rate of change

$$\begin{eqnarray}&&\displaystyle \frac{\dot{{ \mathcal S }}(t)}{{k}_{{\rm{B}}}}=\displaystyle \frac{1}{2}\displaystyle \frac{{\rm{d}}}{{\rm{d}}t}\mathrm{log}| { \mathcal K }(t)| .\end{eqnarray} \tag{ 43 }$$

From equation (38) we then find that

$$\begin{eqnarray}&&\displaystyle \frac{{\dot{{ \mathcal S }}}_{+}(t)}{{k}_{{\rm{B}}}}=\displaystyle \frac{1}{2}\mathrm{Tr}\left[\sigma {\sigma }^{\top }{{ \mathcal K }}^{-1}(t)\right].\end{eqnarray} \tag{ 44 }$$

For finite times, all the entropy fluxes depend on the initial conditions and can be determined from equations (39), (43) and (44).

Let us now focus on the specific setups considered in sections 2.2 and 2.3. For the dimer, we have investigated the PDF for the length of the single bond and thus d = 1. Using our analytical findings with $\sigma {\sigma }^{\top }=4D$ and ${ \mathcal K }(t)=\nu (t)$, we get

$$\begin{eqnarray}&&\displaystyle \frac{{ \mathcal S }(t)}{{k}_{{\rm{B}}}}=1+\displaystyle \frac{1}{2}\mathrm{log}[{\left(2\pi \right)}^{2}\nu (t)]\end{eqnarray} \tag{ 45 }$$

and

$$\begin{eqnarray}&&\displaystyle \frac{\dot{{ \mathcal S }}(t)}{2{{Dk}}_{{\rm{B}}}}=\displaystyle \frac{{\lambda }^{2}(t)}{\nu (t)},\ \ \displaystyle \frac{{\dot{{ \mathcal S }}}_{+}(t)}{2{{Dk}}_{{\rm{B}}}}=\displaystyle \frac{1}{\nu (t)},\ \ \displaystyle \frac{{\dot{{ \mathcal S }}}_{-}(t)}{2{{Dk}}_{{\rm{B}}}}=\displaystyle \frac{{\lambda }^{2}(t)-1}{\nu (t)},\end{eqnarray} \tag{ 46 }$$

where λ(t) is given by equation (18) and the variance reads $\nu (t)=4D{\int }_{0}^{t}{\rm{d}}s\ {\lambda }^{2}(s)$. Similarly, in our analytical investigation of the trimer, we have fixed the angles between the individual bonds and investigated the PDF for the three bond length only, implying that d = 3. Using $\sigma {\sigma }^{\top }=4{DM}$ and Jacobi's formula for the derivative of determinants, we obtain the expressions

$$\begin{eqnarray}&&\displaystyle \frac{{ \mathcal S }(t)}{{k}_{{\rm{B}}}}=\displaystyle \frac{3}{2}+\displaystyle \frac{1}{2}\mathrm{log}[{\left(2\pi \right)}^{3}| { \mathcal K }(t)| ]\end{eqnarray} \tag{ 47 }$$

and

$$\begin{eqnarray}&&\displaystyle \frac{\dot{{ \mathcal S }}(t)}{2{{Dk}}_{{\rm{B}}}}=\mathrm{Tr}\left[\displaystyle \frac{M{\lambda }^{2}(t)}{{ \mathcal K }(t)}\right],\ \ \displaystyle \frac{{\dot{{ \mathcal S }}}_{+}(t)}{2{{Dk}}_{{\rm{B}}}}=\mathrm{Tr}\left[\displaystyle \frac{M}{{ \mathcal K }(t)}\right],\ \ \displaystyle \frac{{\dot{{ \mathcal S }}}_{-}(t)}{2{{Dk}}_{{\rm{B}}}}=\mathrm{Tr}\left[M\displaystyle \frac{{\lambda }^{2}(t)-{ \mathcal I }}{{ \mathcal K }(t)}\right],\end{eqnarray} \tag{ 48 }$$

where λ(t) is given by equation (33) and the covariance matrix ${ \mathcal K }(t)$ reads $4{DM}{\int }_{0}^{t}{\rm{d}}s{\lambda }^{2}(s)$.

The formulas (45)–(48) are valid both in the stable regimes (i) and (ii), where the system at long times relaxes to a stationary time-independent structured state, and in the unstable regime (iii). In the unstable regime, the variance ν(t) and the covariance ${ \mathcal K }(t)$ diverge in time. As a result, the system entropy ${ \mathcal S }(t)$ diverges and the entropy flow ${\dot{{ \mathcal S }}}_{+}(t)$ decays to zero, because the variance of the PDF is so large that the diffusion can hardly further increase it. On the other hand, the rate of entropy change $\dot{{ \mathcal S }}(t)$ and thus also the entropy outflow ${\dot{{ \mathcal S }}}_{-}(t)$ remain finite oscillating functions, as can be seen for the dimer by using the exponential long-time approximation (A.7) for λ(t), and similarly for the trimer.

For the purpose of structure formation, only the regimes (i) and (ii) are of interest, because only then the PDF reaches a time-independent non-equilibrium steady state at long times, i.e. ${\mathrm{lim}}_{t\to \infty }{ \mathcal S }(t)={ \mathcal S }$, ${\mathrm{lim}}_{t\to \infty }\dot{{ \mathcal S }}(t)=0$ and ${\dot{{ \mathcal S }}}_{-}\equiv {\mathrm{lim}}_{t\to \infty }{\dot{{ \mathcal S }}}_{-}(t)=-{\mathrm{lim}}_{t\to \infty }{\dot{{ \mathcal S }}}_{+}(t)$. Let us therefore now evaluate the long-time stationary system entropies S and entropy fluxes ${\dot{{ \mathcal S }}}_{-}(t)$ maintaining the molecular-like structures formed in our model for the dimer and the trimer in these two regimes. Using the asymptotic formulas (24) and (34) for the dimer bond-length stationary variance ν_ss and trimer covariance ${{ \mathcal K }}_{\mathrm{ss}}$, we find from equations (45)–(48)

$$\begin{eqnarray}&&\displaystyle \frac{{S}}{{k}_{{\rm{B}}}}=1+{\rm{l}}{\rm{o}}{\rm{g}}(2\pi )+\displaystyle \frac{1}{2}{\rm{l}}{\rm{o}}{\rm{g}}\left[\displaystyle \frac{2D}{{\omega }_{\tau }}\displaystyle \frac{1+\sin ({\omega }_{\tau }\tau )}{\cos ({\omega }_{\tau }\tau )}\right],\end{eqnarray} \tag{ 49 }$$

$$\begin{eqnarray}&&\displaystyle \frac{{\dot{{S}}}_{-}}{{k}_{{\rm{B}}}}=-\displaystyle \frac{{\omega }_{\tau }\cos \left({\omega }_{\tau }\tau \right)}{1+\sin \left({\omega }_{\tau }\tau \right)},\end{eqnarray} \tag{ 50 }$$

for the dimer and

$$\begin{eqnarray}&&\displaystyle \frac{{S}}{{k}_{{\rm{B}}}}=\displaystyle \frac{3}{2}+\displaystyle \frac{3}{2}{\rm{l}}{\rm{o}}{\rm{g}}(2\pi )+\displaystyle \frac{1}{2}{\rm{l}}{\rm{o}}{\rm{g}}\left[{\left(a-b\right)}^{2}(a+2b)\right],\end{eqnarray} \tag{ 51 }$$

$$\begin{eqnarray}&&\displaystyle \frac{{\dot{{S}}}_{-}}{{k}_{{\rm{B}}}}=-2D{\rm{T}}{\rm{r}}\left[\displaystyle \frac{M}{{{K}}_{{\rm{s}}{\rm{s}}}}\right]=-6D\displaystyle \frac{(a+b)}{(a-b)(a+2b)},\end{eqnarray} \tag{ 52 }$$

for the trimer. In the last two formulas, $a={{ \mathcal K }}_{\mathrm{ss}}(1,1)$ denotes diagonal and $b={{ \mathcal K }}_{\mathrm{ss}}(1,2)$ off-diagonal elements of the covariance matrix. The system entropies (49) and (51) are determined by the width of the PDFs for the bonds. Therefore, they monotonously increase with temperature T = γ D/k_B and with the delay time τ, and exhibit a minimum as functions of the frequencies ω_τ (dimer) and ω (trimer), similarly as the variance ν_ss and the diagonal matrix elements of the covariance matrix ${{ \mathcal K }}_{\mathrm{ss}}$. The quality of the steady-state structures is thus in our model always unfavorably influenced by the delay and, for a given delay, one can tune the frequency in order to minimize this (usually unwanted) effect.

The two entropy fluxes (50) and (52) are negative, highlighting that they correspond to entropy outflows from (or information inflows into) the system. Interestingly enough, the entropy fluxes do not depend on the temperature T (or noise strength D) as already predicted in [11]. This means that the fluxes are discontinuous in the formal limit T → 0 because they must inevitably vanish for zero noise, where the PDF for the system evolution is a δ-function for all times. We plot the stationary entropy fluxes (50) and (52) as functions of the delay τ in figure 6(a) and frequency ω in figure 6(b). Naturally, maintaining a stationary structure in a bigger system (trimer, dashed lines) requires a larger (more negative) entropy flux (or more information) than in the smaller one (dimer, solid lines). The maximum of the fluxes $| {\dot{{ \mathcal S }}}_{-}|$, depicted in figure 6(b), arises as a result of a competition between stronger confinement, corresponding to larger frequencies ω, and gradual destabilization with increasing ωτ, when the system enters the unstable regime (iii). For the figure, we used for simplicity the approximation ω_τ ≈ ω for the dimer.

Zoom In Zoom Out Reset image size

The particles within the molecular structures depicted in figure 4(a) may exchange their positions. Assuming the particles to be distinguishable, different arrangements of the same structure may arise which can be interpreted as different isomers of the same molecule. Their study can provide further insight into the stability properties of our non-equilibrium molecules. In fact, we find that the study for molecules made out of only a few particles is informative also for the phenomenology observed for large particle numbers. While for the purely deterministic motion then isomer transitions only appear for time delays τ in the unstable regime (iii), in a system affected by thermal fluctuations the transitions occur for arbitrary delays. The evaluation of the frequencies of such transitions, which measure the stability of the individual isomers, can thus provide insight into the overall energy landscape responsible for the non-equilibrium structure formation. It is the main topic of transition rate theory [23].

The transition rate ${\kappa }_{A\to B}(t)$ for switching from a conformation A to a conformation B at time t can be (for arbitrary dynamics) found from the mean number of transitions N_A→B(t) from A to B during an infinitesimal time interval Δ t as ${\kappa }_{A\to B}(t)={N}_{A\to B}(t)/{\rm{\Delta }}t$. Alternatively, one can get it from the inverse mean first passage time for changing the two isomers, leading to the same results. In general, the deduced transition rates depend on the initial state of the system and on time and they can be calculated analytically only in few simple situations. While they in principle can straightforwardly be evaluated in simulations, this can take a (forbiddingly) long time if the transition rates are small.

For an analytical treatment, it is more convenient to define the transition rate via the so-called survival probability S(t) that the system has not changed its initial isomer until time t. Our particular problem concerning the transitions between different isomers can be mapped to a particle moving in a high-dimensional energy landscape. We denote by ${S}_{A\to B}(t)$ the survival probability that the system, starting in the conformation A with an absorbing boundary at the top of the barrier to conformation B (and reflecting barriers elsewhere), will remain in the configuration A until time t. The transition rate between the states A and B is then given by ${\kappa }_{A\to B}(t)\,={\dot{S}}_{A\to B}(t)/{S}_{A\to B}(t)$. Hence, if the dynamical equation for the probability distribution for the state of the system with the correct boundary conditions is known, we can determine the transition rate numerically and, in some situations, even analytically.

Considering standard Markovian Langevin dynamics, the asymptotic form ${\mathrm{lim}}_{t\to \infty }{\kappa }_{A\to B}(t)$ of the transition rate can (approximately) be calculated using Kramers' rate theory [22, 23] which was originally developed to describe chemical reaction rates. The approximation works best for a high energy barrier compared to the thermal energy. Kramers' theory was extended to reaction rates for GLEs describing non-Markovian systems. A crucial contribution in this direction came from Grote and Hynes [67] who derived a dynamical correction to Kramers' result. While their analysis was based on a parabolic barrier, Pollak [68] investigated the decay rate of an underdamped particle trapped in a symmetric cusp double well potential obeying the GLE with an arbitrary memory kernel satisfying the fluctuation-dissipation theorem. The time-dependent rate ${\kappa }_{A\to B}(t)$ for driven overdamped systems can be calculated using the recent theory of Bullerjahn et al [24] for forcible molecular bond breaking.

To the best of our knowledge, the literature on the rate theory of time-delayed systems is scarce. The escape from a cubic metastable well under a time-delayed friction was investigated in [69]. Based on their small-delay approximation, Guillouzic et al [70] calculated the transition rate and the mean first passage time for an overdamped Brownian particle in a delayed quartic potential. From an experimental point of view, Curtin et al [71] studied transitions in a bistable system under time-delayed feedback.

Our model does not belong to any of the previously investigated classes of systems. However, for a vanishing delay one can use Kramers' theory, since the system obeys a Markovian overdamped Langevin equation. Moreover, for nonvanishing delays in the stable regimes (i) and (ii), the one-time PDFs for dimer and trimer can be described by standard (time-local) FPEs with time-dependent coefficients, where Bullerjahn's theory applies and where one can evaluate the transition rate numerically. Furthermore, after long times, the coefficients in these FPEs become time independent suggesting that Kramers' theory may be applied also to obtain the long-time form of the transition rates for a nonzero delay.

Although looking promising, all the techniques above are based on the time-local FPE. For non-zero delay, they share one drawback, which may limit their applicability to small delays: the time-local FPE is derived from solutions to the delayed Langevin equations without the absorbing boundary condition. While this represents no problem for diffusion dynamics without delay, it can cause problems in our delayed system. In the following sections, we compare predictions of Kramers' theory, Bullerjahn's theory and direct numerical evaluation of the transition rates from the time-local FPE against BD simulations of the transition rates for dimer and trimer, demonstrating that the rates obtained from the time-local FPE are indeed accurate for small and moderate delays only.

5.1. Dimer

To study transition rates, the simplest configuration of our model is the dimer with two distinguishable particles in one dimension. (Due to rotational symmetry, we cannot distinguish between dimer isomers in two dimensions.) The setting is described by the approximate Langevin equation (15) for the inter-particle distance $r=| {{\bf{r}}}_{1}-{{\bf{r}}}_{2}| \gt 0$. A transition occurs when the two particles exchange their positions, and can be assigned to the moment when the bond length r vanishes. To illustrate the problem, it is useful to extend the domain of the distance variable r such that it is positive for one isomer and negative for the other. For vanishing delay, this redefined signed bond length $\tilde{r}$ then diffuses in the cusped double-well potential $V(\tilde{r})=\gamma \omega {\left(| \tilde{r}| -R\right)}^{2}/2$, depicted in figure 7, with the diffusion coefficient 2D.

Zoom In Zoom Out Reset image size

For a nonzero delay, based on the approximate solution (17) to the Langevin equation (15) assuming $r\in (-\infty ,\infty )$ and x(t) = 0 for t < 0, we have found that the one-time PDF ${P}_{1}={P}_{1}(x,t)={P}_{1}(x,t| {x}_{0},0)$ obeys the FPE (B.6), which reads

$$\begin{eqnarray}&&\displaystyle \frac{\partial }{\partial t}{P}_{1}=\displaystyle \frac{\partial }{\partial x}\left[{\omega }_{\tau }(t)x+2{D}_{\tau }(t)\displaystyle \frac{\partial }{\partial x}\right]{P}_{1}\equiv { \mathcal L }(x,t){P}_{1}.\end{eqnarray} \tag{ 53 }$$

This equation describes diffusion in the harmonic potential

$$\begin{eqnarray}&&V(\tilde{r},t)=\gamma {\omega }_{\tau }(t){\left(| \tilde{r}| -R\right)}^{2}/2,\end{eqnarray} \tag{ 54 }$$

with the time-dependent stiffness γω_τ(t) (given by (B.4)) and the time-dependent diffusion coefficient 2D_τ(t) (given by (B.8)).

The validity of equation (53) for P₁ with natural boundary conditions suggests that one can further employ the analogy between the delayed dynamics and the (effective) Markovian model for calculating the transition rate κ(t) for switching between the isomers. In the Markovian case, the transition rate for surpassing the (effective) energy barrier γω_τ(t)R²/2 at r = 0 to the other isomer can be calculated from equation (53) with an absorbing boundary at x = −R [63]. We now review several methods suitable for this task, and compare the results to BD simulations of the complete model with energy barrier γωR²/2 and delayed dynamics. In order to study the transition rate between the isomers of the dimer analytically, it is enough to consider the dynamics of the system in one of the wells of the potential, i.e. for $x=\tilde{r}-R\gt 0$.

5.1.1. Numerical method

We first consider the situation when the system dwells in the state x(t) = 0 for $t\leqslant 0$ and then starts to diffuse in the time-dependent potential (54). Then, the time-dependent Markovian rate κ_M(t) can be determined from the equation

$$\begin{eqnarray}&&{\partial }_{t}{\tilde{P}}_{{\rm{a}}}(x,t)={ \mathcal L }(x,t){\tilde{P}}_{{\rm{a}}}(x,t)+{\kappa }_{{\rm{M}}}(t){\tilde{P}}_{{\rm{a}}}(x,t),\end{eqnarray} \tag{ 55 }$$

for the normalized PDF ${\tilde{P}}_{{\rm{a}}}(x,t)={P}_{{\rm{a}}}(x,t)/{S}_{{\rm{a}}}(t)$ for the position of the particle surviving in the right well of the cusped potential [72]. Here, P_a(x, t) is the solution to the FPE (53) with absorbing boundary at x = −R and ${S}_{{\rm{a}}}(t)={\int }_{-R}^{\infty }{\rm{d}}{{xP}}_{{\rm{a}}}(x,t)$ is the probability that the particle survives in the right well until time t. Equation (55) follows from equation (53) by inserting the definitions of the PDF ${\tilde{P}}_{{\rm{a}}}$ and of the transition rate

$$\begin{eqnarray}&&{\kappa }_{{\rm{M}}}(t)=-{\dot{S}}_{{\rm{a}}}(t)/{S}_{{\rm{a}}}(t).\end{eqnarray} \tag{ 56 }$$

We solved it numerically using the method presented in [73].

5.1.2. Bullerjahn's method

Alternatively, one can determine the rate approximately using the analytical theory developed by Bullerjahn et al in [24]. Therein, the rate is constructed from the (Gaussian) solution P₁ (20) of the FPE (53) with natural boundary conditions. Specifically, one approximates the probability current

$$\begin{eqnarray}&&j(-R,t)={\dot{S}}_{{\rm{a}}}(t)=-\left[{\omega }_{\tau }(t)x+2{D}_{\tau }(t){\partial }_{x}\right]{P}_{{\rm{a}}}(x,t){| }_{x=-R}\end{eqnarray} \tag{ 57 }$$

across the absorbing boundary by⁷

$$\begin{eqnarray}&&{j}^{* }(-R,t)\equiv -\left[{\omega }_{\tau }(t)x+4{D}_{\tau }(t){\partial }_{x}\right]{P}_{1}(x,t){| }_{x=-R},\end{eqnarray} \tag{ 58 }$$

and the survival probability S_a(t) by

$$\begin{eqnarray}&&{S}_{1}(t)={\int }_{-R}^{\infty }{\rm{d}}x\ {P}_{1}(x,t)=\displaystyle \frac{1}{2}\left[1+\mathrm{Erf}\left(R/\sqrt{2\nu (t)}\right)\right].\end{eqnarray} \tag{ 59 }$$

In the last expression, the symbol Erf denotes the error function and the variance ν(t) is given by equation (22). The approximate Markovian transition rate is then given by

$$\begin{eqnarray}&&{\kappa }_{{\rm{B}}}(t)=-{j}^{* }(-R,t)/{S}_{1}(t).\end{eqnarray} \tag{ 60 }$$

In figure 8, we show the frequency ω_τ(t), the (effective) diffusion coefficient 2D_τ(t), survival probabilities S_a(t) and S₁(t) and the transition rates κ_M(t) and κ_B(t) as functions of time t for parameters in the dynamical regime (i). One can observe that both the parameters ${\omega }_{\tau }(t)=\dot{\lambda }(t)/\lambda (t)$ and ${D}_{\tau }(t)=2D\lambda {\left(t\right)}^{2}+{\omega }_{\tau }(t)\nu (t)$ and the rates saturate with time. The relaxation time of ω_τ(t) is determined by the time in which the Green's function λ(t) approaches the long-time exponential representation (A.7), and the corresponding stationary value, ${\omega }_{\tau }(\infty )=1/{t}_{{\rm{R}}}$, is controlled by the relaxation time t_R for decay of λ(t) to 0, see also figure A1 in appendix A. The effective diffusion coefficient converges to the value ${D}_{\tau }(\infty )=2D[1+\sin ({\omega }_{\tau }\tau )]/\cos ({\omega }_{\tau }\tau )\geqslant 2D$, determined by the stationary variance $\nu (\infty )={\nu }_{\mathrm{ss}}$, see equation (24). The transition rates relax with the relaxation time t_R, similarly as the corresponding PDFs P_a and P₁. The analytical expressions for S₁(t) and κ_B(t) approximate the numerical results for S_a(t) and κ_M(t) best for short times $t\ll {t}_{{\rm{R}}}$, where the PDFs P_a and P₁ are still hardly affected by the different boundary conditions at x = −R. For long times and up to moderate values of time delay, the approximate analytical transition rate κ_B overestimates the corresponding exact rate κ_M, see also figures 9(a) and (b) below. For long delays, the (effective) barrier height over the (effective) thermal energy decreases so that the assumptions of the transition state theory are not valid any more, and κ_B < κ_M, see figure 9(c).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The situation of low barriers can be understood from the behavior at vanishing potential strength ω → 0, when D_τ(t) = 2D and the finite time transition rates κ_M(t) and κ_B(t) can be calculated analytically. Namely, the PDFs ${\tilde{P}}_{{\rm{a}}}$ and P₁ in the definitions (56) and (60) of the rates read ${\tilde{P}}_{{\rm{a}}}(x,t)=\left\{\exp (-{x}^{2}/4{Dt})-\exp \left[-{\left(x+2R\right)}^{2}/4{Dt}\right]\right\}/\sqrt{4\pi {Dt}}$ and ${P}_{1}(x,t)=\exp (-{x}^{2}/4{Dt})/\sqrt{4\pi {Dt}}$ [63] leading to the formulas

$$\begin{eqnarray}&&{\kappa }_{{\rm{M}}}(t)=\displaystyle \frac{2R}{t}\displaystyle \frac{\exp (-{R}^{2}/4{Dt})}{{\int }_{-R}^{\infty }{\rm{d}}x\left\{\exp (-{x}^{2}/4{Dt})-\exp \left[-{\left(x+2R\right)}^{2}/4{Dt}\right]\right\}},\end{eqnarray} \tag{ 61 }$$

$$\begin{eqnarray}&&{\kappa }_{{\rm{B}}}(t)=\displaystyle \frac{2R}{t}\displaystyle \frac{\exp (-{R}^{2}/4{Dt})}{{\int }_{-R}^{\infty }{\rm{d}}x\exp (-{x}^{2}/4{Dt})}.\end{eqnarray} \tag{ 62 }$$

Since the denominator in the expression for the rate κ_M(t) is smaller than that for κ_B(t), we conclude that, for low energy barriers, the inequality ${\kappa }_{{\rm{B}}}(t)\leqslant {\kappa }_{{\rm{M}}}(t)$ holds. On the other hand, for very high barriers, the PDFs ${\tilde{P}}_{{\rm{a}}}$ and P₁ are (almost) identical because the absorbing boundary at x = −R is effectively inaccessible. In such a case, the probability current j(−R, t) is (almost) zero, S(t) ≈ 1, and ${j}^{* }\approx j-2D{\partial }_{x}{P}_{1}{| }_{x=-R}\lt 0$ leads to the inequality ${\kappa }_{{\rm{B}}}(t)\geqslant {\kappa }_{{\rm{M}}}(t)$.

To gain a deeper insight into the behavior of the transition rates, let us consider the stationary regime, $t\to \infty$. In this regime, ${\tilde{P}}_{{\rm{a}}}(-R,\infty )=0$, due to the absorbing boundary condition at x = −R, and ${j}_{1}(x,\infty )\,\equiv -\left[{\omega }_{\tau }(\infty )x+2{D}_{\tau }(\infty ){\partial }_{x}\right]{\tilde{P}}_{1}(x,\infty ){| }_{x=-R}=0$, due to the conservation of probability ${\partial }_{t}{\tilde{P}}_{1}(x,\infty )\,={\partial }_{t}[{P}_{1}(x,\infty )/{S}_{1}(\infty )]=-{\partial }_{x}{j}_{1}(x,\infty )$, where ∂_t ≡ ∂/∂t and ∂_x ≡ ∂/∂x. Then the transition rates ${\kappa }_{{\rm{M}}}(\infty )\,=-{j}_{\mathrm{Da}}$ and ${\kappa }_{{\rm{B}}}(\infty )=-{j}_{{\rm{D}}1}$ are determined solely by the diffusive fluxes ${j}_{\mathrm{Da}}\,\equiv 2{D}_{\tau }(\infty ){\partial }_{x}{\tilde{P}}_{{\rm{a}}}(x,\infty ){| }_{x=-R}$ and ${j}_{{\rm{D}}1}\equiv 2{D}_{\tau }(\infty ){\partial }_{x}{\tilde{P}}_{1}(x,\infty ){| }_{x=-R}$. The smaller the frequency ω, the wider the PDFs ${\tilde{P}}_{1}$ and ${\tilde{P}}_{{\rm{a}}}$. For small ω, the boundary at x = −R is in the region where the PDF ${\tilde{P}}_{1}$ has its maximum and it is also close to the maximum of ${\tilde{P}}_{{\rm{a}}}$. In such a case, the PDF ${\tilde{P}}_{{\rm{a}}}$, which must vanish at x = −R, changes near the boundary faster than ${\tilde{P}}_{1}$, leading to $| {j}_{{\rm{D}}1}| \lt | {j}_{\mathrm{Da}}|$ and ${\kappa }_{{\rm{B}}}(\infty )\leqslant {\kappa }_{{\rm{M}}}(\infty )$, in accord with the argument put forward in the previous paragraph. With increasing ω, the boundary shifts away from the maxima of ${\tilde{P}}_{1}$ and ${\tilde{P}}_{{\rm{a}}}$ towards their tails. Due to the trajectories trapped in the absorbing state [72, 74], the maximum of ${\tilde{P}}_{{\rm{a}}}$ is slightly farther away from the absorbing boundary than the maximum of ${\tilde{P}}_{1}$, and thus, with increasing ω, the tail of ${\tilde{P}}_{{\rm{a}}}$, with small derivative (small j_Da), hits the boundary at x = −R before the corresponding tail of ${\tilde{P}}_{1}$. Hence, for large enough barrier height, the inequality between the rates crosses over to ${\kappa }_{{\rm{B}}}(\infty )\geqslant {\kappa }_{{\rm{M}}}(\infty )$. Finally, for very stiff traps ($\omega \to \infty$), both j_Da and j_D1 vanish and ${\kappa }_{{\rm{M}}}(\infty )={\kappa }_{{\rm{B}}}(\infty )=0$.

As shown in the figure 8(d), the transition rates converge with time to constant values in regime (i), where the limits ${\mathrm{lim}}_{t\to \infty }{\omega }_{\tau }(t)$ and ${\mathrm{lim}}_{t\to \infty }{D}_{\tau }(t)$ exist. Also in regime (ii), the PDF P₁ assumes, after long times, the time-independent stationary form ${P}_{1}=\exp \left(-{x}^{2}/2{\nu }_{\mathrm{ss}}\right)/\sqrt{2\pi {\nu }_{\mathrm{ss}}}$ with the variance ν_ss given by equation (24). This suggests that, also in this regime, the transition rate should saturate at long times. However, both the frequency ω_τ(t) and the diffusion coefficient 2D_τ(t) actually exhibit diverging oscillations caused by the oscillations in the Green's function λ(t) (18), in regime (ii). These divergences cause problems both in the FPE and in the approximate calculation of the rate using Bullerjahn's method. As a consequence, the (effective) Markov description can not be valid in the dynamical regime (ii). Nevertheless, let us now investigate to what extent the long-time transition rate $\kappa =\kappa (\infty )$ obtained from BD simulations of the Langevin equation (15) is captured by the predictions (56) and (60) above.

5.1.3. Long-time behavior and Kramers' method

Assuming that at long times the PDF $\tilde{P}$ is time-independent and the limits ${\mathrm{lim}}_{t\to \infty }{\omega }_{\tau }(t)$ and ${\mathrm{lim}}_{t\to \infty }{D}_{\tau }(t)$ exist, we can rewrite the formula (55) as the eigenvalue problem

$$\begin{eqnarray}&&{ \mathcal L }(x,\infty ){\tilde{P}}_{{\rm{a}}}(x,\infty )={\kappa }_{{\rm{M}}}{\tilde{P}}_{{\rm{a}}}(x,\infty ),\end{eqnarray} \tag{ 63 }$$

for the long time Markovian transition rate ${\kappa }_{{\rm{M}}}={\kappa }_{{\rm{M}}}(\infty )$. We solve this formula numerically using the method described in [72, 73]. The steady state value of the transition rate predicted with Bullerjahn's method reads

$$\begin{eqnarray}&&{\kappa }_{{\rm{B}}}=-{j}^{* }(\infty )/{S}_{1}(\infty )\end{eqnarray} \tag{ 64 }$$

with the survival probability ${S}_{1}(\infty )=\left[1+\mathrm{Erf}\left(R/\sqrt{2}{\nu }_{\mathrm{ss}}\right)\right]/2$ and the probability current ${j}^{* }(\infty )\,={\omega }_{\tau }(\infty )R\exp \left(-{R}^{2}/2{\nu }_{\mathrm{ss}}\right)/\sqrt{2\pi {\nu }_{\mathrm{ss}}}$.

For high barriers, where ${S}_{1}(\infty )\approx 1$, the long time form of Bullerjahn's transition rate coincides with the classical prediction by Kramers [22, 23] for the transition rate for leaving one of the wells of a cusped potential with barrier height E_b,

$$\begin{eqnarray}&&{\kappa }_{{\rm{K}}}=\displaystyle \frac{2}{\sqrt{\pi }{R}^{2}}\displaystyle \frac{{\left(\beta {E}_{{\rm{b}}}\right)}^{3/2}}{\beta \gamma }\exp \left(-\beta {E}_{{\rm{b}}}\right)=\displaystyle \frac{{\omega }_{\tau }(\infty )R}{\sqrt{2\pi {\nu }_{\mathrm{ss}}}}\exp \left(-\displaystyle \frac{{R}^{2}}{2{\nu }_{\mathrm{ss}}}\right)={j}^{* }(\infty ).\end{eqnarray} \tag{ 65 }$$

In the penultimate equality, we used the appropriate inverse thermal energy $\beta =1/2\gamma {D}_{\tau }(\infty )$ and barrier height ${E}_{{\rm{b}}}=\gamma {\omega }_{\tau }(\infty ){R}^{2}/2$, and also the asymptotic form of the diffusion coefficient ${D}_{\tau }(\infty )={\omega }_{\tau }(\infty ){\nu }_{\mathrm{ss}}$, which follows from the condition ${\mathrm{lim}}_{t\to \infty }\lambda (t)=0$, valid in the dynamical regimes (i) and (ii).

5.1.4. Renormalized transition rates

Interestingly, the term ω_τ(t), which causes divergences of the diffusion coefficient and the frequency of the potential in the FPE (53), does not enter the argument of the exponential in the rates κ_B and κ_K. This means that it just determines the kinetic prefactor, as can be also observed directly from the long-time form ${\partial }_{t}{P}_{1}={\omega }_{\tau }(t){\partial }_{x}\left(x+2{\nu }_{\mathrm{ss}}{\partial }_{x}\right){P}_{1}$ of the FPE (53). In the dynamical regime (ii), the kinetic prefactor in equation (65) cannot be correct, due to the diverging oscillations in the time-dependent frequency ω_τ(t). Nevertheless, the exponential term seems to be reasonable, and thus it is tempting to use in the prefactor of the transition rates simply ω, instead of the problematic ${\omega }_{\tau }(\infty )$. This substitution gives the correct pre-exponential factor of the rate for vanishing delay τ = 0, where ${\omega }_{\tau }(\infty )=\omega$ and $2{D}_{\tau }(\infty )=2D$. We denote the rates with the renormalized prefactor as

$$\begin{eqnarray}&&{\tilde{\kappa }}_{{\rm{x}}}=\omega {\kappa }_{{\rm{x}}}/{\omega }_{\tau }(\infty )\end{eqnarray} \tag{ 66 }$$

with x = M, B or K indicating Markov, Bullerjahn, or Kramers, respectively.

The necessity to change the kinetic prefactor in the rates stems from the fact that, although the absorbing boundary condition we used in equation (63) is correct for Markov dynamics (τ = 0), it can not be precisely valid for the time-delayed dynamics (τ > 0). To see this, it is enough to realize that the delayed system arriving at the boundary at a time t does not feel the energy barrier E_b = γ ω R²/2, but the barrier with energy $\gamma \omega {\left[r(t-\tau )\right]}^{2}/2$.

In figures 9(a)–(c), we compare the various analytical predictions κ_x and ${\tilde{\kappa }}_{{x}}$, x = M, B or K, for the long-time transition rate (lines) with the asymptotic rate κ (symbols) calculated from BD simulations of equation (15) using the inverse first passage time for reaching the absorbing boundary at r = 0. For a broad range of parameters fulfilling ${k}_{B}T\ll V(-R,0)$, we have found that the rate κ can be predicted reasonably well only for values of τ in the dynamical regime (i) (ωτ < 1/e ≈ 0.37) and in the first part of the dynamical regime (ii) (ωτ < π/2 ≈ 1.57). The rate κ is best approximated by the expression ${\tilde{\kappa }}_{{\rm{M}}}$ obtained numerically from equation (63) with ω substituted for ${\omega }_{\tau }(\infty )$ in the operator ${ \mathcal L }$. From the analytical expressions the re-scaled Bullerjahn expression ${\tilde{\kappa }}_{{\rm{B}}}$ (see equations (64) and (66)) works best. However, in the figure we have used parameters leading to a high barrier E_b, and thus Kramer's and Bullerjahn's predictions, κ_K and κ_B, almost coincide. As a consequence, the line for ${\tilde{\kappa }}_{{\rm{K}}}$ in the figure overlaps with ${\tilde{\kappa }}_{{\rm{B}}}$, similarly as the lineκ_K (suppressed in the figure for better readability) overlaps with κ_B.

5.1.5. Delay-dependence of κ

In figure 9(d), we also show the behavior of the rate κ in the parameter regime (iii) inaccessible to the analytical and numerical formulas due to the diverging oscillation. The simulated transition rate in figures 9(a)–(d) is first approximately exponential and thus convex in τ, then its curvature changes to concave and it runs through a maximum and, finally, the rate starts to decrease. The value of τ where the curvature changes sign coincides with the boundary π/2ωτ ≈ 1.57 between the dynamical regimes (ii) and (iii). Interestingly, no qualitative change of κ(τ) is observed at the boundary ωτ = 1/e ≈ 0.37 between the regimes (i) and (ii). It is tempting to attribute, the (approximate) exponential increase of the rate with τ in regimes (i) and (ii) to the increase of the steady-state variance ν_ss, which is given as the ratio of the effective energy barrier $V(-R,\infty )$ and the effective thermal energy $\gamma /2{D}_{\tau }(\infty )$. Although this explanation may work well for small delays, it breaks down for values of τ in the second half of the regime (ii), where the actual rate κ is no longer well approximated by our analytical and numerical predictions. This means that the identification of the parameters γ ω_τ(t) and 2D_τ(t) with the effective potential stiffness and the effective diffusion coefficient, respectively, suggested by the effective Markov equation (53), is reasonable for relatively small values of τ only. In the dynamical regime (iii), the particle undergoes oscillations with an amplitude that increases both with τ and t. The corresponding transition rate, obtained from the BD simulations, thus decreases with the delay τ as a result of oscillations leading away from the transition boundary at r = −R. Note that, in this regime, the stationary transition rate actually does not exist, since the amount of time spent distant from the boundary increases with t (so that the transition rate decreases with t), where t is the duration of the simulation.

5.2. Trimer

Consider the trimer depicted in figure 10 with the three distinguishable particles labeled by the numbers 1–3. We can count the particles either clockwise or anticlockwise and thus two different isomers can form in two dimensions. As in the case of the dimer discussed above, transitions between the two isomers occur with a transition rate depending on the diffusion coefficient, the coupling strength, and the equilibrium spring length.

Zoom In Zoom Out Reset image size

There are several ways how the clockwise isomer may turn into the anticlockwise one and vice versa. For example, the particles 1 and 2 can switch their positions, or the particle 2 can migrate from above the line connecting particles 1 and 3 to below that line, to name a few. The transition rate for hopping between the two isomers is then given as a sum of the transition rates for all the possible realizations of the transition. In order to make an analytical prediction for the transition rate, we choose the coordinate frame in such a way that the x-axis always points from particle 1 to particle 3 (see figure 10). Then all possible transitions between the two isomers boil down to a single event when particle 2 crosses the x-axis. In particular, this also includes the transitions due to exchanging particles 1 and 3. In this case, the direction of the x-axes changes and thus the particle 2 effectively moves to its other side. In the following, we estimate the long-time transition rate $\kappa =\kappa (\infty )$ for the isomer transition in the steady state by means of Kramers' theory and compare it to BD simulations.

In the BD simulations, we have evaluated the rate κ using the angle φ between the abscissas $| 12|$ and $| 13|$ and a neutral region of width ${\rm{\Delta }}x=\sqrt{3/16}R$ as exemplarily shown in figure 10. A neutral state $\left|0\right\rangle$ is introduced to avoid over counting due to fluctuations of φ around 0 and π. It is occupied if the smallest height of the triangle formed by the three particles, is smaller Δ x/2. i.e. either if $| \varphi | \lt {\phi }_{\gt }\equiv \max$ $\left[\arcsin ({\rm{\Delta }}x/2{r}_{12}),\arcsin ({\rm{\Delta }}x/2{r}_{13})\right]$ or $| \varphi | \gt {\phi }_{\lt }\equiv \arccos ({\rm{\Delta }}x/2{r}_{12})+\arccos ({\rm{\Delta }}x/2{r}_{13})$. For φ ∈ [ϕ_> , ϕ_< ] the system is said to be in the clockwise state $\left|1\right\rangle$, while for φ ∈ [−ϕ_> , −ϕ_<] it is in the counter-clockwise state $\left|-1\right\rangle$. To calculate the transition rate, we have counted the number of transitions between the states $\left|\pm 1\right\rangle$ during a specific simulation time window, where the transition occurred if the system underwent the sequence of states $\left|1\right\rangle \to \left|0\right\rangle \to \left|-1\right\rangle$ or $\left|-1\right\rangle \to \left|0\right\rangle \to \left|1\right\rangle$.

The resulting transition rate κ as a function of the time delay τ is depicted in figure 11(a). Therein, one can observe the four dynamical regimes described in section 3. For small delays the rate is approximately exponential, then the curvature of $\kappa (\tau )$ changes from positive to negative, after which the derivative of the rate starts to increase due to the appearance of the fourth dynamical regime where the particle hops between the individual minima of the potential, and for large delays in the unstable regime (iii) the rate drops, while the particle exhibits diverging fluctuations. The transition rate is qualitatively similar to that obtained for the dimer in figure 9. The only difference is that for the dimer we have not observed the fourth dynamical regime, because we have calculated the rate using the first-passage time method, which is insensitive to the potential shape beyond the boundary.

Zoom In Zoom Out Reset image size

For an approximate analytical treatment, we map the situation to a one-dimensional transition problem, to which we apply Kramers' theory. Specifically, we focus on the distance r_b of the particle 2 from the x-axis and construct an appropriate effective energy barrier ${E}_{{\rm{b}},\mathrm{eff}}$ and diffusion coefficient D_eff. In the steady state, the three particles are most likely found at the vertices of an equilateral triangle. Using the notation of section 2.3, we express the vector ${{\bf{r}}}_{{\bf{b}}}$ from the particle 2 to the center of the abscissa $| 13|$ as r_b = r₁₂ + r₃₁/2. Hence the Gaussian white noise ${{\boldsymbol{\eta }}}_{{\rm{b}}}(t)$ corresponding to the coordinate r_b is given by

$$\begin{eqnarray}&&\sqrt{4D}\left({{\boldsymbol{\eta }}}_{1}(t)-{{\boldsymbol{\eta }}}_{2}(t)\right)+\displaystyle \frac{1}{2}\sqrt{4D}\left({{\boldsymbol{\eta }}}_{3}(t)-{{\boldsymbol{\eta }}}_{1}(t)\right)\equiv \sqrt{2\left(3D\right)}{{\boldsymbol{\eta }}}_{{\rm{b}}}(t),\end{eqnarray} \tag{ 67 }$$

where we have used the relations (4). Based on these considerations, we approximate the effective diffusion coefficient by D_eff = 3D.

For the corresponding effective energy barrier, we make the ansatz ${E}_{{\rm{b}},\mathrm{eff}}=C(\tau )\gamma \omega {R}^{2}/12$. Namely, we first consider the minimum value E_b = kR²/6 = γωR²/12 of the energy difference between a configuration with all the three particles aligned (V = kR²/6) and the configuration where the particles form an equilateral triangle with side length R (V = 0), where V is given by equation (2). The (possibly delay-dependent) unknown dimensionless prefactor C(τ) accounts for the time delay and all additional ways the particle 2 may take in the multidimensional energy profile in order to pass from $\left|1\right\rangle$ to $\left|-1\right\rangle$. The resulting Kramers' rate for the transition from $\left|1\right\rangle$ to $\left|-1\right\rangle$ thus reads

$$\begin{eqnarray}&&{\kappa }_{\mathrm{TST}}=a(\tau )\exp \left\{-\displaystyle \frac{C(\tau )\omega {R}^{2}}{36D}\right\},\end{eqnarray} \tag{ 68 }$$

where a(τ) is a further unknown prefactor that may depend on all parameters of the model (see, for example, equation (65) for Kramers' rate in a cusped potential).

In order to test the formula (68), we have fitted it to the transition rate κ obtained from BD simulations as a function of the diffusion coefficient D. The results for different values of the time delay are shown in figure 11(b) and the corresponding values of the coefficients C(τ) and a(τ), obtained from the fits, are given in table 1. The presented results prove that the transition rate exhibits exponential increase with the diffusion coefficient for a relatively broad range of values of the time delay, and thus the D-dependence of κ can be relatively well described by the Kramers-type ansatz (65).

Table 1. The phenomenological parameters C(τ) and a(τ) for five values of the time delay τ corresponding to the three different dynamical regimes of the trimer isomerization. The presented values were obtained by fitting the formula (68) to the BD data shown in figure 11(b).

Regime	τ (s)	C(τ) (1)	a(τ) (s−1)
(i) Exponential decay	0	3.52 ± 0.08	1.16 ± 0.09
(i) Exponential decay	0.1	3.45 ± 0.01	1.10 ± 0.01
(i) Exponential decay	0.3	3.65 ± 0.07	1.69 ± 0.10
(ii) Damped oscillations	0.5	3.41 ± 0.07	2.07 ± 0.12
(iii) Exponential divergence	1.0	1.48 ± 0.05	1.23 ± 0.05

Our attempts to include the effects of the delay in the effective diffusion coefficient and energy barrier, analogously to the approach described in section 5.1, where we made the replacements $\omega \to {\omega }_{\tau ,\mathrm{ss}}$ and $D\to {D}_{\tau ,\mathrm{ss}}$, did not lead to a significant change of the value for C(τ). We thus conclude that the deviations of C(τ) from unity are mainly caused by the assumption that the particle will almost in all cases cross the axis at the minimum of the potential energy. Further improvement would thus require a multidimensional transition rate theory, which is clearly beyond the scope of the present paper.

So far, we have considered only the interactions involving a given positive delay time. But how robust are our analytical findings? Do they critically hinge on (possibly artificial) model details and break down upon some minor variation of the model definition? It turns out that most of the presented results can be directly applied also to other models with delay or memory. To see this, note that the central equation (B.1) is equivalent to the GLE

$$\begin{eqnarray}&&\dot{x}=-\omega {\int }_{0}^{t}{\rm{d}}t^{\prime} \phi (t-t^{\prime} )x(t^{\prime} )+\sigma \eta (t)\end{eqnarray} \tag{ 69 }$$

with positive frequency ω > 0 and the memory kernel given by ϕ(t) = δ(t − τ), τ > 0. Deriving a time-local Langevin equation from the GLE (69) with arbitrary ϕ(t) along the lines of appendix B.1, we obtain equation (B.3) with λ(t) being the Green's function for equation (69), i.e. solving (69) with the initial condition λ(0) = 1 and λ(t) = 0, t < 0. Therefore, all our results that do not depend on the specific form of the Green's function can be readily generalized to arbitrary ϕ(t) after substituting the Green's function (A.5) for ϕ(t) = δ(t − τ) by the Green's function corresponding to the chosen ϕ(t). In the rest of this section, we review some paradigmatic examples of memory kernels ϕ(t), to provide readers with a set of examples of the potential generalizations we have in mind.

The simplest generalization of systems with the memory kernel ϕ(t) = δ(t − τ) are systems with multiple different time delays with the memory kernel

$$\begin{eqnarray}&&\phi (t)=\displaystyle \sum _{i=1}^{N}{\omega }_{i}\delta (t-{\tau }_{i}).\end{eqnarray} \tag{ 70 }$$

Properties of these systems are studied in [35].

A slightly unusual but interesting variation that makes sense in the context of active matter employs a negative delay. The individual active agents may react to a future state of their neighborhood which they predict in the basis of its present state. For idealized systems capable of a perfect prediction, the GLE (or, equivalently, the linear SDDE) contains the memory kernel ϕ(t) = δ(t + τ), τ > 0 and it can be solved using the strategy described in appendix A. The resulting Green's function

$$\begin{eqnarray}&&\lambda (t)=\displaystyle \sum _{k=0}^{\infty }\displaystyle \frac{{\left(-\omega \right)}^{k}}{k!}{\left(t+k\tau \right)}^{k}\end{eqnarray} \tag{ 71 }$$

is the so called Bruwier series [75] that is convergent for $| \omega | \lt | e\tau |$. Alternatively, the series (71) can be written in the form [76]

$$\begin{eqnarray}&&\lambda (t)=\displaystyle \frac{{{\rm{e}}}^{{st}}}{1-\tau s},\end{eqnarray} \tag{ 72 }$$

where s is the absolute value of the smallest root of the equation ρ = −ωe^τρ.

More realistic predictive systems might instead only have an imperfect knowledge of the future position and anticipate a position ${x}_{\mathrm{pre}}(t+\tau )\ne x(t+\tau )$. Therefore, we reformulate the deterministic part of the GLE (69) as

$$\begin{eqnarray}&&\dot{x}(t)=-\omega {x}_{\mathrm{pre}}(t+\tau ).\end{eqnarray} \tag{ 73 }$$

One of the reasonable strategies for predicting x_pre is to use the linear extrapolation

$$\begin{eqnarray}&&{x}_{\mathrm{pre}}(t+\tau )=x(t)+\tau \dot{x}(t),\end{eqnarray} \tag{ 74 }$$

which is identical to a small delay expansion of x(t + τ). The equation of motion (73) then assumes the time local form

$$\begin{eqnarray}&&(1+\omega \tau )\dot{x}(t)=-\omega x(t)\end{eqnarray} \tag{ 75 }$$

which is solved by x(t) = x₀λ(t) with the exponentially decaying Green's function

$$\begin{eqnarray}&&\lambda (t)=\exp \left[-\displaystyle \frac{\omega }{1+\omega \tau }(t-{t}_{0})\right].\end{eqnarray} \tag{ 76 }$$

The rescaled frequency ${\omega }_{{\rm{r}}}=\omega /(1+\omega \tau )$ decreases with increasing delay and thus the resulting dynamics in general exhibits slower relaxation and larger fluctuations (variance) than a system with vanishing time delay. Note that for conventional time delays (now corresponding to τ < 0 in equation (73)), the presented first order approximation predicts dynamics with rescaled frequency ω_r increasing with increasing time delay (for ωτ < −1). However, such dynamics would lead to smaller variance ν_ss for non-zero delays than for a vanishing delay, which contradicts our exact result (24), highlighting the limited practical significance of the naive small-delay expansion, as already pointed out in [30].

The most frequently used generic form of the memory kernel is the exponential

$$\begin{eqnarray}&&\phi (t)=b\exp (-{bt})\end{eqnarray} \tag{ 77 }$$

which is obtained, for example, after integrating out the momentum in the Langevin equation for position of an underdamped harmonic oscillator. This memory kernel leads to the corresponding Green's function

$$\begin{eqnarray}&&\lambda (t)=\exp \left(-\displaystyle \frac{{bt}}{2}\right)\left[\cos ({\rm{\Omega }}t)+\displaystyle \frac{b}{2{\rm{\Omega }}}\sin ({\rm{\Omega }}t)\right],\end{eqnarray} \tag{ 78 }$$

where ${\rm{\Omega }}=\sqrt{{\omega }^{2}-{b}^{2}/4}$, which is reminiscent of the Green's function (A.7) for the system with conventional time delay [ϕ(t) = δ(t − τ)]. The difference is that, while the Green's function (78) always decays exponentially with time, the Green's function (A.7) allows also for negative relaxation times and thus an exponential increase with time. Properties of the Green's function (78), as well as those corresponding to a power-law memory, are discussed in more detail in [77].

Inspired by the surging interest in self-organized active matter and, more specifically, the experiments of Khadka et al [11], we considered N Brownian particles interacting via time-delayed harmonic interactions and confined to a plane, as depicted in figure 1 in section 2. The system is described by the set (3) of 2N nonlinear delayed Langevin equations and hence its dynamics is non-Markovian. At long times, the particles form highly symmetric dynamical molecular-like structures, depicted in figure 4(a) in section 3, which become increasingly compact for large N.

We have analyzed small systems of N = 2 (dimer) and N = 3 (trimer) particles analytically finding molecules with nearest-neighbor distance given by the equilibrium spring length R. To this end, we linearized the corresponding Langevin equations around the zero-temperature steady-state configurations, or, equivalently, around the minimum of the potential energy (2). The linearized Langevin equations could be solved analytically, leading to Gaussian stationary probability densities with delay-dependent effective parameters. In the appendices, we provide analytical expressions for mean values, covariance matrix and time-correlation matrix for a multidimensional system of linear delayed Langevin equations. For the dimer and trimer, we have compared our analytical predictions with BD simulations of the complete model (3). We have found good quantitative agreement in the parameter regimes where the system evolves relatively close to its minimum energy configuration, and good qualitative agreement otherwise (see section 2).

Our analytical results for the dimer and trimer imply that these structures are stable only for small enough values of the product kτ, where k denotes the stifness of the potential. More precisely, these systems converge either exponentially or by exponentially-damped oscillations to corresponding steady states, or they exhibit exponentially diverging oscillations. Our analysis of systems with $N\geqslant 2$ by BD simulations, described in section 3, reveals that these dynamical regimes are stable beyond the linearization approximation and for an arbitrary number of particles. Specifically, we have found that the stability actually extends to larger values of kτ than predicted from the linearized equations, the critical value of the product kτ decaying approximately as 1/N. Therefore, larger systems are more unstable than smaller ones, and the dependence of the stability on the particle number almost vanishes after rescaling the potential stiffness as k → k/N. We conjecture that these instabilities are induced by the chosen form of the interaction which has infinite range and diverges with increasing inter-particle distance. In contrast, the model with constant forces, considered in [11], did not lead to unstable behavior.

Interpreting the inter-particle interactions as an action of a feedback control mechanism, we have, in section 4, used our analytical results for the dimer and trimer to evaluate the amount of entropy extracted by the feedback from (or information injected to) the system in order to maintain the non-equilibrium structures. Interestingly enough, the entropy fluxes do not depend on the noise amplitude D and hence they are discontinuous at D = 0, where the steady-state structures are stable without feedback and thus the entropy fluxes vanish.

Assuming the particles to be distinguishable, the steady-state structures (molecules) can form different isomers. Their transition dynamics can provide rich additional insight into the energy landscape underlying the non-equilirbium structure formation. For the dimer and trimer, we have investigated how and when the transitions between the individual isomers can be described by transition state theory. For the dimer, we have applied our analytical results, based on the time-convolutionless transform leading to the time-local FPE (53), to construct several analytical approximations for the transition rate using Kramers' theory [22, 23] and Bullerjahn's theory [24]. We have also calculated the transition rate from the FPE numerically. Finally, we have compared the obtained predictions to results of BD simulations of the full problem. While the FPE gives the exact value of the transition rate for vanishing delay, our results show that the obtained rates agree with the true ones for small and moderate values of the delay only. We conjecture that this is caused by the fact that the classical absorbing boundary used in our numerical and analytical evaluation of the transition rate can not be used for larger values of τ. Concerning the analytical results, the best agreement with the true rates was obtained by the Bullerjahn's formula (64) with effective barrier height and diffusion coefficient taken from the time-local FPE (53) and the prefactor rescaled according to equation (66). In the case of the trimer, we have confirmed by BD simulations that the transition rate increases exponentially with the noise strength D even for longer delays and thus Kramers' or Bullerjahn's type predictions can be used also in this case. We plan to further investigate suitable absorbing boundary conditions for delayed systems to predict (at least numerically) transition rates also for large delays.

Finally in section 6, we have considered the robustness of our analytical results with respect to details of the realization of the delay. We demonstrated that most of the presented equations can be used also for systems with memory kernels different from that for discrete time delays, i.e. ϕ(t) = δ(t − τ). It is enough to substitute the Green's function λ(t) (A.5) corresponding the delayed Langevin equation by the Green's function corresponding to the memory kernel of interest. We reviewed some paradigmatic memory kernels and provided an outlook on the differences and similarities of the corresponding Green's functions. A more detailed study is left for future work.

As a further extension of our work, it would be interesting to consider physically more realistic interactions that vanish at large distances. Furthermore, we plan to investigate the reaction of the studied system to an external perturbation. Of particular interest could be the propagation and decay behavior of a local perturbation through the system, especially in case of large numbers of particles. Last but not least, we aim to investigate the behavior of the studied system under the action of an additional deterministic time-dependent driving and study the corresponding stochastic dynamics and thermodynamics.

We acknowledge funding by Deutsche Forschungsgemeinschaft (DFG) via SPP 1726/1 and support from the German Research Foundation (DFG) and Universität Leipzig within the program of Open Access Publishing. VH gratefully acknowledges support by the Humboldt foundation and by the Czech Science Foundation (project No. 17-06716S). DG acknowledges funding by International Max Planck Research Schools (IMPRS). Furthermore, we thank Thomas Ihle and Sarah Loos for discussion.

In this appendix we solve the multi-dimensional linear delay differential equation (LDDE)

$$\begin{eqnarray}&&\dot{{\bf{x}}}(t)=-\omega {\bf{x}}(t-\tau ),\end{eqnarray} \tag{ A.1 }$$

where ω is a positive semi-definite matrix with real entries and x(t) is a column vector. Laplace transformation of this equation leads to the formula

$$\begin{eqnarray}&&s\tilde{{\bf{x}}}(s)-{{\bf{x}}}_{0}=-\omega {{\rm{e}}}^{-s\tau }\left[\tilde{{\bf{x}}}(s)+{\int }_{-\tau }^{0}{\rm{d}}t\ {{\rm{e}}}^{-{st}}{\bf{x}}(t)\right]\end{eqnarray} \tag{ A.2 }$$

with x₀ ≡ x(0). The solution of this equation for the Laplace transformed variable reads

$$\begin{eqnarray}\begin{array}{rcl}\tilde{{\bf{x}}}(s) & = & {\displaystyle \int }_{0}^{\infty }{\rm{d}}t\ {{\rm{e}}}^{-{st}}{\bf{x}}(t)={\left(s{ \mathcal I }+\omega {{\rm{e}}}^{-s\tau }\right)}^{-1}\left[{{\bf{x}}}_{0}-\omega {{\rm{e}}}^{-s\tau }{\displaystyle \int }_{-\tau }^{0}{\rm{d}}t\ {{\rm{e}}}^{-{st}}{\bf{x}}(t)\right]\\ & = & \displaystyle \sum _{k=0}^{\infty }\displaystyle \frac{{\left(-\omega \right)}^{k}}{{s}^{k+1}}{{\rm{e}}}^{-ks\tau }\left[{{\bf{x}}}_{0}-\omega {{\rm{e}}}^{-s\tau }{\displaystyle \int }_{-\tau }^{0}{\rm{d}}t\ {{\rm{e}}}^{-{st}}{\bf{x}}(t)\right],\end{array}\end{eqnarray} \tag{ A.3 }$$

where ${ \mathcal I }$ denotes the identity matrix. In the last step, we have expanded the inverse matrix using the Neumann series. The inverse Laplace transform of the ratio ${{\rm{e}}}^{-s\tau }/{s}^{k+1}$ is given by ${\left(t-\tau \right)}^{k}\theta (t-\tau )/k!$ [78] and thus the formula (A.3) can be inverted. Finally, the solution of equation (A.1) is given by

$$\begin{eqnarray}&&{\bf{x}}(t)=\lambda (t){{\bf{x}}}_{0}-\omega {\int }_{-\tau }^{0}{\rm{d}}s\ \lambda (t-s-\tau ){\bf{x}}(s),\end{eqnarray} \tag{ A.4 }$$

where

$$\begin{eqnarray}&&\lambda (t)=\displaystyle \sum _{k=0}^{\infty }\displaystyle \frac{{\left(-\omega \right)}^{k}}{k!}{\left(t-k\tau \right)}^{k}\theta (t-k\tau )\end{eqnarray} \tag{ A.5 }$$

is a matrix-valued function which solves equation (A.1) with the initial condition x(t) = 0 for all t < 0 and ${\bf{x}}(0)={ \mathcal I }$. In the present paper, we always assume that the system is initialized at time t at position x₀ with a special history, namely x(t) = 0 for all t < 0, allowing us to simplify equation (A.4) to x(t) = λ(t)x₀.

The only fixed point of the DDE (A.1) is x(t) = 0. In order to investigate its stability, it is useful to present an alternative solution of the LDDE (A.1) using the exponential ansatz ${\bf{x}}(t)\propto \exp \left(-\alpha t\right)$. Inserting this ansatz into equation (A.1) leads to the equation $\alpha =\omega \exp \left(\alpha \tau \right)$ for the matrix α. Except for some notable exceptions [79], the solution of this equation is given by

$$\begin{eqnarray}&&\alpha =-\displaystyle \frac{1}{\tau }W\left(-\tau \omega \right),\end{eqnarray} \tag{ A.6 }$$

where W denotes the matrix valued Lambert W function [80]. The Lambert W function is a multivalued complex function. The long-time behavior of solutions to equation (A.1) and thus also the stability of its fixed point are determined by the branch of W yielding the largest real parts of the eigenvalues of the matrix α. The corresponding values of the Lambert W function strongly depend on the reduced delay τω.

For example, in one dimension, where ω is a positive real number, the branch of the Lambert W function with the largest real part exhibits three qualitatively different regimes as a function of τω leading to three different dynamical regimes of solutions

$$\begin{eqnarray}&&x(t)\propto {\mathfrak{R}}\left[\exp (-\alpha t)\right]=\exp (-t/{t}_{{\rm{R}}})\cos \left({\rm{\Omega }}t\right),\end{eqnarray} \tag{ A.7 }$$

$1/{t}_{{\rm{R}}}={\mathfrak{R}}(\alpha )$, ${\rm{\Omega }}=\left|{\mathfrak{I}}(\alpha )\right|$, to equation (A.1). The boundaries between these regimes can be determined analytically [25, 27, 30]: (i) $0\lt \tau \omega \leqslant 1/e$, where α is real and positive and x(t) decays exponentially to 0 with lifetime t_R; (ii) 1/e < τω < π/2, where α is complex with a positive real part, producing exponentially damped oscillations of x(t) with frequency Ω and lifetime t_R; and (iii) π/2 < τ ω, where α is complex with a negative real part, corresponding to exponentially diverging oscillations of x(t) with frequency Ω. For τ ω = π/2, 1/t_R = 0 and λ(t) oscillates with frequency Ω without any decay.

In the panel (a) of figure A1, we show that for long times the Green's function λ(t) for equation (A.1) is well approximated by the exponential solution (A.7). The above described dynamical regimes are reflected in the behavior of the decay rate 1/t_R and the frequency Ω of oscillations, plotted in the panels (c) and (d), respectively. The panel (b) of the figure shows the steady auto correlation of x(t) calculated in appendix C, which is also well described by the formula (A.7).

Zoom In Zoom Out Reset image size

B.1. One dimension

Let us now solve equation (A.1) with an additional noise term. For simplicity, we present the detailed derivation first in one dimension. We thus want to solve the equation

$$\begin{eqnarray}&&\dot{x}(t)=-\omega x(t-\tau )+\sigma \eta (t),\end{eqnarray} \tag{ B.1 }$$

where η(t) is a white noise fulfilling $\left\langle \eta (t)\right\rangle =0$ and $\left\langle \eta (t)\eta (t^{\prime} )\right\rangle =\delta (t-t^{\prime} )$. Due to the time delay, this is a time-nonlocal (and consequently non-Markovian) Langevin equation, which, however, can be transformed into a time-local Langevin equation with a colored noise via the so-called time-convolutionless transform [35]. The time-local equation can then be used to derive the FPE for the PDFs for x(t). In order to do that, we write the formal solution of equation (B.1) as

$$\begin{eqnarray}&&x(t)={x}_{0}\lambda (t)+\sigma {\int }_{0}^{t}{\rm{d}}s\ \lambda (t-s)\eta (s)+\psi (t),\end{eqnarray} \tag{ B.2 }$$

where λ(t) is given by equation (A.5) and $\psi (t)=-\omega {\int }_{-\tau }^{0}{\rm{d}}s\ \lambda (t-\tau -s)x(s)$ is determined by the initial condition x(t) for t < 0. As we show in the next section, the formula (B.2) can already be used for the calculation of the time-correlation function for x(t). Here, we differentiate the solution (B.2) which yields the time-local Langevin equation

$$\begin{eqnarray}&&\dot{x}(t)=-{\omega }_{\tau }(t)x(t)+b(t)+\sigma \xi (t).\end{eqnarray} \tag{ B.3 }$$

Due to the time-nonlocal nature of equation (B.1), the potential in the time-local equation (B.3) possess the time-dependent stiffness

$$\begin{eqnarray}&&{\omega }_{\tau }(t)\equiv -\displaystyle \frac{\dot{\lambda }(t)}{\lambda (t)},\end{eqnarray} \tag{ B.4 }$$

and the time-varying position of the minimum $-b(t)/{\omega }_{\tau }(t)$, with $b(t)\equiv {\omega }_{\tau }(t)\psi (t)+\dot{\psi }(t)$ vanishing for the special initial condition x(t) = 0 for all t < 0. Furthermore, equation (B.3) includes the Gaussian colored noise $\xi (t)=\lambda (t)\tfrac{{\rm{d}}}{{\rm{d}}t}{\int }_{0}^{t}{\rm{d}}s\tfrac{\lambda (t-s)}{\lambda (t)}\eta (s)$ which satisfies $\langle \xi (t)\rangle =0$ and

$$\begin{eqnarray}&&\langle \xi (t)\xi (t^{\prime} )\rangle =\lambda (t)\lambda (t^{\prime} )\displaystyle \frac{{\rm{d}}}{{\rm{d}}t^{\prime} }{\int }_{0}^{t^{\prime} }{\rm{d}}s\ \displaystyle \frac{{\rm{d}}}{{\rm{d}}t}\displaystyle \frac{\lambda (t-s)\lambda (t^{\prime} -s)}{\lambda (t)\lambda (t^{\prime} )},\ \ t^{\prime} \lt t.\end{eqnarray} \tag{ B.5 }$$

While Markov processes are completely determined by the the transition probability density ${P}_{1}(x,t| {x}_{0},{t}_{0})$ for going from the initial state x₀ at time t₀ to the final state x at time t, non-Markov processes in general require a full hierarchy of joint probability densities. Nevertheless, similarly to the Markovian case, the Gaussian non-Markov process (B.3) is completely determined by the joint probability distribution ${P}_{2}(x,t;x^{\prime} ,t^{\prime} | {x}_{0},0)$ [35].

The FPEs for the one- and two-time probability distributions ${P}_{1}(x,t| {x}_{0},0)=\langle \delta [x-x(t)]\rangle$ and ${P}_{2}(x,t;x^{\prime} ,t^{\prime} | {x}_{0},0)=\langle \delta [x-x(t)]\delta [x^{\prime} -x(t^{\prime} )]\rangle$, where the averages are taken over all realizations of the process $x(t)$ departing from state x₀ at time 0, are found to be

$$\begin{eqnarray}&&\displaystyle \frac{\partial }{\partial t}{P}_{1}=\displaystyle \frac{\partial }{\partial x}\left[{\omega }_{\tau }(t)x-b(t)+2{D}_{\tau }(t)\displaystyle \frac{\partial }{\partial x}\right]{P}_{1},\end{eqnarray} \tag{ B.6 }$$

$$\begin{eqnarray}&&\displaystyle \frac{\partial }{\partial t}{P}_{2}=\displaystyle \frac{\partial }{\partial x}\left[{\omega }_{\tau }(t)x-b(t)+c(t,t^{\prime} )\displaystyle \frac{\partial }{\partial x}+2{D}_{\tau }(t)\displaystyle \frac{\partial }{\partial x}\right]{P}_{2}.\end{eqnarray} \tag{ B.7 }$$

Similarly as the trap stiffness ω_τ(t), also the effective diffusion coefficient, corresponding to the time-local description, is time dependent if τ > 0. It reads

$$\begin{eqnarray}\begin{array}{rcl}2{D}_{\tau }(t) & \equiv & \displaystyle \frac{{\sigma }^{2}\lambda {\left(t\right)}^{2}}{2}\displaystyle \frac{{\rm{d}}}{{\rm{d}}t}{\displaystyle \int }_{0}^{t}{\rm{d}}s\displaystyle \frac{\lambda {\left(s\right)}^{2}}{\lambda {\left(t\right)}^{2}}={\sigma }^{2}\left(\displaystyle \frac{\lambda {\left(t\right)}^{2}}{2}+{\omega }_{\tau }(t){\displaystyle \int }_{0}^{t}{\rm{d}}s\lambda {\left(s\right)}^{2}\right)\\ & = & {\sigma }^{2}\lambda {\left(t\right)}^{2}/2+{\omega }_{\tau }(t)\nu (t)\end{array}\end{eqnarray} \tag{ B.8 }$$

with the variance $\nu (t)={\sigma }^{2}{\int }_{0}^{t}{\rm{d}}s\lambda {\left(s\right)}^{2}$ and $c(t,t^{\prime} )\equiv {\sigma }^{2}\lambda (t)\tfrac{{\rm{d}}}{{\rm{d}}t}{\int }_{0}^{t^{\prime} }{\rm{d}}s\lambda (t-s)\lambda (t^{\prime} -s)/\lambda (t)$.

Because of the oscillatory nature of λ(t) in the dynamical regimes (ii) and (iii), the coefficients ω_τ, b, c and D_τ in the FPEs (B.6) and (B.7) change their signs and they can even diverge. These divergences, however, always mutually balance each other such that the solutions of the FPEs, as given by the equations (19) and (20), are always reasonable [35].

B.2. Higher dimensions

Let us now consider the problem

$$\begin{eqnarray}&&\dot{{\bf{x}}}(t)=-\omega {\bf{x}}(t-\tau )+\sigma {\boldsymbol{\eta }}(t),\end{eqnarray} \tag{ B.9 }$$

with general matrices ω and σ and the vector ${\boldsymbol{\eta }}(t)$ of white noises fulfilling $\left\langle {\boldsymbol{\eta }}(t)\right\rangle =0$ and $\left\langle {\eta }_{i}(t){\eta }_{j}(t^{\prime} )\right\rangle \,={\delta }_{{ij}}\delta (t-t^{\prime} )$. Since this system of Langevin equations is linear, the one- and two-time probability distributions ${P}_{1}({\bf{x}},t| {{\bf{x}}}_{0},0)$ and ${P}_{2}({\bf{x}},t;{\bf{x}}^{\prime} ,t^{\prime} | {{\bf{x}}}_{0},0)$ for ${\bf{x}}(t)$ defined in the preceding section must be Gaussian [63] as in the one-dimensional case and also the corresponding FPEs can be derived along similar lines as in one dimension. Instead of deriving these equations, we now provide a simpler alternative derivation of the properties of the Gaussian distribution

$$\begin{eqnarray}&&{P}_{1}({\bf{x}},t| {{\bf{x}}}_{0},0)=\displaystyle \frac{1}{\sqrt{{\left(2\pi \right)}^{3}\det { \mathcal K }(t)}}\exp \left\{-\displaystyle \frac{1}{2}{\left({\bf{x}}-{\boldsymbol{\mu }}(t)\right)}^{\top }\cdot { \mathcal K }{\left(t\right)}^{-1}\cdot \left({\bf{x}}-{\boldsymbol{\mu }}(t)\right)\right\}\end{eqnarray} \tag{ B.10 }$$

based solely on the formal solution of the Langevin system (B.9)

$$\begin{eqnarray}&&{\bf{x}}(t)=\lambda (t){{\bf{x}}}_{0}+{\int }_{0}^{t}{\rm{d}}s\ \lambda (t-s)\sigma {\boldsymbol{\eta }}(s)\end{eqnarray} \tag{ B.11 }$$

with the initial condition x(t) = 0 for all t < 0 and x(0) = x₀ and with the Green's function λ(t) of the Langevin system given by equation (A.5).

The mean value ${\boldsymbol{\mu }}(t)=\langle {\bf{x}}(t)\rangle$ and the elements ${{ \mathcal K }}_{{ij}}(t)$ of the covariance matrix ${ \mathcal K }(t)=\left\langle {\bf{x}}(t){{\bf{x}}}^{\top }(t)\right\rangle \,-\left\langle {\bf{x}}(t)\right\rangle \left\langle {{\bf{x}}}^{\top }(t)\right\rangle$ defining the PDF (B.10) can be obtained by inserting x_i(t) from equation (B.11) into the definitions and averaging over the noise η(t). The results are

$$\begin{eqnarray}&&{\boldsymbol{\mu }}(t)=\lambda (t){{\bf{x}}}_{0}\end{eqnarray} \tag{ B.12 }$$

and

$$\begin{eqnarray}&&{ \mathcal K }(t)={\int }_{0}^{t}{\rm{d}}s\ \lambda (t-s)\sigma {\sigma }^{\top }{\lambda }^{\top }(t-s).\end{eqnarray} \tag{ B.13 }$$

These formulas can be generalized straightforwardly to arbitrary initial conditions, where x(t) for $t\leqslant 0$ is drawn from some probability distribution $P[{\bf{x}}(t),t\leqslant 0]$. Then the formal solution of the system (B.1) reads

$$\begin{eqnarray}&&{\bf{x}}(t)={\bf{y}}(t)+{\int }_{0}^{t}{\rm{d}}s\ \lambda (t-s)\sigma {\boldsymbol{\eta }}(s),\end{eqnarray} \tag{ B.14 }$$

$$\begin{eqnarray}&&{\bf{y}}(t)=\lambda (t){\bf{x}}(0)-\omega {\int }_{-\tau }^{0}{\rm{d}}s\,\lambda (t-s-\tau ){\bf{x}}(s).\end{eqnarray} \tag{ B.15 }$$

The mean value ${\boldsymbol{\mu }}(t)$ is given by ${\boldsymbol{\mu }}(t)={\left\langle {\bf{y}}(t)\right\rangle }_{{\bf{x}}(t\leqslant 0)}$, and ${ \mathcal K }(t)={\left\langle {\bf{y}}(t){\left[{\bf{y}}(t)\right]}^{\top }\right\rangle }_{{\bf{x}}(t\leqslant 0)}+{\int }_{0}^{t}{\rm{d}}s\ \lambda (t-s)\sigma {\sigma }^{\top }{\lambda }^{\top }(t-s)$. The averages ${\left\langle \bullet \right\rangle }_{{\bf{x}}(t\leqslant 0)}$ above are taken with respect to the PDF $P[{\bf{x}}(t),t\leqslant 0]$. The long-time behavior of the covariance matrix (B.13) is studied in the next section of this appendix.

The coefficients in the Gaussian two-time PDF ${P}_{2}({\bf{x}},t;{\bf{x}}^{\prime} ,t^{\prime} | {{\bf{x}}}_{0},0)$ can be obtained in a similar manner as in appendix B.2. Here, we calculate only the stationary space-time correlation matrix ${ \mathcal C }(t)\,={\mathrm{lim}}_{s\to \infty }\left[\left\langle {\bf{x}}(s+t){{\bf{x}}}^{\top }(s)\right\rangle -\left\langle {\bf{x}}(s+t)\right\rangle \left\langle {{\bf{x}}}^{\top }(s)\right\rangle \right]={\mathrm{lim}}_{s\to \infty }\left\langle {\bf{x}}(s+t){{\bf{x}}}^{\top }(s)\right\rangle$ which exists only if ${\mathrm{lim}}_{s\to \infty }\left\langle {\bf{y}}(s)\right\rangle \,={\mathrm{lim}}_{s\to \infty }\left\langle {\bf{x}}(s)\right\rangle =0$. Its matrix elements can be calculated in an analogous fashion as the elements of ${ \mathcal K }(t)$. The result

$$\begin{eqnarray}&&{ \mathcal C }(t)=\mathop{\mathrm{lim}}\limits_{s\to \infty }{\int }_{0}^{s}{\rm{d}}s^{\prime} \ \lambda (s+t-s^{\prime} )\sigma {\sigma }^{\top }{\lambda }^{\top }(s-s^{\prime} )\end{eqnarray} \tag{ C.1 }$$

can be evaluated numerically. It is possible to rewrite it in a simpler form. Taking the derivative of ${ \mathcal C }(t)$ with respect to t reveals that for t > 0 the correlation matrix obeys the same DDE as the Green's function λ(t), i.e.

$$\begin{eqnarray}&&\dot{{ \mathcal C }}(t)=-\omega { \mathcal C }(t-\tau ),\ \ t\gt 0.\end{eqnarray} \tag{ C.2 }$$

The restriction t > 0 for validity of this equation comes from discontinuity (and thus non-differentiability) of λ(t) at t = 0. The solution to equation (C.2) is given by equation (A.4) and hence the stationary space-time correlation matrix can be written as

$$\begin{eqnarray}&&{ \mathcal C }(t)=\lambda (t){{ \mathcal C }}_{0}-\omega {\int }_{-\tau }^{0}{\rm{d}}s\ \lambda (t-s-\tau ){ \mathcal C }(s),\end{eqnarray} \tag{ C.3 }$$

where ${{ \mathcal C }}_{0}={ \mathcal C }(0){\mathrm{lim}}_{t\to \infty }{ \mathcal K }(t)$ is given by the long time limit of the covariance matrix. The stationary correlation matrix is thus solely determined by the unknown initial condition ${ \mathcal C }(t)$, t ∈ [−τ, 0]. Fortunately, this initial condition can be calculated using the approach of Frank et al [31] who calculated the time-correlation function ${ \mathcal C }(t)$ for t ∈ [0, τ] in one-dimension (see also [32]). They utilized the symmetry ${ \mathcal C }(t)={ \mathcal C }(-t)$ following from equation (C.1) to rewrite the DDE (C.2) as $\dot{{ \mathcal C }}(t)=-\omega { \mathcal C }(\tau -t)$ for t ∈ (0, τ). Taking the derivative of this equation and using equation (C.2) yields the second order ordinary differential equation

$$\begin{eqnarray}&&\ddot{{ \mathcal C }}(t)=-{\omega }^{2}{ \mathcal C }(t),\ \ t\in (0,\tau )\end{eqnarray} \tag{ C.4 }$$

for the initial condition ${ \mathcal C }(t)={ \mathcal C }(-t)$, t ∈ [−τ, 0]. The solution of this equation reads ${ \mathcal C }(t)={{ \mathcal C }}_{0}\cos (\omega t)\,+{\dot{{ \mathcal C }}}_{0}{\omega }^{-1}\sin (\omega t)$, where we still need to determine the unknown coefficients ${{ \mathcal C }}_{0}={ \mathcal C }(0)$ and ${\dot{{ \mathcal C }}}_{0}={\mathrm{lim}}_{t\to 0+}\dot{{ \mathcal C }}(t)$. To this end, we need to evaluate independently ${ \mathcal C }(t)$ and/or $\dot{{ \mathcal C }}(t)$ for two times t∈(0, τ). Specifically, we show at the end of this appendix that ${ \mathcal C }(\tau )=0.5{\omega }^{-1}\sigma {\sigma }^{\top }$ and ${\dot{{ \mathcal C }}}_{0}=-0.5\sigma {\sigma }^{\top }$. Using these results, the final expression for the correlation matrix for t ∈ [−τ, τ] reads

$$\begin{eqnarray}&&{ \mathcal C }(t)={{ \mathcal C }}_{0}\cos (\omega t)-0.5\sigma {\sigma }^{\top }{\omega }^{-1}\sin (\omega | t| )\end{eqnarray} \tag{ C.5 }$$

with the initial value

$$\begin{eqnarray}&&{ \mathcal C }(0)={{ \mathcal C }}_{0}=\mathop{\mathrm{lim}}\limits_{s\to \infty }{ \mathcal K }(s)=\displaystyle \frac{1}{2}\left[{\omega }^{-1}\sigma {\sigma }^{\top }+\sigma {\sigma }^{\top }{\omega }^{-1}\sin (\omega \tau )\right]{\cos }^{-1}(\omega \tau )\end{eqnarray} \tag{ C.6 }$$

given by the stationary value of the covariance matrix. The whole time-dependence of ${ \mathcal C }(t)$ for $t\geqslant 0$ is thus described by the formulas (C.3), (C.5), and (C.6). The correlation matrix for negative times then follows from the symmetry ${ \mathcal C }(t)={ \mathcal C }(-t)$. Finally, let us note than in the one-dimensional case, where ω and σ stand for real numbers and, more generally, if the matrices ω⁻¹ and $\sigma {\sigma }^{\top }$ commute, we can rewrite equation (C.6) as

$$\begin{eqnarray}&&{{ \mathcal C }}_{0}=\displaystyle \frac{\sigma {\sigma }^{\top }}{2\omega }\left[{ \mathcal I }+\sin (\omega \tau )\right]{\cos }^{-1}(\omega \tau ),\end{eqnarray} \tag{ C.7 }$$

where ${ \mathcal I }$ denotes the identity matrix. An example of the time-correlation function for a one dimensional system is depicted in figure A1(b) in appendix A.

The expression for C(τ) can be obtained by multiplying equation (B.9) by ${{\bf{x}}}^{\top }(s)$, averaging the result over the noise, using the assumed stationarity of the process implying the formula $\dot{{ \mathcal C }}(0)={\rm{d}}\left[{\mathrm{lim}}_{s\to \infty }\left\langle {\bf{x}}(s){{\bf{x}}}^{\top }(s)\right\rangle \right]/{\rm{d}}s=0$, and applying the symmetry ${ \mathcal C }(t)={ \mathcal C }(-t)$. The result is ${ \mathcal C }(\tau )={\omega }^{-1}\sigma \left\langle {\boldsymbol{\eta }}(s){{\bf{x}}}^{\top }(s)\right\rangle =0.5{\omega }^{-1}\sigma {\sigma }^{\top }$, where the last equality comes after inserting the formal solution (B.11) for x(t) into the average, using the covariance of the noise, and noticing that in the resulting integral we integrate over half of the emerging δ-function only. The expression for ${\dot{{ \mathcal C }}}_{0}$ then comes simply from equation (C.2), which is invalid for t = 0, but can be used for t arbitrarily close to 0 from the right, and using the result for ${ \mathcal C }(\tau )$.