Brownian Molecules Formed by Delayed Harmonic Interactions

TL;DR

This paper investigates how delayed harmonic interactions among Brownian particles lead to complex, molecule-like structures, analyzing their stability, entropy, and transition dynamics through analytical and simulation methods.

Contribution

It provides the first analytical solutions for small systems and develops a new transition-state theory for understanding stability and transitions in delayed stochastic systems.

Findings

Particles form unstable molecule-like structures with increasing delay and interaction strength.

Analytical solutions for N ≤ 3 particles are derived.

A new approximate transition-state theory characterizes structural transitions.

Abstract

A time-delayed response of individual living organisms to information exchanged within flocks or swarms leads to the emergence of complex collective behaviors. A recent experimental setup by (Khadka et al 2018 Nat. Commun. 9 3864), employing synthetic microswimmers, allows to emulate and study such behavior in a controlled way, in the lab. Motivated by these experiments, we study a system of N Brownian particles interacting via a retarded harmonic interaction. For $N \leq 3$ , we characterize its collective behavior analytically, by solving the pertinent stochastic delay-differential equations, and for $N>3$ by Brownian dynamics simulations. The particles form molecule-like non-equilibrium structures which become unstable with increasing number of particles, delay time, and interaction strength. We evaluate the entropy and information fluxes maintaining these structures and, to quantitatively characterize their stability, develop an approximate time-dependent transition-state theory to characterize transitions between different isomers of the molecules. For completeness, we include a comprehensive discussion of the analytical solution procedure for systems of linear stochastic delay differential equations in finite dimension, and new results for covariance and time-correlation matrices.

Tables

Table 1. Table 1: The phenomenological parameters C ​ ( τ ) 𝐶 𝜏 C(\tau) and a ​ ( τ ) 𝑎 𝜏 a(\tau) for five values of the time delay τ 𝜏 \tau 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 Fig.  11 b).

Regime	$\tau$ [s]	$C\left(\tau \right)$ [1]	$a\left(\tau \right)$ [s-1]
(i) exponential decay	0	3.52 $\pm$ 0.08	1.16 $\pm$ 0.09
(i) exponential decay	0.1	3.45 $\pm$ 0.01	1.10 $\pm$ 0.01
(i) exponential decay	0.3	3.65 $\pm$ 0.07	1.69 $\pm$ 0.10
(ii) damped oscillations	0.5	3.41 $\pm$ 0.07	2.07 $\pm$ 0.12
(iii) exponential divergence	1.0	1.48 $\pm$ 0.05	1.23 $\pm$ 0.05

Equations

$$V_2[r_{ij}(t-\tau )]=\frac{k}{2}{\left[r_{ij}(t-\tau )-R\right]}^2,$$

$$V=\frac{1}{2}\sum _{(i,j)}{V}_2(r_{ij}),$$

$${\dot{\mathbf{r}}}_i(t)=-\frac{k}{\gamma }\sum _{j\ne i}\left[r_{ij}(t-\tau )-R\right]{\mathbf{e}}_{ij}(t-\tau )+\sqrt{2D_0}{\mathbf{\eta }}_i(t),i=1,\dots ,N.$$

$$⟨{\mathbf{\eta }}_i(t)⟩=0,⟨{\eta }_i^{\alpha }(t_1){\eta }_j^{\beta }(t_2)⟩={\delta }_{ij}{\delta }_{\alpha \beta }{\delta }_{\mathrm{D}}(t_1-t_2).$$

$${\dot{\mathbf{r}}}_i(t)=-\frac{k}{\gamma }\frac{{\sum }_{j\ne i}\mathrm{sign}\left[r_{ij}(t-\tau )-R\right]{\mathbf{e}}_{ij}(t-\tau )}{\left|{\sum }_{j\ne i}\mathrm{sign}\left[r_{ij}(t-\tau )-R\right]{\mathbf{e}}_{ij}(t-\tau )\right|}+\sqrt{2D_0}{\mathbf{\eta }}_i(t),$$

$${\dot{\mathbf{r}}}_{\mathrm{c}}(t)=\sqrt{2D_c}{\mathbf{\eta }}_{\mathrm{c}}(t),$$

$$P_N^{\mathrm{c}}(\mathbf{r},t)=\sqrt{\frac{N}{4\pi Dt}}\exp \left[-N\frac{(\mathbf{r}-{\mathbf{r}}_0{)}^2}{4Dt}\right].$$

$${\dot{\mathbf{r}}}_1(t)$$

$${\dot{\mathbf{r}}}_2(t)$$

$$\dot{\mathbf{r}}(t)=-\omega \left[r(t-\tau )-R\right]\mathbf{e}(t-\tau )+\sqrt{4D}{\mathbf{\eta }}_{\mathrm{r}}(t)$$

$$\omega \equiv 2k/\gamma$$

$$\dot{r}(t)=$$

$$\dot{\phi }(t)=$$

$$\dot{\phi }(t)=\sqrt{\frac{4D}{R^2}}{\eta }_{\mathrm{r}}^{\phi }(t).$$

$$\dot{r}(t)=-{\omega }_{\tau }\left[r(t-\tau )-R\right]+\sqrt{4D}{\eta }_{\mathrm{r}}^r(t),$$

$$x(t)=r(t)-R,$$

$$x(t)=x_0\lambda (t)+\sqrt{4D}{\int }_0^t\mathrm{d}s\lambda (t-s){\eta }_{\mathrm{r}}^r(s),$$

$$\lambda (t)=\sum _{k=0}^{\infty }\frac{{\left(-{\omega }_{\tau }\right)}^k}{k!}{\left(t-k\tau \right)}^k\theta (t-k\tau )$$

$$\displaystyle P_{1}(x,t|x_{0},0)=\frac{1}{\sqrt{2\pi\nu(t)}}\exp\bigg{\{}-\frac{1}{2}\bigg{(}\frac{x-\mu(t)}{\sqrt{\nu(t)}}\bigg{)}^{2}\bigg{\}},$$

$$P_2(x,t;x^{\prime },t^{\prime }|x_0,0)=\frac{1}{\sqrt{4{\pi }^2\nu (t)\nu (t^{\prime })[1-w^2(t^{\prime },t)]}}$$

$$\displaystyle\qquad\times\exp\bigg{\{}\frac{1}{2[1-w^{2}(t^{\prime},t)]}\bigg{[}\bigg{(}\frac{x-\mu(t)}{\sqrt{\nu(t)}}\bigg{)}^{2}+\bigg{(}\frac{x^{\prime}-\mu(t^{\prime})}{\sqrt{\nu(t^{\prime})}}\bigg{)}^{2}$$

$$\displaystyle\qquad-2w(t,t^{\prime})\bigg{(}\frac{x-\mu(t)}{\sqrt{\nu(t)}}\bigg{)}\bigg{(}\frac{x^{\prime}-\mu(t^{\prime})}{\sqrt{\nu(t^{\prime})}}\bigg{)}\bigg{]}\bigg{\}},$$

$$\mu (t)$$

$$\nu (t)$$

$$w(t,t^{\prime })$$

$${\nu }_{\mathrm{ss}}=\underset{t\to \infty }{\lim }\nu (t)=\frac{2D}{{\omega }_{\tau }}\frac{1+\sin ({\omega }_{\tau }\tau )}{\cos ({\omega }_{\tau }\tau )}.$$

$${\nu }_{\mathrm{ss}}\ll R^2$$

$${\dot{\mathbf{r}}}_1(t)$$

$${\dot{\mathbf{r}}}_2(t)$$

$${\dot{\mathbf{r}}}_3(t)$$

$$\begin{gathered} {\dot{\mathbf{r}}}_{12}(t)=-\omega \{\left[r_{12}(t-\tau )-R\right]{\mathbf{e}}_{12}(t-\tau )+\frac{1}{2}\left[r_{13}(t-\tau )-R\right]{\mathbf{e}}_{13}(t-\tau ) \\ +\frac{1}{2}\left[r_{32}(t-\tau )-R\right]{\mathbf{e}}_{32}(t-\tau )\}+\sqrt{2D}\left[{\mathbf{\eta }}_1(t)-{\mathbf{\eta }}_2(t)\right], \end{gathered}$$

$$\displaystyle\dot{x}_{\alpha}(t)=-\omega\bigg{[}x_{\alpha}(t-\tau)+\frac{1}{4}x_{\alpha+\mathrm{I}}(t-\tau)+\frac{1}{4}x_{\alpha+\mathrm{II}}(t-\tau)\bigg{]}$$

$$\mathcal{A}=\left(\begin{array}{ccc}1& -1& 0\\ 1& 0& -1\\ 0& -1& 1\end{array}\right),$$

$$\dot{\mathbf{x}}(t)=-\omega M\mathbf{x}(t-\tau )+\sqrt{2D}\mathbf{\xi }(t),$$

$$M=\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}_{\mathrm{I}}^1& -e_{\mathrm{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}_{\mathrm{I}}^2& -e_{\mathrm{I}}^2& 0\\ e_{\mathrm{II}}^2& 0& -e_{\mathrm{II}}^2\\ 0& -e_{\mathrm{III}}^2& e_{\mathrm{III}}^2\end{array}\right),$$

Full text

††: New J. Phys.

• August 2017

Keywords: Stochastic Delay Differential Equation, Active Matter, Feedback Driving, Non-Markovian Dynamics, Fokker-Planck Equation, Structure Formation, Transition-State Theory

1 Introduction

1.1 Feedback systems

From the synchronized response of a flock of starlings [1] avoiding an attack of a predator to the formation of colonies of living bacteria [2, 3], the surging field of active matter provides a wide range of fascinating phenomena. Its ultimate aim is to develop a microscopic understanding of the behavior of large numbers of interacting, active and energy consuming “agents” [4, 5], with a focus on emergent collective behavior [6]. Most of the quantitative models, such as the Vicsek model [7], neglect the finite speed of information transmission between the individual particles. However, recent studies [8, 9, 10, 11] have shown that a time delay in the interaction may significantly affect the system dynamics. Moreover, experimental realizations mimicking natural interacting systems require implementing the non-physical interactions, such as a reaction of a bird to its environment, via a feedback loop [10, 11, 12, 13]. Finite processing of the information in the feedback loop then inevitably introduces time delay into the system dynamics.

Current (mainly optical) micromanipulation techniques allow to realize such feedback systems on microscale [11, 13, 14, 15, 16, 17, 18, 19] . Many [13, 18, 19] of these techniques are based on spherical Janus particles [20, 21] with hemispheres coated with different materials in order to excite surface flows to propel them actively upon illumination or in presence of other energy sources (e.g. chemical fuel added to the solvent). In order to steer these particles, one usually has to wait until the rotational diffusion reorients them towards the desired location. This issue was resolved by the setup introduced by Khadka et al. [11] based on Brownian particles, symmetrically decorated by gold nanoparticles, that thermophoretically self-propel in the direction determined by the position of the laser focus on their circumference. In the feedback experiment, the particles are tracked with a camera with finite exposure time and the position of the heating laser is determined by positions of the particles in the previous frame. The setup allows to create arbitrary time-delayed interactions in the many-body system. In Ref. [11], an interaction leading to constant absolute values of velocities of the individual particles was considered.

In the present paper, we theoretically analyze a system similar to that considered by Khadka et al. [11], but with harmonic interactions between the individual particles. Similarly to the case of Ref. [11], the two-dimensional $N$-particle system is described by a set of $2N$ coupled non-linear stochastic delay differential equations. For small enough values of the delay, highly symmetric non-equilibrium molecular-like structures form after a transient period, which fluctuate due to thermal noise. The resulting structures strongly differ from the molecules created with the constant-velocity protocol studied in Ref. [11], which oscillated, even for vanishing noise amplitude, due to the nonzero delay. Another difference between the two realizations is that our setup leads, for large delays, to oscillations with amplitudes exponentially increasing in time, while, in the setup of Khadka et al., they are always bounded. The specific form of the interaction considered in our setup moreover allows us to linearize the underlying set of stochastic delay differential equations and to study many aspects of the model behavior analytically. For dimer ($N=2$) and trimer ($N=3$), we use the linearized model to calculate properties of the resulting Gaussian probability distributions for the bond length, namely its mean values, covariance matrix, and time-correlation matrix. We verify the validity of these results by Brownian dynamics (BD) simulations of the complete model. Moreover, we use the BD simulations to show that the behavior of larger molecules ($N>3$) is qualitatively the same as that of the dimer and the trimer, with the difference that the critical value of the delay, beyond which the molecules become unstable, decreases inversely in the particle number $N$. If we would scale the interaction strength by the particle number, as it is common in toy models of many-body systems, the critical value of the delay would thus be constant. If we label the individual particles, we can distinguish between several isomers of the respective molecules according to their ordering. In the course of time, the noise induces jumps of a given molecule between the individual isomers. We utilize our analytical results for the dimer and for the trimer to evaluate the corresponding transition rates using Kramers’ theory [22, 23] and a more recent theory by Bullerjahn et al. [24]. We compare the results with the rates calculated from our BD simulations and identify a useful formula for the transition rate that provides good predictions for small and moderate values of the delay.

1.2 Stochastic delay differential equations

In general, delay differential equations (DDE’s) [25, 26] may generate rich dynamics [27]. Their solutions may converge to fixed points or limit cycles, behave chaotically, and exhibit multistability [28]. For systems affected by noise, the DDEs are generalized to stochastic delay differential equations (SDDE) [29], which exhibit non-Markovian dynamics. Due to delay-induced temporal correlations, the corresponding Fokker-Planck equation (FPE) cannot be written in a closed form [30, 31, 32]. Instead, one obtains an infinite hierarchy of coupled FPEs for the $n$-time joint probability densities for which generally no finite closure is known.

For non-linear systems, there are three established approximate approaches how to tackle the infinite hierarchy: (i) the so called small delay approximation [30], which employs a Taylor expansion in the delay to make the equations time local; (ii) also, if the delayed terms in the SDDE are small so that the system dynamics is almost Markovian, a perturbation theory can be applied, leading to closed FPEs for the individual joint probability densities [33]; (iii) a closed equation for the 1-time probability density valid in the steady state can be obtained by linearization of all equations of the FPE hierarchy except for the first one [32].

So far, the only exactly solved problem is a one-dimensional linear stochastic delay equation with Gaussian white noise, whose $n$-time probability densities are given by multivariate Gaussians. Its stationary solution and the conditions for its existence were discussed in Refs. [34, 30, 31]. Recently, a full time-dependent solution for 1- and 2-time probability densities was found by Giuggioli et al. [35]. Employing the so-called time-convolutionless transform introduced in the 1970s [36, 37, 38, 39, 40], these authors transformed the non-Markovian linear delayed Langevin equation (LE) into a time-local form. Afterwards, they utilized this result in a derivation of analytically solvable time-local FPEs for 1- and 2-time probability densities.

Even though an analytical treatment is thus rather complicated, there is a great interest in understanding DDE’s and SDDE’s, due to their broad range of applications. Prominent examples are found in population dynamics [41, 42], where the delay results from maturation times, economics [43, 44, 45, 46] when the limited reaction times of the market participants matters, or engineering [47]. In biology, finite transition times can play a significant role in physiological systems [48, 49, 50] and neural networks [51, 52, 53, 28]. Recently [54, 55, 56], first efforts were also made to incorporate a time delay into the language of stochastic thermodynamics [57, 58] in order to evaluate energy and entropy fluxes in time-delayed stochastic system.

1.3 Outline

The rest of the paper is structured as follows. In Sec. 2 we first introduce the general model and formulate the underlying equations of motion in terms of SDDEs. After appropriate linearization, we study them analytically, considering both transient and stationary properties of the probability distributions for “bond” lengths in “molecules” self-assembling from two and three particles. Larger systems are studied in Sec. 3. In Sec. 4, we apply the obtained results for evaluation of the entropy outflux (or information influx) from the system due to the feedback maintaining the non-equilibrium stationary structures. In order to obtain a more quantitative characterization of the stability of the non-equilibrium molecules, we utilize our analytical description of the dimer and trimer for analytical and numerical investigation of the isomer transitions and back up the results by BD simulations in Sec. 5. In order to assess the robustness of our findings, Sec. 6 addresses the role of the functional form of the memory kernel considering negative delays and exponential memories. We summarize our findings and conclude in Sec. 7. Most of the technical details are given in A – C. In A, we review the known results concerning the solution of systems of LDDEs. In B, we show how to generalize these results for linear SDDEs. Finally, we apply the obtained results in C for the calculation of the time-correlation matrix and the covariance matrix for systems of linear SDDEs.

2 Stochastic Dynamics

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

[TABLE]

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

Altogether, the particles diffuse in the composity potential

[TABLE]

where the summation runs over all pairs $(i,j)$, so that the $i$th particle is driven by the force $\mathbf{F}_{i}=-\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-\tau$, here $x_{i}$ and $y_{i}$ denote the Cartesian coordinates of the position vector $\mathbf{r}_{i}$ in the past. In effect, the $N$-particle system therefore obeys the set of nonlinear delayed Langevin equations

[TABLE]

The unit vector $\mathbf{e}_{ij}=\mathbf{r}_{ij}/r_{ij}$ points from particle $j$ to particle $i$ and the diffusion coefficient $D_{0}=(\beta\gamma)^{-1}$ is related to the inverse temperature $\beta=1/k_{\rm B}T$ and the friction coefficient $\gamma$ via the Einstein relation ($k_{\rm B}$ denotes the Boltzmann constant). The vectors $\bm{\eta}_{i}$ comprise independent Gaussian white noises satisfying the relations

[TABLE]

The numbers $\alpha$ and $\beta$ label the components of the vector $\bm{\eta}_{i}$, while $i$ refers to the specific particle. Note that the noise and the friction in Eq. (3) are related by the fluctuation-dissipation theorem [59] for a vanishing delay $\tau=0$ only, and that the system is always out of thermodynamic equilibrium for $\tau>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 Fig. 4 a) in Sec. 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 $\beta$, 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 Eq. (4) that leads to non-trivial non-equilibrium steady states. Consequently, the Boltzmann distribution as stated above has to be modified. For the simplest case of a dimer one finds an (approximate) Gaussian distribution as discussed later while for larger molecules the description is 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

[TABLE]

$i=1,\dots,N$, with ${\rm sign}(x)$ denoting the signum function, was discussed earlier in [11]. The main difference from our setting (3) is that, in Eq. (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 Eq. (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 Ref. [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 Sec. 3.

2.1 Center of Mass

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

[TABLE]

where $\bm{\eta}_{\mathrm{c}}\equiv\sum_{i=1}^{N}\bm{\eta}_{i}/\sqrt{2}$ denotes Gaussian white noise satisfying Eq. (4) (with the labels $i,j$ replaced by ${\rm 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 $\mathbf{r}_{0}$, the probability density function (PDF) for $\mathbf{r}_{\rm c}$ reads

[TABLE]

2.2 Dimer

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

[TABLE]

where we have used the abbreviations $\mathbf{e}(t)\equiv\mathbf{e}_{12}(t)$ and $r(t)\equiv r_{12}(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 $\mathbf{r}=\mathbf{r}_{12}=\mathbf{r}_{1}-\mathbf{r}_{2}$ which obeys the equation of motion

[TABLE]

with frequency

[TABLE]

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

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

[TABLE]

where $\varphi(t,t-\tau)=\varphi(t)-\varphi(t-\tau)$ denotes the change of orientation of the vector $\mathbf{e}(t-\tau)$ during time $\tau$. Above, we used the formulas $\mathbf{r}\equiv r\mathbf{e}$, $\dot{\mathbf{r}}\equiv\dot{r}\mathbf{e}+r\dot{\varphi}\mathbf{e}_{\varphi}$ and $\bm{\eta}_{\rm r}\equiv\eta_{\rm r}^{r}\mathbf{e}+\eta_{\rm r}^{\varphi}\mathbf{e}_{\varphi}$.

From symmetry considerations, it follows that the stationary PDF for the orientation must be constant in $\varphi$. 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 Eq. (12) is strong enough such that the terms proportional to $[r(t-\tau)-R]/R$ can safely be neglected independently of the noise strength, the formula (13) assumes the form 111In Eq. (12), we set $r=[(r-R)/R+1]R$, expand it in $(r-R)/R$, and neglect all terms proportional to $(r-R)/R$.

[TABLE]

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

[TABLE]

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

Equation (15) is a linear stochastic delay differential equation which can be solved analytically for $r\in(-\infty,\infty)$. In our setting, $r\geq 0$ and thus we should solve Eq. (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 Eq. (15) on the full real axis should approximate well the desired solution on the positive half-line. The solution of Eq. (15) for $r\in(-\infty,\infty)$ and $t\geq 0$ in terms of deviations of the bond length from its equilibrium length (which we call as shifted bond length),

[TABLE]

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

[TABLE]

where the dimensionless Green’s function

[TABLE]

is obtained by solving Eq. (15) with $D=0$ for the specfic initial conditions $\lambda(t)=0$ for $t<0$ and $\lambda(0)=1$. The symbol $\theta(x)$ in Eq. (18) stands for the Heaviside step function. For an arbitrary initial condition $x(t)$, $t<0$ and $x(0)=x_{0}$, the expression $x_{0}\lambda(t)$ in Eq. (17) must be substituted by $x_{0}\lambda(t)-\omega_{\tau}\int_{-\tau}^{0}\mathrm{d}s\ \lambda(t-\tau-s)x(s)$. Based on the value of the reduced delay $\omega_{\tau}\tau$ which is a dimensionless measure for the relevance of the delay relative to the natural relaxation time, the Green’s function $\lambda(t)$ in Eq. (17) exhibits three different dynamical regimes discussed in detail in A, in Fig. 12 and also below: (i) monotonic exponential decay to zero for short delays $0\leq\omega_{\tau}\tau\leq 1/e$, (ii) oscillatory exponential decay to zero for intermediate (‘resonant’) delays $1/e\leq\omega_{\tau}\tau\leq\pi/2$, and (iii) oscillatory exponential divergence for long delays $\omega_{\tau}\tau>\pi/2$.

The stochastic process $x(t)$ in Eq. (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 $\delta(x)$ for $t<0$ and $\delta(x-x_{0})$ at $t=0$ read

[TABLE]

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

[TABLE]

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_{0}$ at time [math] is found at time $t$ somewhere in the interval $(x,x+{\rm d}x)$. 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+{\rm d}x)$ under the condition that it has started at time $t=0$ at $x_{0}$. The one-time PDF $P_{1}$ possesses the same structure as the corresponding PDF for $\tau=0$ [63]. The non-Markov character of the process (15) with nonzero delay manifests itself in the fact that the two-time PDF $P_{2}$ cannot be constructed from the one-time PDF $P_{1}$, while this is always possible for a Markov process.

The FPEs corresponding to the PDFs (19) and (20) are given by Eqs. (91) and (92) in B, respectively. Interestingly enough, the diffusion and drift terms in the FPEs are given by the natural scales $2D$ and $\omega_{\tau}$ only in the limit $\tau\to 0$. Moreover, both coefficients acquire a time dependence, determined by the function $\lambda(t)$. Specifically, the diffusion and drift coefficients in Eq. (91) for $P_{1}$ 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 222The case of $\lambda(t)=0$, where these coefficients diverge, is discussed in more detail in Sec. 5.1.. While the drift coefficient in Eq. (92) for $P_{2}$ is also given by $\omega_{\tau}(t)$, the diffusion coefficient reads $2D_{\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_{2}$ can be constructed from $P_{1}$ in the standard way for Markov processes only for $\tau=0$, where both diffusion coefficients coincide.

For $\tau=0$ the PDF $P_{1}$ 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 $\mu(t)=x_{0}\lambda(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)\neq\omega$ and $D_{\tau}(\infty)\neq 2D$, see Fig. 8 in Sec. 5.1. This means that, due to the delay, the diffusion coefficient and frequency in the Gaussian Boltzmann distribution for the approximate model have to be replaced by their effective counterparts. The corresponding 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 $\tau$ 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<x^{2}(t)\right>-\left<x(t)\right>^{2}}$ converges to the stationary value [34, 30, 31]

[TABLE]

The derivation of this formula is given in C, where we also derive an analytical expression for the stationary time correlation function $C(t)=\lim_{s\to\infty}\left<x(s)x(s+t)\right>$. Note that the variance (24) diverges upon entering the unstable regime (iii) for $\omega_{\tau}\tau\rightarrow\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 Eq. (15) from Eqs. (12) and (13) is reasonable because the PDF for $r$ is relatively sharply peaked around the mean bond length $R$. As we already noted, Eq. (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 $\varphi(t,t-\tau)$ in Eqs. (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

[TABLE]

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 Fig. 2 a). In Fig. 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 $\nu_{\mathrm{ss}}/R^{2}$, the analytical result underestimates the correct value. This is because the bond length in the BD simulation is obtained from Eq. (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 $\nu(t)$ approximates very well the value obtained from BD simulations for large enough $R$, as shown in Fig. 2 c). In Fig. 2 d), we depict the monotonous rapid divergence of the stationary variance (24) as the time delay $\tau$ 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 $\omega_{\tau}$. The best localized structure is obtained for the frequency fulfilling the formula $\cos(\omega_{\tau}\tau)=\omega_{\tau}\tau$, i.e. $\omega_{\tau}\approx 0.74/\tau$ and thus $2k\tau/\gamma\approx 0.74\exp(2D\tau/R^{2})$. For the corresponding value 0.74 of the reduced delay $\omega_{\tau}\tau$ the system is in the dynamical regime (ii) with exponentially decaying oscillations.

2.3 Trimer

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

[TABLE]

For the relative coordinates $\mathbf{r}_{12}(t)=\mathbf{r}_{1}(t)-\mathbf{r}_{2}(t)$ we obtain

[TABLE]

where $\omega=2k/\gamma$ and similarly for $\dot{\mathbf{r}}_{13}(t)$ and $\dot{\mathbf{r}}_{32}(t)$. To get analytical results for bond lengths $r_{ij}(t)=|\mathbf{r}_{ij}(t)|$, we multiply the formulas for $\dot{\mathbf{r}}_{ij}(t)$ by the corresponding unit vectors $\mathbf{e}_{ij}(t)=\mathbf{r}_{ij}(t)/|\mathbf{r}_{ij}(t)|$ and linearize the resulting equations. To this end, we need to deal with the expressions $\mathbf{e}_{\alpha}(t-\tau)\cdot\mathbf{e}_{\beta}(t)$, $\alpha,\beta=\mathrm{I},\dots,\mathrm{III}$, where we introduced roman numbers as a shorthand indexing $\mathrm{I}\equiv 12,$ $\mathrm{II}\equiv 13,$ $\mathrm{III}\equiv 32$. For a vanishing delay $\tau=0$, $\mathbf{e}_{\alpha}\cdot\mathbf{e}_{\alpha}=1$ and the scalar products $\mathbf{e}_{\alpha}\cdot\mathbf{e}_{\beta}$ describe the angles of the triangle formed by the three particles (see Fig. 1 in Sec. 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 $\pi/3$, leading to the relations $\mathbf{e}_{\mathrm{I}}\cdot\mathbf{e}_{\mathrm{II}}=\mathbf{e}_{\mathrm{I}}\cdot\mathbf{e}_{\mathrm{III}}=-\mathbf{e}_{\mathrm{II}}\cdot\mathbf{e}_{\mathrm{III}}\approx 1/2+O[(r_{\alpha}-r_{\beta})/R]$ with a correction that is on the order of $(r_{\mathrm{I}}-r_{\mathrm{II}})/R$ for $\mathbf{e}_{\mathrm{I}}\cdot\mathbf{e}_{\mathrm{II}}$ and similarly for the other scalar products. The linearized equation for the relative coordinate $x_{\alpha}\equiv r_{\alpha}-R$ thus reads

[TABLE]

where the lower index $\alpha\equiv\alpha\mod\mathrm{III}$ is considered as periodic with the period 3, i.e., $x_{\mathrm{IV}}\equiv x_{\mathrm{I}}$ and $x_{\mathrm{V}}\equiv x_{\mathrm{II}}$. Similarly as in the case of the dimer, Eq. (27) describes the dynamics of $x_{\alpha}(t)$ well for large equilibrium bond lengths $R$ and for time delays small compared to reorientation times of the unit vectors $\mathbf{e}_{\alpha}$.

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

[TABLE]

for the three-dimensional column vector $\mathbf{x}(t)=\left[x_{\mathrm{I}}(t),x_{\mathrm{II}}(t),x_{\mathrm{III}}(t)\right]^{\intercal}$. In Eq. (28), the noise vector $\bm{\xi}(t)$ is given by $\bm{\xi}(t)\equiv A^{1}\bm{\eta}^{1}(t)+A^{2}\bm{\eta}^{2}(t)$ with the auxiliary noise vectors $\bm{\eta}^{j}(t)\equiv[\eta_{1}^{j},\eta_{2}^{j},\eta_{3}^{j}]^{\intercal}$, $j=1,2$, containing the $j$th components of the original noises $\bm{\eta}_{i}(t)$, $i=1,\dots,3$. From the system (27) follows that the matrices $M$, $A^{1}$ and $A^{2}$ read

[TABLE]

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

[TABLE]

Due to the linearity of Eq. (28) and Gaussianity of the noise, the Green’s function for the one-time PDF $P_{1}(\mathbf{r},t|\mathbf{r}_{0},0)$ is Gaussian [37], and determined by the mean value $\bm{\mu}(t)=\left<\mathbf{x}(t)\right>$ and the covariance matrix $\mathcal{K}(t)=\left<\mathbf{x}(t)\mathbf{x}^{\intercal}(t)\right>-\bm{\mu}(t)\left[\bm{\mu}(t)\right]^{\intercal}$. We review in detail the derivation of these functions in A and B.2.

For the initial condition $\mathbf{x}(t)=0$, $t<0$ and $\mathbf{x}(0)=\mathbf{x}_{0}$ we get

[TABLE]

where $\lambda(t)$ denotes the Green’s function for Eq. (28) given by

[TABLE]

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

[TABLE]

where $\mathcal{I}$ denotes the identity matrix. This formula follows from the results of C after substituting the matrices $\omega$ and $\sigma\sigma^{\intercal}$ from the formula (107) in the appendix by $\omega M$ and $4DM$. In the appendix, we also derive an analytical expression for the stationary time correlation matrix $C(t)=\lim_{s\to\infty}\left<\mathbf{x}(s)\mathbf{x}^{\intercal}(s+t)\right>$. The regime of stability $3\omega\tau/2<\pi/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 Fig. 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 $\omega$ and delay $\tau$ 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 $\tau$ and exhibit a minimum as functions of $\omega$, 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.

3 Structure Formation

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 $\omega\tau>\pi/2$ and the trimer for $\omega\tau>\pi/3$. Let us now discuss how general the presented findings are.

The stationary average bond length $\mu$ can be determined by minimizing the potential energy $V=\frac{1}{2}\sum_{(i,j)}V_{2}(r_{ij})$ with the two-particle potential $V_{2}$ given by Eq. (1). Minimizing the potential in our two-dimensional geometry yields the highly symmetric molecule-like structures shown in Fig. 4 a). Due to the confinement to 2d, the global minimum of the potential corresponding to $\mu=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 $\omega\tau$ 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 $\omega\tau$ causes exponentially diverging oscillations.

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 Fig. 4 b) and c), respectively. In the regime (i), we observe that larger systems relax faster than those with smaller $N$. Furthermore, in Fig. 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 $\omega\tau$. More precisely, all the curves in Fig. 4 c) are plotted using the same parameters except for $N$ and, while the curves for $N\leq 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 $(\omega\tau)^{\star}$ determining the threshold between the regimes (ii) and (iii) for $N=2,...,10$. In Fig. 5, we show its rescaled value

[TABLE]

where the coefficient $2/\pi$ 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 [$(\omega\tau)^{\star}=\pi/2$] and trimer [$(\omega\tau)^{\star}=\pi/3$]. However, the analytical results would imply that the constant $C$ equals to 2, which is smaller than the value $C\approx 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 extensive scaling $c_{\mathrm{crit}}\propto 1/N$ of the stability of the system could (almost) completely be removed by rescaling the potential stiffness according to $k\rightarrow k/N$ as it is often done in thermodynamical system with long-range interactions [64].

To better understand this behavior, let us first consider the approximate analytical model described in Sec. 2.2. Imagine a particle in a harmonic potential centered around $x=0$, which is initially located at $x(0)=x_{0}>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\in[0,\tau]$ such that $x(t)=x_{0}$ for all $t\in[0,\tau]$. In the subsequent time interval $t\in[\tau,2\tau]$, the particle experiences the force $F=-\omega_{\tau}x_{0}$ pushing it towards the opposite wall of the trap. For times $t>2\tau$ the force starts changing dynamically according to the earlier position at time $t-\tau$. The particle keeps its direction of motion until it reaches the position $x_{1}$ where the force changes its sign. For large delays, the particle may stop significantly later than crossing the minimum, so that $x_{1}<0$ and $|x_{1}|>|x_{0}|$. A similar process then repeats when the particle returns back, with the difference that now it reaches a maximum position $x_{2}>0$, $|x_{2}|>|x_{1}|>|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 Fig. 7 in Sec. 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 $\omega\tau$ 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 Ref. [11], our analysis thus reveals two qualitative differences. First, the structures formed in the quasi-constant velocity model oscillate for arbitrary nonzero delay $\tau$, while in the harmonic model these oscillations appear only if $\omega\tau$ 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.

4 Entropy Fluxes

The investigated model, much as the model discussed in Ref. [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 Ref. [11], the analysis of entropy flows in the supplementary material therein was performed for vanishing delay only. Using the approximate Gaussian PDFs found in Secs. 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({\bf x})$ and thus with a time constant configurational entropy $\mathcal{\mathcal{S}}=-k_{\rm B}\int{\rm d}\mathbf{x}P(\mathbf{x})\log P(\mathbf{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 Eqs. (3) and (4) for $\tau>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 Ref. [56]. The entropy production $\dot{\mathcal{S}}_{\rm F}=q_{\rm H}/T$ is the housekeeping heat flux $q_{\rm 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 $\eta_{\rm 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 Refs. [11, 65]. The center of mass coordinate of the system is not affected by the feedback and diffuses freely (see Sec. 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(\mathbf{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}\mathbf{x}P(\mathbf{x},t)\log P(\mathbf{x},t)$ can formally be written as

[TABLE]

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 $\mathbf{x}(t)$ describing the bonds obeys the generalized Langevin equation $\dot{\mathbf{x}}(t)=\mathbf{F}[\mathbf{x}(t),\mathbf{x}(t-\tau)]+\sigma\eta(t)$, where $\eta$ denotes a zero mean Gaussian white noise with the covariance matrix $\left<\eta_{i}(t)\eta_{j}(t^{\prime})\right>=\delta_{ij}\delta(t-t^{\prime})$, the dynamical equation for $P(\mathbf{x},t)$ can be written in the form [30, 66]

[TABLE]

In this equation, the fist term on the right stands for the diffusive spreading of the PDF. The term $\mathcal{L}[\mathbf{x},t]$ corresponds to the time-delayed force $\mathbf{F}[\mathbf{x}(t-\tau)]$ 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, 66] for more details about its structure.

Inserting Eq. (37) into the formal time derivative $\dot{\mathcal{S}}(t)=-k_{\rm B}\int{\rm d}\mathbf{x}[\partial_{t}P(\mathbf{x},t)]\log P(\mathbf{x},t)$ of the system entropy ${\mathcal{S}}(t)$, we find that

[TABLE]

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({\bf 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, 67] and use Eq. (39) to define the stochastic (position dependent) entropy flux

[TABLE]

The average flux (39) then follows as the average $\dot{\mathcal{S}}_{-}(t)=\left<\dot{s}_{-}(\mathbf{x},t)\right>$ either over the PDF $P(\mathbf{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

[TABLE]

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

[TABLE]

leading to the rate of change

[TABLE]

From Eq. (38) we then find that

[TABLE]

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

Let us now focus on the specific setups considered in Secs. 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^{\intercal}=4D$ and ${\mathcal{K}}(t)=\nu(t)$, we get

[TABLE]

and

[TABLE]

where $\lambda(t)$ is given by Eq. (18) and the variance reads $\nu(t)=4D\int_{0}^{t}\mathrm{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^{\intercal}=4DM$ and Jacobi’s formula for the derivative of determinants, we obtain the expressions

[TABLE]

and

[TABLE]

where $\lambda(t)$ is given by Eq. (33) and the covariance matrix ${\mathcal{K}}(t)$ reads $4DM\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 $\nu(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 (85) for $\lambda(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. $\lim_{t\to\infty}{\mathcal{S}}(t)={\mathcal{S}}$, $\lim_{t\to\infty}\dot{\mathcal{S}}(t)=0$ and $\dot{\mathcal{S}}_{-}\equiv\lim_{t\to\infty}\dot{\mathcal{S}}_{-}(t)=-\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 $\nu_{\rm ss}$ and trimer covariance ${\mathcal{K}}_{\rm ss}$, we find from Eqs. (45) – (48)

[TABLE]

for the dimer and

[TABLE]

for the trimer. In the last two formulas, $a={\mathcal{K}}_{\rm ss}(1,1)$ denotes diagonal and $b={\mathcal{K}}_{\rm 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=\gamma D/k_{\rm B}$ and with the delay time $\tau$, and exhibit a minimum as functions of the frequencies $\omega_{\tau}$ (dimer) and $\omega$ (trimer), similarly as the variance $\nu_{\rm ss}$ and the diagonal matrix elements of the covariance matrix ${\mathcal{K}}_{\rm 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 Ref. [11]. This means that the fluxes are discontinuous in the formal limit $T\to 0$ because they must inevitably vanish for zero noise, where the PDF for the system evolution is a $\delta$-function for all times. We plot the stationary entropy fluxes (50) and (52) as functions of the delay $\tau$ in Fig. 6 a) and frequency $\omega$ in Fig. 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 Fig. 6 b), arises as a result of a competition between stronger confinement, corresponding to larger frequencies $\omega$, and gradual destabilization with increasing $\omega\tau$, when the system enters the unstable regime (iii). For the dimer, it is located at the maximum of the variance (24) defined by the equation $\cos(\omega_{\tau}\tau)=\omega_{\tau}\tau$ since $\dot{\mathcal{S}}_{-}=-\nu_{\rm SS}/2D$. This result is intuitive since more localized distributions spread due to the diffusion faster than broader ones. For the figure, we used for simplicity the approximation $\omega_{\tau}\approx\omega$ for the dimer.

5 Transition Rates for Isomer Transformations

The particles within the molecular structures depicted in Fig. 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 $\tau$ 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\to B}(t)$ from $A$ to $B$ during an infinitesimal time interval $\Delta t$ as $\kappa_{A\to B}(t)=N_{A\to B}(t)/\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\rightarrow 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\rightarrow B}(t)=\dot{S}_{A\rightarrow B}(t)/S_{A\rightarrow 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 $\lim_{t\to\infty}\kappa_{A\to B}(t)$ of the transition rate can (approximately) be calculated using Kramers’ rate theory [23, 22] 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 generalized Langevin equations (GLE’s) describing non-Markovian systems. A crucial contribution in this direction came from Grote and Hynes [68] who derived a dynamical correction to Kramers’ result. While their analysis was based on a parabolic barrier, Pollak [69] 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 Ref. [70]. Based on their small-delay approximation, Guillouzic et al. [71] 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. [72] 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=|\mathbf{r}_{1}-\mathbf{r}_{2}|>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(|\tilde{r}|-R)^{2}/2$, depicted in Fig. 7, with the diffusion coefficient $2D$.

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 (91), which reads

[TABLE]

This equation describes diffusion in the harmonic potential

[TABLE]

with the time-dependent stiffness $\gamma\omega_{\tau}(t)$ [given by (89)] and the time-dependent diffusion coefficient $2D_{\tau}(t)$ [given by (93)].

The validity of Eq. (53) for $P_{1}$ 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 $\kappa(t)$ for switching between the isomers. In the Markovian case, the transition rate for surpassing the (effective) energy barrier $\gamma\omega_{\tau}(t)R^{2}/2$ at $r=0$ to the other isomer can be calculated from Eq. (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 $\gamma\omega R^{2}/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>0$.

Numerical Method

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

[TABLE]

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  [73]. Here, $P_{\rm a}(x,t)$ is the solution to the FPE (53) with absorbing boundary at $x=-R$ and $S_{\rm a}(t)=\int_{-R}^{\infty}\mathrm{d}xP_{\rm a}(x,t)$ is the probability that the particle survives in the right well until time $t$. Equation (55) follows from Eq. (53) by inserting the definitions of the PDF $\tilde{P}_{\rm a}$ and of the transition rate

[TABLE]

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

Bullerjahn’s Method

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

[TABLE]

across the absorbing boundary by 333Rescaling the diffusion coefficient by factor 2 corrects for the part of the diffusive flux that can not return to the system due to the absorbing boundary, see Ref. [24] for details.

[TABLE]

and the survival probability $S_{\rm a}(t)$ by

[TABLE]

In the last expression, the symbol $\mathrm{Erf}$ denotes the error function and the variance $\nu(t)$ is given by Eq. (22). The approximate Markovian transition rate is then given by

[TABLE]

In Fig. 8, we show the frequency $\omega_{\tau}(t)$, the (effective) diffusion coefficient $2D_{\tau}(t)$, survival probabilities $S_{\rm a}(t)$ and $S_{1}(t)$ and the transition rates $\kappa_{\rm M}(t)$ and $\kappa_{\rm 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(t)^{2}+\omega_{\tau}(t)\nu(t)$ and the rates saturate with time. The relaxation time of $\omega_{\tau}(t)$ is determined by the time in which the Green’s function $\lambda(t)$ approaches the long-time exponential representation (85), and the corresponding stationary value, $\omega_{\tau}(\infty)=1/t_{\rm R}$, is controlled by the relaxation time $t_{\rm R}$ for decay of $\lambda(t)$ to [math], see also Fig. 12 in A. The effective diffusion coefficient converges to the value $D_{\tau}(\infty)=2D[1+\sin(\omega_{\tau}\tau)]/\cos(\omega_{\tau}\tau)\geq 2D$, determined by the stationary variance $\nu(\infty)=\nu_{\mathrm{ss}}$, see Eq. (24). The transition rates relax with the relaxation time $t_{\rm R}$, similarly as the corresponding PDFs $P_{\rm a}$ and $P_{1}$. The analytical expressions for $S_{1}(t)$ and $\kappa_{\rm B}(t)$ approximate the numerical results for $S_{\rm a}(t)$ and $\kappa_{\rm M}(t)$ best for short times $t\ll t_{\rm R}$, where the PDFs $P_{\rm a}$ and $P_{1}$ 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 $\kappa_{B}$ overestimates the corresponding exact rate $\kappa_{M}$, see also Figs. 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 $\kappa_{B}<\kappa_{M}$, see Fig. 9 c).

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

[TABLE]

Since the denominator in the expression for the rate $\kappa_{\rm M}(t)$ is smaller than that for $\kappa_{\rm B}(t)$, we conclude that, for low energy barriers, the inequality $\kappa_{\rm B}(t)\leq\kappa_{\rm M}(t)$ holds. On the other hand, for very high barriers, the PDFs $\tilde{P}_{\rm a}$ and $P_{1}$ 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)\approx 1$, and $j^{*}\approx j-2D\partial_{x}P_{1}|_{x=-R}<0$ leads to the inequality $\kappa_{\rm B}(t)\geq\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+2D_{\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 $\partial_{t}\equiv\partial/\partial t$ and $\partial_{x}\equiv\partial/\partial x$. Then the transition rates $\kappa_{\rm M}(\infty)=-j_{\rm Da}$ and $\kappa_{\rm B}(\infty)=-j_{\rm D1}$ are determined solely by the diffusive fluxes $j_{\rm Da}\equiv 2D_{\tau}(\infty)\partial_{x}\tilde{P}_{\rm a}(x,\infty)|_{x=-R}$ and $j_{\rm D1}\equiv 2D_{\tau}(\infty)\partial_{x}\tilde{P}_{\rm 1}(x,\infty)|_{x=-R}$. The smaller the frequency $\omega$, the wider the PDFs $\tilde{P}_{1}$ and $\tilde{P}_{\rm a}$. For small $\omega$, 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}_{\rm 1}$, leading to $|j_{\rm D1}|<|j_{\rm Da}|$ and $\kappa_{\rm B}(\infty)\leq\kappa_{\rm M}(\infty)$, in accord with the argument put forward in the previous paragraph. With increasing $\omega$, 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 [73, 75], 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 $\omega$, the tail of $\tilde{P}_{\rm a}$, with small derivative (small $j_{\rm 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)\geq\kappa_{\rm M}(\infty)$. Finally, for very stiff traps ($\omega\to\infty$), both $j_{\rm Da}$ and $j_{\rm 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 $\lim_{t\to\infty}\omega_{\tau}(t)$ and $\lim_{t\to\infty}D_{\tau}(t)$ exist. Also in regime (ii), the PDF $P_{1}$ assumes, after long times, the time-independent stationary form $P_{1}=\exp\left(-x^{2}/2\nu_{\rm ss}\right)/\sqrt{2\pi\nu_{\rm ss}}$ with the variance $\nu_{\rm ss}$ given by Eq. (24). This suggests that, also in this regime, the transition rate should saturate at long times. However, both the frequency $\omega_{\tau}(t)$ and the diffusion coefficient $2D_{\tau}(t)$ actually exhibit diverging oscillations caused by the oscillations in the Green’s function $\lambda(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.

Long-Time Behavior and Kramers’ Method

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

[TABLE]

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

[TABLE]

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 [23, 22] for the transition rate for leaving one of the wells of a cusped potential with barrier height $E_{\rm b}$,

[TABLE]

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 $\lim_{t\to\infty}\lambda(t)=0$, valid in the dynamical regimes (i) and (ii).

Renormalized Transition Rates

Interestingly, the term $\omega_{\tau}(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 $\kappa_{\mathrm{B}}$ and $\kappa_{\mathrm{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 Eq. (65) cannot be correct, due to the diverging oscillations in the time-dependent frequency $\omega_{\tau}(t)$. Nevertheless, the exponential term seems to be reasonable, and thus it is tempting to use in the prefactor of the transition rates simply $\omega$, instead of the problematic $\omega_{\tau}(\infty)$. This substitution gives the correct pre-exponential factor of the rate for vanishing delay $\tau=0$, where $\omega_{\tau}(\infty)=\omega$ and $2D_{\tau}(\infty)=2D$. We denote the rates with the renormalized prefactor as

[TABLE]

with $\mathrm{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 Eq. (63) is correct for Markov dynamics ($\tau=0$), it can not be precisely valid for the time-delayed dynamics ($\tau>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_{\rm b}=\gamma\omega R^{2}/2$, but the barrier with energy $\gamma\omega[r(t-\tau)]^{2}/2$.

In Figs. 9 a) – c), we compare the various analytical predictions $\kappa_{\mathrm{x}}$ and $\tilde{\kappa}_{\mathrm{x}}$, $x=$ M, B or K, for the long-time transition rate (lines) with the asymptotic rate $\kappa$ (symbols) calculated from BD simulations of Eq. (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 $\kappa$ can be predicted reasonably well only for values of $\tau$ in the dynamical regime (i) ($\omega\tau<1/{\rm e}\approx 0.37$) and in the first part of the dynamical regime (ii) ($\omega\tau<\pi/2\approx 1.57$). The rate $\kappa$ is best approximated by the expression $\tilde{\kappa}_{\mathrm{M}}$ obtained numerically from Eq. (63) with $\omega$ substituted for $\omega_{\tau}(\infty)$ in the operator $\mathcal{L}$. From the analytical expressions the re-scaled Bullerjahn expression $\tilde{\kappa}_{\mathrm{B}}$ [see Eqs. (64) and (66)] works best. However, in the figure we have used parameters leading to a high barrier $E_{\rm b}$, and thus Kramer’s and Bullerjahn’s predictions, $\kappa_{\mathrm{K}}$ and $\kappa_{\mathrm{B}}$, almost coincide. As a consequence, the line for $\tilde{\kappa}_{\mathrm{K}}$ in the figure overlaps with $\tilde{\kappa}_{\mathrm{B}}$, similarly as the line $\kappa_{\mathrm{K}}$ (suppressed in the figure for better readability) overlaps with $\kappa_{\mathrm{B}}$.

Delay-Dependence of $\kappa$

In Fig. 9 d), we also show the behavior of the rate $\kappa$ in the parameter regime (iii) inaccessible to the analytical and numerical formulas due to the diverging oscillation. The simulated transition rate in Figs.  9 a) – d) is first approximately exponential and thus convex in $\tau$, then its curvature changes to concave and it runs through a maximum and, finally, the rate starts to decrease. The value of $\tau$ where the curvature changes sign coincides with the boundary $\pi/2\omega\tau\approx 1.57$ between the dynamical regimes (ii) and (iii). Interestingly, no qualitative change of $\kappa(\tau)$ is observed at the boundary $\omega\tau=1/{\rm e}\approx 0.37$ between the regimes (i) and (ii). It is tempting to attribute, the (approximate) exponential increase of the rate with $\tau$ in regimes (i) and (ii) to the increase of the steady-state variance $\nu_{\rm ss}$, which is given as the ratio of the effective energy barrier $V(-R,\infty)$ and the effective thermal energy $\gamma/2D_{\tau}(\infty)$. Although this explanation may work well for small delays, it breaks down for values of $\tau$ in the second half of the regime (ii), where the actual rate $\kappa$ is no longer well approximated by our analytical and numerical predictions. This means that the identification of the parameters $\gamma\omega_{\tau}(t)$ and $2D_{\tau}(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 $\tau$ only. In the dynamical regime (iii), the particle undergoes oscillations with an amplitude that increases both with $\tau$ and $t$. The corresponding transition rate, obtained from the BD simulations, thus decreases with the delay $\tau$ 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 Fig. 10 with the three distinguishable particles labeled by the numbers 1 to 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.

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 Fig. 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 $\kappa$ using the angle $\varphi$ between the abscissas $|12|$ and $|13|$ and a neutral region of width $\Delta x=\sqrt{3/16}R$ as exemplarily shown in Fig. 10. A neutral state $\left|0\right>$ is introduced to avoid over counting due to fluctuations of $\varphi$ around [math] and $\pi$ . It is occupied if the smallest height of the triangle formed by the three particles, is smaller $\Delta x/2$. i.e. either if $|\varphi|<\phi_{>}\equiv\max\left[\arcsin(\Delta x/2r_{12}),\arcsin(\Delta x/2r_{13})\right]$ or $|\varphi|>\phi_{<}\equiv\arccos(\Delta x/2r_{12})+\arccos(\Delta x/2r_{13})$. For $\varphi\in[\phi_{>},\phi_{<}]$ the system is said to be in the clockwise state $\left|1\right>$, while for $\varphi\in[-\phi_{>},-\phi_{<}]$ it is in the counter-clockwise state $\left|-1\right>$. To calculate the transition rate, we have counted the number of transitions between the states $\left|\pm 1\right>$ during a specific simulation time window, where the transition occurred if the system underwent the sequence of states $\left|1\right>\to\left|0\right>\to\left|-1\right>$ or $\left|-1\right>\to\left|0\right>\to\left|1\right>$.

The resulting transition rate $\kappa$ as a function of the time delay $\tau$ is depicted in Fig. 11 a). Therein, one can observe the four dynamical regimes described in Sec. 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 Fig. 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.

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_{\mathrm{b}}$ of the particle 2 from the $x$-axis and construct an appropriate effective energy barrier $E_{\mathrm{b,eff}}$ and diffusion coefficient $D_{\mathrm{eff}}$. In the steady state, the three particles are most likely found at the vertices of an equilateral triangle. Using the notation of Sec. 2.3, we express the vector $\mathbf{r_{b}}$ from the particle 2 to the center of the abscissa $|13|$ as $\mathbf{r}_{b}=\mathbf{r}_{12}+\mathbf{r}_{31}/2$. Hence the Gaussian white noise $\bm{\eta}_{\mathrm{b}}(t)$ corresponding to the coordinate $\mathbf{r}_{b}$ is given by

[TABLE]

where we have used the relations (4). Based on these considerations, we approximate the effective diffusion coefficient by $D_{\mathrm{eff}}=3D$.

For the corresponding effective energy barrier, we make the ansatz $E_{\mathrm{b,eff}}=C(\tau)\gamma\omega R^{2}/12$. Namely, we first consider the minimum value $E_{\mathrm{b}}=kR^{2}/6=\gamma\omega R^{2}/12$ of the energy difference between a configuration with all the three particles aligned ($V=kR^{2}/6$) and the configuration where the particles form an equilateral triangle with side length $R$ ($V=0$), where $V$ is given by Eq. (2). The (possibly delay-dependent) unknown dimensionless prefactor $C(\tau)$ 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>$ to $\left|-1\right>$. The resulting Kramers’ rate for the transition from $\left|1\right>$ to $\left|-1\right>$ thus reads

[TABLE]

where $a(\tau)$ is a further unknown prefactor that may depend on all parameters of the model [see, for example, Eq. (65) for Kramers’ rate in a cusped potential].

In order to test the formula (68), we have fitted it to the transition rate $\kappa$ obtained from BD simulations as a function of the diffusion coefficient $D$. The results for different values of the time delay are shown in Fig. 11 b) and the corresponding values of the coefficients $C(\tau)$ and $a(\tau)$, obtained from the fits, are given in Tab. 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 $\kappa$ can be relatively well described by the Kramers-type ansatz (65).

Our attempts to include the effects of the delay in the effective diffusion coefficient and energy barrier, analogously to the approach described in Sec. 5.1, where we made the replacements $\omega\rightarrow\omega_{\tau,\mathrm{ss}}$ and $D\rightarrow D_{\tau,\mathrm{ss}}$, did not lead to a significant change of the value for $C(\tau)$. We thus conclude that the deviations of $C(\tau)$ 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.

6 Extensions to other memory kernels

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 (86) is equivalent to the generalized Langevin equation (GLE)

[TABLE]

with positive frequency $\omega>0$ and the memory kernel given by $\phi(t)=\delta(t-\tau)$, $\tau>0$. Deriving a time-local Langevin equation from the GLE (69) with arbitrary $\phi(t)$ along the lines of B.1, we obtain Eq. (88) with $\lambda(t)$ being the Green’s function for Eq. (69), i.e. solving (69) with the initial condition $\lambda(0)=1$ and $\lambda(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 $\phi(t)$ after substituting the Green’s function (83) for $\phi(t)=\delta(t-\tau)$ by the Green’s function corresponding to the chosen $\phi(t)$. In the rest of this section, we review some paradigmatic examples of memory kernels $\phi(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 $\phi(t)=\delta(t-\tau)$ are systems with multiple different time delays with the memory kernel

[TABLE]

Properties of these systems are studied in Ref. [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 $\phi(t)=\delta(t+\tau)$, $\tau>0$ and it can be solved using the strategy described in A. The resulting Green’s function

[TABLE]

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

[TABLE]

where $s$ is the absolute value of the smallest root of the equation $\rho=-\omega e^{\tau\rho}$.

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

[TABLE]

One of the reasonable strategies for predicting $x_{\mathrm{pre}}$ is to use the linear extrapolation

[TABLE]

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

[TABLE]

which is solved by $x(t)=x_{0}\lambda(t)$ with the exponentially decaying Green’s function

[TABLE]

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 $\tau<0$ in Eq. (73) ], the presented first order approximation predicts dynamics with rescaled frequency $\omega_{\rm r}$ increasing with increasing time delay (for $\omega\tau<-1$). One might conclude that such dynamics would lead to a smaller variance $\nu_{ss}$ for non-zero delays than for a vanishing delay, which contradicts our exact result (24). However, performing the first-order expansion in the equation $\dot{x}(t)=-kx(t-\tau)+\sigma\eta(t)$ including the noise term, one finds that both the rescaled frequency and diffusion coefficients are smaller than in the original equation and the resulting variance equals the first-order expansion of the exact result given in Eq. (24). For details and further discussion of Taylor expansions in SDDE, see Ref. [30].

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

[TABLE]

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

[TABLE]

where $\Omega=\sqrt{\omega^{2}-b^{2}/4}$, which is reminiscent of the Green’s function (85) for the system with conventional time delay [$\phi(t)=\delta(t-\tau)$]. The difference is that, while the Green’s function (78) always decays exponentially with time, the Green’s function (85) 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 Ref. [78].

7 Conclusion and outlook

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 Fig. 1 in Sec. 2. The system is described by the set (3) of $2N$ non-linear delayed Langevin equations and hence its dynamics is non-Markovian. At long times, the particles form highly symmetric dynamical molecular-like structures, depicted in Fig. 4 a) in Sec. 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 Brownian dynamics 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 Sec. 2).

Our analytical results for the dimer and trimer imply that these structures are stable only for small enough values of the product $k\tau$, 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\geq 2$ by BD simulations, described in Sec. 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\tau$ than predicted from the linearized equations, the critical value of the product $k\tau$ 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\to 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 Ref. [11], did not lead to unstable behavior.

Interpreting the inter-particle interactions as an action of a feedback control mechanism, we have, in Sec. 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 Fokker-Planck equation (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 $\tau$. 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 Eq. (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 Sec. 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. $\phi(t)=\delta(t-\tau)$. It is enough to substitute the Green’s function $\lambda(t)$ (83) 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.

Acknowledgments

We acknowledge funding by Deutsche Forschungsgemeinschaft (DFG) via SPP 1726/1. 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.

Appendix A Solution of the Noiseless Problem

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

[TABLE]

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

[TABLE]

with $\mathbf{x}_{0}\equiv\mathbf{x}(0)$. The solution of this equation for the Laplace transformed variable reads

[TABLE]

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 $e^{-s\tau}/s^{k+1}$ is given by $(t-\tau)^{k}\theta(t-\tau)/k!$ [79] and thus the formula (81) can be inverted. Finally, the solution of Eq. (79) is given by

[TABLE]

where

[TABLE]

is a matrix-valued function which solves Eq. (79) with the initial condition $\mathbf{x}(t)=0$ for all $t<0$ and $\mathbf{x}(0)=\mathcal{I}$. In the present paper, we always assume that the system is initialized at time $t$ at position $\mathbf{x}_{0}$ witha special history, namely $\mathbf{x}(t)=0$ for all $t<0$, allowing us to simplify Eq. (82) to $\mathbf{x}(t)=\lambda(t)\mathbf{x}_{0}$.

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

[TABLE]

where $W$ denotes the matrix valued Lambert $W$ function [81]. The Lambert $W$ function is a multivalued complex function. The long-time behavior of solutions to Eq. (79) 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 $\alpha$. The corresponding values of the Lambert $W$ function strongly depend on the reduced delay $\tau\omega$.

For example, in one dimension, where $\omega$ 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 $\tau\omega$ leading to three different dynamical regimes of solutions

[TABLE]

$1/t_{\rm R}=\Re(\alpha)$, $\Omega=\left|\Im(\alpha)\right|$, to Eq. (79). The boundaries between these regimes can be determined analytically [25, 27, 30]: (i) $0<\tau\omega\leq 1/e$, where $\alpha$ is real and positive and $x(t)$ decays exponentially to [math] with lifetime $t_{\rm R}$; (ii) $1/e<\tau\omega<\pi/2$, where $\alpha$ is complex with a positive real part, producing exponentially damped oscillations of $x(t)$ with frequency $\Omega$ and lifetime $t_{\rm R}$; and (iii) $\pi/2<\tau\omega$, where $\alpha$ is complex with a negative real part, corresponding to exponentially diverging oscillations of $x(t)$ with frequency $\Omega$. For $\tau\omega=\pi/2$, $1/t_{\rm R}=0$ and $\lambda(t)$ oscillates with frequency $\Omega$ without any decay.

In the panel a) of Fig. 12, we show that for long times the Green’s function $\lambda(t)$ for Eq. (79) is well approximated by the exponential solution (85). The above described dynamical regimes are reflected in the behavior of the decay rate $1/t_{\rm R}$ and the frequency $\Omega$ 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 C, which is also well described by the formula (85).

Appendix B Solution with Noise

B.1 One Dimension

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

[TABLE]

where $\eta(t)$ is a white noise fulfilling $\left<\eta(t)\right>=0$ and $\left<\eta(t)\eta(t^{\prime})\right>=\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 Fokker-Planck equation for the PDFs for $x(t)$. In order to do that, we write the formal solution of Eq. (86) as

[TABLE]

where $\lambda(t)$ is given by Eq. (83) and $\psi(t)=-\omega\int_{-\tau}^{0}\mathrm{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 (87) can already be used for the calculation of the time-correlation function for $x(t)$. Here, we differentiate the solution (87) which yields the time-local Langevin equation

[TABLE]

Due to the time-nonlocal nature of Eq. (86), the potential in the time-local equation (88) possess the time-dependent stiffness

[TABLE]

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, Eq. (88) includes the Gaussian colored noise $\xi(t)=\lambda(t)\frac{\mathrm{d}}{\mathrm{d}t}\int_{0}^{t}\mathrm{d}s\frac{\lambda(t-s)}{\lambda(t)}\eta(s)$ which satisfies $\langle\xi(t)\rangle=0$ and

[TABLE]

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_{0}$ at time $t_{0}$ 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 (88) 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_{0}$ at time [math], are found to be

[TABLE]

Similarly as the trap stiffness $\omega_{\tau}(t)$, also the effective diffusion coefficient, corresponding to the time-local description, is time dependent if $\tau>0$. It reads

[TABLE]

with the variance $\nu(t)=\sigma^{2}\int_{0}^{t}\mathrm{d}s\lambda(s)^{2}$ and $c(t,t^{\prime})\equiv\sigma^{2}\lambda(t)\frac{\mathrm{d}}{\mathrm{d}t}\int_{0}^{t^{\prime}}\mathrm{d}s\lambda(t-s)\lambda(t^{\prime}-s)/\lambda(t)$.

Because of the oscillatory nature of $\lambda(t)$ in the dynamical regimes (ii) and (iii), the coefficients $\omega_{\tau}$, $b$, $c$ and $D_{\tau}$ in the FPEs (91) and (92) 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 Eqs. (19) and (20), are always reasonable [35].

B.2 Higher Dimensions

Let us now consider the problem

[TABLE]

with general matrices $\omega$ and $\sigma$ and the vector $\bm{\eta}(t)$ of white noises fulfilling $\left<\bm{\eta}(t)\right>=0$ and $\left<\eta_{i}(t)\eta_{j}(t^{\prime})\right>=\delta_{ij}\delta(t-t^{\prime})$. Since this system of Langevin equations is linear, the one- and two-time probability distributions $P_{1}(\mathbf{x},t|\mathbf{x}_{0},0)$ and $P_{2}(\mathbf{x},t;\mathbf{x}^{\prime},t^{\prime}|\mathbf{x}_{0},0)$ for $\mathbf{x}(t)$ defined in the preceding section must be Gaussian [63] as in the one-dimensional case and also the corresponding Fokker-Planck equations 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

[TABLE]

based solely on the formal solution of the Langevin system (94)

[TABLE]

with the initial condition $\mathbf{x}(t)=0$ for all $t<0$ and $\mathbf{x}(0)=\mathbf{x}_{0}$ and with the Green’s function $\lambda(t)$ of the Langevin system given by Eq. (83).

The mean value $\bm{\mu}(t)=\langle\mathbf{x}(t)\rangle$ and the elements $\mathcal{K}_{ij}(t)$ of the covariance matrix $\mathcal{K}(t)=\left<\mathbf{x}(t)\mathbf{x}^{\intercal}(t)\right>-\left<\mathbf{x}(t)\right>\left<\mathbf{x}^{\intercal}(t)\right>$ defining the PDF (95) can be obtained by inserting $\mathbf{x}_{i}(t)$ from Eq. (96) into the definitions and averaging over the noise $\eta(t)$. The results are

[TABLE]

and

[TABLE]

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

[TABLE]

The mean value $\bm{\mu}(t)$ is given by $\bm{\mu}(t)=\left<\mathbf{y}(t)\right>_{\mathbf{x}(t\leq 0)}$, and $\mathcal{K}(t)=\left<\mathbf{y}(t)\left[\mathbf{y}(t)\right]^{\intercal}\right>_{\mathbf{x}(t\leq 0)}+\int_{0}^{t}\mathrm{d}s\ \lambda(t-s)\sigma\sigma^{\intercal}\lambda^{\intercal}(t-s)$. The averages $\left<\bullet\right>_{\mathbf{x}(t\leq 0)}$ above are taken with respect to the PDF $P[\mathbf{x}(t),t\leq 0]$. The long-time behavior of the covariance matrix (98) is studied in the next section of this appendix.

Appendix C Time-Correlation Matrix and Stationary Covariance Matrix

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

[TABLE]

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 delay differential equation as the Green’s function $\lambda(t)$, i.e.

[TABLE]

The restriction $t>0$ for validity of this equation comes from discontinuity (and thus non-differentiability) of $\lambda(t)$ at $t=0$. The solution to Eq. (102) is given by Eq. (82) and hence the stationary space-time correlation matrix can be written as

[TABLE]

where $\mathcal{C}_{0}=\mathcal{C}(0)\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\in[-\tau,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\in[0,\tau]$ in one-dimension (see also Ref. [32]). They utilized the symmetry $\mathcal{C}(t)=\mathcal{C}(-t)$ following from Eq. (101) to rewrite the delay differential equation (102) as $\dot{\mathcal{C}}(t)=-\omega\mathcal{C}(\tau-t)$ for $t\in(0,\tau)$. Taking the derivative of this equation and using Eq. (102) yields the second order ordinary differential equation

[TABLE]

for the initial condition $\mathcal{C}(t)=\mathcal{C}(-t)$, $t\in[-\tau,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}=\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\in(0,\tau)$. Specifically, we show at the end of this appendix that $\mathcal{C}(\tau)=0.5\omega^{-1}\sigma\sigma^{\intercal}$ and $\dot{\mathcal{C}}_{0}=-0.5\sigma\sigma^{\intercal}$. Using these results, the final expression for the correlation matrix for $t\in[-\tau,\tau]$ reads

[TABLE]

with the initial value

[TABLE]

given by the stationary value of the covariance matrix. The whole time-dependence of $\mathcal{C}(t)$ for $t\geq 0$ is thus described by the formulas (103), (105), and (106). 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 $\omega$ and $\sigma$ stand for real numbers and, more generally, if the matrices $\omega^{-1}$ and $\sigma\sigma^{\intercal}$ commute, we can rewrite Eq. (106) as

[TABLE]

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

The expression for $C(\tau)$ can be obtained by multiplying Eq. (94) by $\mathbf{x}^{\intercal}(s)$, averaging the result over the noise, using the assumed stationarity of the process implying the formula $\dot{\mathcal{C}}(0)=\mathrm{d}\left[\lim_{s\to\infty}\left<\mathbf{x}(s)\mathbf{x}^{\intercal}(s)\right>\right]/\mathrm{d}s=0$, and applying the symmetry $\mathcal{C}(t)=\mathcal{C}(-t)$. The result is $\mathcal{C}(\tau)=\omega^{-1}\sigma\left<\bm{\eta}(s)\mathbf{x}^{\intercal}(s)\right>=0.5\omega^{-1}\sigma\sigma^{\intercal}$, where the last equality comes after inserting the formal solution (96) for $\mathbf{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 $\delta$-function only. The expression for $\dot{\mathcal{C}}_{0}$ then comes simply from Eq. (102), 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)$.