Brownian Molecules Formed by Delayed Harmonic Interactions

A time-delayed response of individual living organisms to information exchange within flocks or swarms leads to the formation of complex collective behaviors. A recent experimental setup by Khadka et al., Nat. Commun. 9, 3864 (2018), employing synthetic microswimmers, allows to realize 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 via linear stochastic delay-differential equations, and for $N>3$ by Brownian dynamics simulations. The particles form non-equilibrium molecule-like structures which become unstable with increasing number of particles, delay time, and interaction strength. We evaluate the entropy fluxes in the system and 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.


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 systems of interacting, active and energy consuming particles [4,5]. A special focus is devoted to emergent collective behavior [6]. Most of the corresponding 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 Janus particles [20,21] that are made from two half-spheres from different materials and propel in the direction of one of them if illuminated by a laser, or supplied by another energy source. In order to steer these particles, one thus 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, nevertheless, in Ref. [11], only 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 time-delayed harmonic interactions between the individual particles. Similarly to the case of Ref. [11], the considered two-dimensional N -particle system is described by a set of 2N coupled non-linear stochastic delay differential equations that induce non-Markovian dynamics. For small enough values of the delay, the particles form, after a transient period, highly symmetric non-equilibrium molecularlike structures that fluctuate due to thermal noise. The resulting structures strongly differ from the molecules created by constant velocities in Ref. [11] that oscillated due to the nonzero delay even for vanishing noise amplitude. Another difference between the two settings 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 allows us to linearize the underlying set of stochastic delay differential equations and, to a large extend, study the model 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 the mean values, covariance matrix, and time-correlation matrix. We verify the validity of these results by Brownian dynamics (BD) simulations of the full 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 value of the delay for which the system cannot form stable molecules decreases with increasing number of particles as 1/N . If we label the individual particles, we can distinguish between several isomers of the respective molecules according to different orderings of the particles. In the course of time, the noise induces jumps of a given molecule between the individual isomers. We use our analytical results for the dimer and for the trimer to evaluate the corresponding transition rates on the bases of Kramers' theory [22,23] and also of the more recent Bullerjahn's theory [24]. We compare the obtained results with the rates calculated using BD simulations and identify a useful formula for the transition rate that provides a good prediction for small and moderate values of the delay.

Stochastic delay differential equations
From a mathematical point of view, delay differential equations (DDE's) [25,26] in general generate rich dynamics [27]. The solutions may converge to fix points or limit cycles, behave chaotically, and they can even exhibit multistability [28]. For systems affected by noise, the DDEs are generalized to stochastic delay differential equations (SDDE) [29] inducing 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 with generally unknown closure of the hierarchy.
For non-linear systems there are three established approximate approaches how to tackle the infinite hierarchy: (i) The so called small delay approximation [30] is based on a Taylor expansion in the delay, making the problem time local. (ii) Assuming that the terms in the SDDE containing time delay are small and thus 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 Gaussian distributions. 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. Afterward, 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 rather complicated, there is a great interest in understanding DDE's and SDDE's, due to their broad range of applications. Prominent examples are population dynamics [41,42], where the delay results from maturation times, economics [43,44,45,46] according to the limited reaction times of the market participants, or engineering [47]. In biology, finite transition times can play a significant role, for example in physiological systems [48,49,50] and neural networks [51,52,53,28]. These effects are loosely related to the hereby concerned field of active matter as they may induce finite reaction times of the individual living agents. Recently, first efforts [54,55,56] 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.

Outline
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 performing an appropriate linearization of the underlying SDDEs, we study analytically both transient and stationary properties of the probability distributions for bond length in systems of two and three particles. The stability of larger systems with respect to delay is studied in Sec. 3. In Sec. 4, we evaluate the entropy out-flux from the system due to the feedback maintaining the non-equilibrium stationary structures. In Sec. 5 we utilize our analytical description of the dimer and trimer for analytical and numerical investigation of the transition rates for hopping between the individual isomers of the stationary structures, and we compare the results with transition rates obtained from BD simulation of the systems in question. In Sec. 6, we explain how the presented results can be extended for systems with arbitrary memory kernels. In particular, we focus on the cases of negative delays and exponential memories. We summarize our findings and conclude in Sec. 7. Most of the technical details are given in Appendix A -Appendix C. In Appendix A, we review the known results concerning the solution of systems of LDDEs. In Appendix B, we show how to generalize these results for the solution of systems of linear SDDEs. Finally, we apply the obtained results in Appendix C for the calculation of the time-correlation matrix and the covariance matrix for systems of linear SDDEs.

Stochastic Dynamics
We consider a two-dimensional system of N overdamped Brownian particles coupled via time-delayed pair interactions via the harmonic potential In panel a), we show stochastic trajectories of N = 18 Brownian particles bound by delayed harmonic forces. At long times and for small enough delays, the particles form molecular-like vibrating structures, while for large delays, they exhibit exponentially diverging oscilations. In b), we depict the model as a system of N = 3 Brownian particles interconnected by ideal springs with stiffness k and equilibrium lengths R, whose response to deformation is time-delayed by evaluating r ij at an ealier time t − τ .
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 − τ ) = |r i (t − τ ) − r j (t − τ )| is the distance between the particles i and j located at positions r i and r j . Clearly, the naive picture of springs representing the interactions poses a valid description only for a vanishing time delay τ .
Thus, the particles diffuse in the composed potential where the summation runs over all pairs (i, j). The ith particle is driven by the force x i and y i denoting the Cartesian coordinates of the position vector r i . In total, the N particle system obeys the set of nonlinear delayed Langevin equationṡ where the unit vector e ij = r ij /r ij points from particle j to particle i and the diffusion coefficient D 0 = (βγ) −1 is related to the inverse temperature β = 1/k B T and the friction coefficient γ via the Einstein relation (k B denotes the Boltzmann constant). The vectors η i comprise independent Gaussian white noises satisfying the relations The numbers α and β label the components of the vector η i , while i refers to the specific particle. The noise and the friction in Eq. (3) are related by the fluctuation-dissipation theorem [59] for a vanishing delay τ = 0, so that for τ > 0 the system is always out of thermodynamic equilibrium [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. Within the equilibrium model that obeys the fluctuationdissipation theorem, the stationary structures would simply be described by the Boltzmann distribution exp(−βV )/Z with potential V , inverse temperature β, and normalization Z. However, the physical situation at hand, where the delay is interpreted as a result of a feedback control mechanism and thus it is independent of the noise, requires the more involved description with Eq. (4) that leads to non-trivial nonequilibrium steady states.
The stationary structures attained 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 oscillations on its return to the energy minimum. 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. A similar system with a quasi-constant force between the particles (constant upto a change of sign at distance R) was discussed earlier in [11]. The benefit of assuming the harmonic potential is that it allows an extensive analytical treatment. 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.

Center of Mass
Similarly as for the dynamics considered in Ref. [11], the center of mass coordinate r c ≡ ( n i=1 r i ) /N of the system obeys the Langevin equatioṅ where η c ≡ N i=1 η i / √ 2 denotes Gaussian white noise satisfying Eq. (4) (with the labels i, j replaced by c) and the diffusion coefficient D c = D/N . Regardless of the interactions, the center of mass performs ordinary Brownian motion and, assuming the center of mass is at time t = 0 located at the point r 0 , the probability density function (PDF) for r c reads

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: where we have used the abbreviations e(t) ≡ e 12 (t) and r(t) ≡ 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 r = r 12 = r 1 − r 2 which obeys the equation of motioṅ with frequency ω ≡ 2k/γ (10) and η r ≡ (η 1 − η 2 ) / √ 2 describing Gaussian white noise satisfying Eq. (4) with the vector components i, j replaced by the label r.
Projecting Eq. (9) on the direction of the bond at time t [by multiplication with e(t) = (cos ϕ(t), sin ϕ(t))] and on the direction perpendicular to the bond [by multiplying it with e ϕ = (− sin ϕ(t), cos ϕ(t))], we obtain the equationṡ where ϕ(t, t−τ ) = ϕ(t)−ϕ(t−τ ) denotes the change of orientation of the vector e(t−τ ) during time τ . Above, we used the formulas r ≡ re,ṙ ≡ṙe + rφe ϕ and η r ≡ η r r e + η ϕ r e ϕ . From symmetry considerations, it follows that the stationary PDF for the orientation must be constant in ϕ. To gain analytical insight into the dynamics and PDF of the bond-length r, we linearize the coupled Langevin equations (11) and (12). If the angle dependent stiffness 2k/γ cos [ϕ(t, t − τ )] in Eq. (11) is strong enough such that r(t) can be approximated by R and the term proportional to [r(t − τ ) − R]/r(t) can safely be neglected independently of the noise strength, the formula (12) for the angle assumes the formφ The corresponding transition PDF (Green's function) for change of the orientation by Using this function, we average Eq. (11) over ϕ(t, t − τ ) obtaininġ where ω τ = ω exp(−2Dτ /R 2 ) is the natural relaxation rate. Note that the same formula with ω τ substituted by ω is obtained by simply assuming that the change of the orientation ϕ(t, t − τ ) of the bond per delay time τ is small, i.e. for 2τ D/R 2 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. (14) in more detail around Eqs. (24) and (25) below.
Equation (14) is a linear stochastic delay differential equation which can be solved analytically for r ∈ (−∞, ∞). In our setting, r ≥ 0 and thus we should solve Eq. (14) with the reflecting boundary at the origin. However, since we have assumed that |r(t − τ ) − R|/r(t) 1, we already work in the regime where r only seldom significantly deviates from R and thus the solution of Eq. (14) on the full real axis should approximate well the desired solution on the positive half-line. The solution of Eq. (14) for r ∈ (−∞, ∞) and t ≥ 0 in terms of deviation of the bond length from its equilibrium length, can be derived by several methods. We review two of them (time-convolutionless transform and Gaussian ansatz) for a general multidimensional linear SDDE in Appendix A-Appendix 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. (14) for r ∈ (−∞, ∞) and t ≥ 0 reads where the Green's function solves Eq. (14) with D = 0, x(t) = 0 for t < 0 and x(0) = 1. The symbol θ(x) in Eq. (17) stands for the Heaviside function. For an arbitrary initial condition x(t), t < 0 and x(0) = x 0 , the expression x 0 λ(t) in Eq. (16) must be substituted by Based on the value of the dimensionless delay ω τ τ which is a measure for the relevance of the delay relative to the natural relaxation time the Green's function λ(t) for the linear SDDE (16) exhibits three different dynamical regimes discussed in detail in Appendix A, in Fig. A1 and also below: (i) monotonic exponential decay for short delays 0 ≤ ω τ τ ≤ 1/e, (ii) oscillatory exponential decay for intermediate ('resonant') delays 1/e ≤ ω τ τ ≤ π/2, and (iii) oscillatory exponential divergence for long delays ω τ τ > π/2. The stochastic process x(t) in Eq. (16) arises as a linear combination of white noises and thus the corresponding PDFs must be Gaussian. Indeed, we find that one-and twotime conditional PDFs for x(t) with the initial condition δ(x) for t < 0 and δ(x − x 0 ) at t = 0 read where t ≥ t ≥ 0 and denote the mean of the shifted bond-length (15), its variance, and normalized time correlation, respectively. The function P 1 (x, t|x 0 , 0)dx stands for the probability that the system which departs with certainty from state x 0 at time 0 is found at time t somewhere in the interval (x, x + dx). Similarly, P 2 (x, t; x , t |x 0 , 0)dxdx denotes the probability that the (shifted) bond length is in the interval (x , x + dx ) at time t and at a later time t in (x, x + dx) 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 τ = 0 [63]. The non-Markovianity of the process (14) 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 (18) and (19) are given by Eqs. (B.6) and (B.7) in Appendix B, respectively. Interestingly enough, the diffusion and drift terms in the FPEs are given by the natural scales 2D and ω τ only in the limit τ → 0. Moreover, both coefficients acquire a time dependence, determined by the function λ(t). Specifically, the diffusion and drift coefficients in Eq. (B.6) for P 1 read D τ (t) = Dλ 2 (t)d t 0 dsλ 2 (s)/λ 2 (t) /dt and ω τ (t) = −λ(t)/λ(t), respectively ‡. While the drift coefficient in Eq. (B.7) for P 2 is also given by ω τ (t), the diffusion coefficient reads 2D τ (t) + 4Dλ(t) t 0 ds d [λ(t − s)λ(t − s)/λ(t)] /dt. This difference in diffusion ‡ The case of λ(t) = 0, where these coefficients diverge, is discussed in more detail in Sec. 5.1. coefficients is the reason why the PDF P 2 can be constructed from P 1 in the standard way for Markov processes only for τ = 0, where both diffusion coefficients coincide. For τ = 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, ∞; x 0 , 0) ∝ exp (−ωx 2 /2D). For a nonzero delay in the regimes (i) and (ii), i.e., when the noiseless solution (17) and thus the mean value µ(t) = x 0 λ(t) of x converges to 0 for t → ∞, the system relaxes into a time-independent non-equilibrium steady state Fig. 8. This state is characterized by a nonzero entropy production rate [56]. For long time delays, no stationary state exists. In comparison to the analogous setting with a piece-wise constant force [11], this destabilization for long delay times τ is a new feature.
In the regimes (i) and (ii), the variance ν(t) = x 2 (t) − x(t) 2 converges to the stationary value [34,30,31] The derivation of this formula is given in Appendix C, where we also derive an analytical expression for the stationary time correlation function C(t) = lim s→∞ x(s)x(s + t) .
Note that the variance (23) diverges entering the unstable regime (iii) for ω τ τ → π/2 upon. The formula (23) finally allows us to specify the regime ν ss R 2 where the approximations [r(t − τ ) − R] /r(t) ≈ 0 and r(t) ≈ R used in the derivation of Eq. (14) from Eqs. (11) and (12) are reasonable because the PDF for r is relatively sharply peaked around the mean bond length R. As we already noted, Eq. (14) is also valid in the case 2τ D/R 2 1 when the bond rotates only slightly in each delay interval and thus the angle ϕ(t, t − τ ) in Eqs. (11) and (12) is small. However, in this case we need to additionally assume that ν ss < R 2 in order to ensure that the error caused by considering the wrong boundary condition at r = −R is small. Altogether, the used approximation is expected to give good results if one of the two conditions 2τ D/R 2 1, and ν ss < R 2 , is fulfilled. An example of the stochastic evolution of the dimer obtained from BD simulations of the exact system (11) and (12) 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 (20) for different values of the equilibrium length R. As expected, the approximate analytical formula (20) describes well the exact result for large enough R satisfying at least one of the inequalities (24) and (25). For smaller values of R, the analytical result underestimates the correct value. This is because the bond length in the BD simulation is obtained from Eq. (11) 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 (21) for the bond length variance ν(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 (23) as the time delay τ reaches the unstable regime (iii). This means that the delay tends to delocalize structures. On the other hand, the variance can be optimized as a function of the frequency ω τ . The best localized structure is obtained for the frequency fulfilling the formula cos(ω τ τ ) = ω τ τ , i.e. ω τ ≈ 0.74/τ and thus 2kτ /γ ≈ 0.74 exp(2Dτ /R 2 ). For the corresponding value 0.74 of the dimensionless delay ω τ τ the system is in the dynamical regime (ii) with exponentially decaying oscillations.

Trimer
Let us now consider the system composed of three particles. Then, Eq. (3) gives the system of six coupled equations of motion: For the relative coordinates r 12 (t) = r 1 (t) − r 2 (t) we obtaiṅ where ω = 2k/γ and similarly forṙ 13 (t) andṙ 32 (t). To get analytical results for bond lengths r ij (t) = |r ij (t)|, we multiply the formulas forṙ ij (t) by the corresponding unit vectors e ij (t) = r ij (t)/|r ij (t)| and linearize the resulting equations. To this end, we need to deal with the expressions e i (t − τ ) · e j (t), i, j = 1, . . . , 3, where we introduced a shorthand indexing 1 ≡ 12, 2 ≡ 13, 3 ≡ 32. For a vanishing delay τ = 0, e i · e i = 1 and the scalar products e i · e j describe the angles of the triangle formed by the three particles (see Fig. 1). We find that up to the leading order in the equilibrium bond length R the triangle is equilateral and thus the internal angles are π/3, leading to the relations e 1 · e 2 = e 1 · e 3 = −e 2 · e 3 = 1/2 + O(1/R) with a correction that is on the order of 1/R. The linearized equation for the relative coordinate x i ≡ r i − R thus readṡ where the lower index i ≡ i mod 3 is considered as periodic with the period 3, i.e., Similarly as in the case of the dimer, Eq. (27) describes the dynamics of x i (t) well for large equilibrium bond lengths R and for time delays small compared to reorientation times of the unit vectors e i . For an analytical treatment, it is advantageous to rewrite the system (27) in the matrix formẋ for the three-dimensional column vector . In Eq. (28), the noise vector ξ(t) is given by ξ(t) ≡ A 1 η 1 (t) + A 2 η 2 (t) with the auxiliary noise vectors η j (t) ≡ [η j 1 , η j 2 , η j 3 ] , j = 1, 2, containing the jth components of the original noises η i (t), i = 1, . . . , 3. From the system (27) follows that the matrices M , A 1 and A 2 read where e j i denote jth component of the two-dimensional unit vector e i . The timecorrelations between the three components of the noise vector ξ(t) are not mutually independent and read Due to the linearity of Eq. (28) and Gaussianity of the noise, the Green's function for the one-time PDF P 1 (r, t|r 0 , 0) is Gaussian [37], determined the mean value µ(t) = x(t) and the covariance matrix . We review in detail the derivation of these functions in Appendix A and Appendix B.2.
For the initial condition x(t) = 0, t < 0 and x(0) = x 0 we get where λ(t) denotes the Green's function for Eq. (28) given by Multiplying λ(t) by the vector [1, 1, 1] , using the formula M [1, 1, 1] = 3/2[1, 1, 1] and comparing the result to the one-dimensional Green's function (17), we find that λ(t) undergoes with increasing t (i) monotonous exponential decay to 0 for 0 < 3ωτ /2 ≤ 1/e, (ii) oscillatory exponential decay to 0 for 1/e < 3ωτ /2 < π/2 and (iii) oscillatory exponential divergence for π/2 < 3ωτ /2. In the regimes (i) and (ii) the stationary value of the covariance matrix reads where I denotes the identity matrix. This formula follows from the results of Appendix C after substituting the matrices ω and σσ from the formula (C.7) in the appendix by ωM and 4DM . In the appendix, we also derive an analytical expression for the stationary time correlation matrix C(t) = lim s→∞ x(s)x (s + t) . The regime of stability 3ωτ /2 < π/2 of the trimer can also be determined from the condition that the matrix cos (ωτ M ) is not singular, i.e., its determinant cos 2 (3ωτ /2) cos (3ωτ ) is nonzero. Due to the symmetry of the problem, all diagonal elements of the matrix 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 full 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 with the frequency ω and with the delay τ is similar to the behavior of the variance (23) for the dimer. Specifically, the diagonal elements of K ss monotonously increase (the PDF for the bond lengths become broader) with τ and exhibit a minimum as functions of ω, opening a possibility to optimize the width of the bond length PDF. The off-diagonal elements of K ss monotonously increase (the individual bonds of the trimer become more correlated) both with the delay and with the frequency.

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 µ(∞) = R, respectively, and (iii) an exponential divergence. The performed BD simulations confirmed that for large R, when the full model is well described by the approximate analytical formulas, these regimes can indeed be observed also in the full model (3). Furthermore, the analytical study predicted that the dimer is in the unstable regime (iii) for ωτ > π/2 and the trimer for ωτ > π/3. Let us now discuss how general the presented findings are.
The stationary average bond length µ can be determined by minimizing the potential energy V = 1 2 (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 µ = R is accessible only for the dimer (N = 2) and the trimer (N = 3). For larger molecules, the average bond length decreases as a result of the infinite range of the potential. The system asymptotically  relaxes to the depicted structures if the noise D vanishes and the reduced delay time ωτ is small enough such that the system is in the dynamical regime (i) or (ii). Nonzero D leads to fluctuations around the asymptotic structures and large ωτ causes exponentially diverging oscillations.
We have solved the full 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 ωτ . More precisely, all the curves in Fig. 4 c) are plotted using the same parameters except for N and, while the curves for N ≤ 3 are in the regime (ii), the curves for N > 3 correspond to the regime (iii). These observations are in accord with our analytical findings for the dimer and trimer.
By analyzing the mean bond length at late times, we have evaluated the critical dimensionless delay (ωτ ) determining the threshold between the regimes (ii) and (iii) for N = 2, ..., 10. In Fig. 5, we show the its rescaled value where the coefficient 2/π is introduced for the comparison to the approximate result for the dimer. We find that the stability factor is well described by C/N as suggested by the approximate analytical results for the dimer [(ωτ ) = π/2] and trimer [(ωτ ) = π/3]. However, the analytical results would imply that the constant C equals to 2, which is smaller than the value C ≈ 3 obtained from the full model. The actual dimer and trimer are thus more stable than their linearized versions considered in our analytical study.
To better understand this behavior, 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, the particle does not feel any force in the time interval t ∈ [0, τ ] such that x(t) = x 0 for all t ∈ [0, τ ]. In the subsequent time interval t ∈ [τ, 2τ ], the particle experiences the force F (t) = −ω τ x 0 pushing it towards the opposite wall of the trap. For times t > 2τ the force starts changing dynamically according to the earlier position at time t − τ . The particle keeps its direction of motion until it reaches the position x 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. The described situation with the harmonic potential corresponds to the approximate analytical model considered in Sec. 2.2.
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. The potential is now realized to contain a more extended low-energy region compared to the purely harmonic case, which effectively leads to a decreased dimensionless delay ωτ and smaller amplitude of the oscillations. 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 τ , while in the harmonic model these oscillations appear only if ωτ is large enough. Second, the amplitude of the oscillations in the quasi-constant velocity model is always constant in time, while the osculations in the harmonic model either vanish with time, if the system is in the regime (ii), or explode with time, if the system is in the regime (iii). The behavior observed in the harmonic model can be traced back to the increase of the force with the particle distance and thus we expect an analogous behavior also for other systems with time-delayed forces increasing with the distance.

Entropy Fluxes
The investigated model, much as the model discussed in Ref. [11], is inspired by selforganized 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 the feedback, i.e., without the time-delayed harmonic interactions (2), the particles would spread diffusively and the system entropy would increase accordingly. The feedback control utilizes the information about the particle positions to drive the system into a non-equilibrium steady state with a time independent PDF P (x) and thus with a time constant configurational entropy S = −k B dxP (x) log P (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 −Ṡ − 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. Evaluating the stationary entropy productioṅ S F due to the feedback control mechanism and the stationary entropy productionṠ D due to breaking the fluctuation-dissipation theorem in Eqs. (3) and (4) for τ > 0, one can define the feedback efficiency as the ratio η F = −Ṡ − /(Ṡ F +Ṡ D ). The entropy productioṅ S D can be calculated along the lines of Ref. [56]. The entropy productionṠ F = q H /T is the housekeeping heat flux q H flowing to the bath at temperature T , due to the overall operation of the feedback device, divided by the bath temperature. It clearly depends on the specific technical realization of the feedback. In all known relevant realizations of the feedback in microscopic systems [11,13,16,18,19], the housekeeping heat flux is very large compared to the "functional" energy fluxes in the controlled system, resulting in a largeṠ F compared to −Ṡ − , so that the efficiency η F of such devices is usually negligibly small.
To evaluate the entropy flow due to the feedback (the time-delayed harmonic interaction) in the present setup, we proceed along the similar lines as in Refs. [11,64]. 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 (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 Shanon entropy S(t) = −k B dxP (x, t) log P (x, t) can formally be written aṡ whereṠ + (t) stands for the positive entropy flowing into the system due to the diffusive spreading of the particles andṠ − (t) corresponds to the outflow of entropy due to the feedback.
Assuming that the stochastic dynamics of the column vector x(t) describing the bonds obeys the generalized Langevin equationẋ(t) = F[x(t), x(t − τ )] + ση(t), where η denotes a zero mean Gaussian white noise with the covariance matrix η i (t)η j (t ) = δ ij δ(t − t ), the dynamical equation for P (x, t) can be written in the form [30,65] In this equation, the fist term on the right stands for the diffusive spreading of the PDF. The term L[x, t] corresponds to the time-delayed force F[x(t − τ )] in the Langevin equation and thus it describes the effect of the feedback. Its concrete form is not relevant for the discussion below and thus we refer to the works [30,65] for more details about its structure.
Inserting Eq. (37) into the formal time derivativeṠ(t) = −k B dx[∂ t P (x, t)] log P (x, t) of the system entropy S(t), we find thaṫ The last equation allows us to calculate the amount of entropy taken out of the system due to the feedback per unit time,Ṡ − (t), from the PDF P (x, t) without knowing the explicit form of the operator L.
Let us now evaluate the three entropy fluxes (36), (38) and (39) for a general d-dimensional Gaussian PDF where |K(t)| denotes determinant of the covariance matrix K(t). The corresponding entropy S(t) reads leading to the rate of changeṠ From Eq. (38) we then find thaṫ For finite times, all the entropy fluxes depend on the initial conditions and can be determined from Eqs. (39), (42) and (43).
Let us now focus on the specific system considered in this paper. Using our analytical findings for the dimer with σσ = 4D and K(t) = ν(t), we find that where λ(t) is given by Eq. (17) and the variance reads ν(t) = 4D t 0 ds λ 2 (s). Similarly, for the trimer we obtain using σσ = 4DM and Jacobi's formula for the derivative of determinants, the expressions If not specified otherwise we used the parameters D = 1 µm 2 /s, τ = 1 s and ω = 1 1/s. For the dimer, we used the approximation ω τ = ω, which is accurate for 'large' equilibrium spring length R.
where λ(t) is given by Eq. (33) and the covariance matrix K(t) reads 4DM t 0 dsλ 2 (s). The formulas (44) - (47) are valid both in the stable regimes (i) and (ii), where the system at long times relaxes to a stationary time-independent structured state, and in the unstable regime (iii). In the unstable regime, the variance ν(t) and the covariance K(t) diverge in time. As a result, the system entropy S(t) diverges and the entropy floẇ 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Ṡ(t) and thus also the entropy outflowṠ − (t) remain finite oscillating functions, as can be seen for the dimer by using the exponential long-time approximation (A.7) for λ(t), and similarly for the trimer.
For the purpose of structure formation, only the regimes (i) and (ii) are of interest, because only then the PDF reaches a time-independent non-equilibrium steady state at long times, i.e. lim t→∞ S(t) = S, lim t→∞Ṡ (t) = 0 andṠ − ≡ lim t→∞Ṡ− (t) = − lim t→∞Ṡ+ (t). Let us now evaluate the long-time stationary system entropies S and entropy fluxesṠ − (t) maintaining the molecular-like structures formed in our model for the dimer and the trimer. Using the asymptotic formulas (23) and (34) for the dimer bond-length stationary variance ν ss and trimer covariance K ss , we find from Eqs. (44) - for the dimer and for the trimer. In the last two formulas, a = K ss (1, 1) denotes diagonal and b = K ss (1, 2) off-diagonal elements of the covariance matrix. The system entropies (48) and (50) are determined by the width of the PDFs for the bonds. Therefore, they monotonously increase with temperature T = γD/k B and with the delay time τ , and exhibit a minimum as functions of the frequencies ω τ (dimer) and ω (trimer), similarly as the variance ν ss and the diagonal matrix elements of the covariance matrix K 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 (49) and (51) 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 → 0 because they must inevitably vanish for zero noise, where the PDF for the system evolution is a δ-function for all times. We plot the stationary entropy fluxes (49) and (51) as functions of the delay τ in Fig. 6 a) and frequency ω 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 |Ṡ − |, depicted in Fig. 6 b), arises as a result of a competition between stronger confinement, corresponding to larger frequencies ω, and gradual destabilization with increasing ωτ , when the system enters the unstable regime (iii). For the figure, we used for simplicity the approximation ω τ ≈ ω for the dimer.

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. While for the purely deterministic motion these transitions only appear for time delays τ in the unstable regime (iii), in a system affected by thermal fluctuations the transitions occur for arbitrary delays. Evaluation of frequencies of such transitions, which measure the stability of the individual isomers, is the main topic of transition rate theory [23].
The transition rate κ A→B (t) for switching from a conformation A to a conformation B at time t can be (for arbitrary dynamics) found from the mean number of transitions N A→B (t) from A to B during an infinitesimal time interval ∆t as κ A→B (t) = N A→B (t)/∆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 can easily be evaluated in simulations, this can take a 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→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 κ A→B (t) =Ṡ A→B (t)/S A→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→∞ κ A→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 large energy barrier compared to the thermal energy. Kramers' theory was later 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 [66] who derived a dynamical correction to Kramers' result. While their analysis was based on a parabolic barrier, Pollak [67] 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 κ A→B (t) for driven overdamped systems can be calculated using the recent theory of Bullerjahn et. al [24] for 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. [68]. Based on their small-delay approximation, Guillouzic et al. [69] 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. [70] 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 enough for small and moderate delays only.

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 different isomers in two dimensions). The setting is described by the approximate Langevin equation (14) for the inter-particle distance r = |r 1 − 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. This redefined signed bond lengthr then diffuses in the cusped double-well potential depicted in Fig. (7), with the diffusion coefficient 2D τ (t). 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 =r − R > 0. Based on the approximate solution (16) to the Langevin equation (14) assuming r ∈ (−∞, ∞) and x(t) = 0 for t < 0, we have found that the one-time PDF P 1 = P 1 (x, t) = P 1 (x, t|x 0 , 0) obeys the FPE (B.6), which reads This equation describes diffusion in the harmonic potential γω τ (t)x 2 /2 with the time- 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 Markovian dynamics for calculating the transition rate κ(t) for switching between the isomers. In the Markovian case, the transition rate for surpassing the energy barrier at r = 0 to the other isomer can be calculated from Eq. (53) with absorbing boundary at x = −R [63]. Let us now review several methods for the determination of transition rates for Markovian dynamics described by Eq. (53) and compare the results to BD simulations of the full model.

Numerical Method
We first consider the situation when the system dwells in the state x(t) = 0 for t ≤ 0 and then it starts to diffuse in the time-dependent potential (52). Then, the time-dependent Markovian rate κ M (t) can be determined from the equation for the normalized PDFP a (x, t) = P a (x, t)/S a (t) for position of the particle surviving in the right well of the cusped potential [71]. Here, P a (x, t) is the solution to the FPE (53) with absorbing boundary at x = −R and S a (t) = ∞ −R dxP a (x, t) is the probability that the particle survives in the right well until time t. Eq. (54) follows from Eq. (53) immediately after using the definitions of the PDFP a and of the transition rate The problem (54) can be solved numerically and we did that using the method presented in Ref. [72].
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 (19) of the FPE (53) with natural boundary conditions. Specifically, one approximates the probability current across the absorbing boundary by and the survival probability S a (t) by In the last expression, the symbol Erf denotes the error function and the variance ν(t) is given by Eq. (21). The approximate Markovian transition rate is then given by In Fig. 8, we show the frequency ω τ (t), the effective diffusion coefficient 2D τ (t), survival probabilities S a (t) and S 1 (t) and the transition rates κ M (t) and κ B (t) as functions of time t for parameters in the dynamical regime (i). One can observe that both the parameters ω τ (t) =λ(t)/λ(t) and D τ (t) = 2Dλ(t) 2 + ω τ (t)ν(t) and the rates saturate with time. The relaxation time of ω τ (t) is determined by the time in which the Green's function λ(t) approaches the long-time exponential representation (A.7), and the corresponding stationary value, ω τ (∞) = 1/t R , is controlled by the relaxation time t R for decay of λ(t) to 0, see also Fig. A1. The effective diffusion coefficient converges to the value D τ (∞) = 2D[1 + sin(ω τ τ )]/ cos(ω τ τ ) ≥ 2D, determined by the stationary variance ν(∞) = ν ss , see Eq. (23). The transition rates relax with the relaxation time t R , similarly as the corresponding PDFs P a and P 1 . The analytical expressions for S 1 (t) and κ B (t) approximate the numerical results for S a (t) and κ M (t) best for short times t t R , where the PDFs P 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 κ B overestimates the corresponding exact rate κ M , see also Figs. 9 a) and b). For large time delays, the effective barrier height over the effective thermal energy decreases so that the assumptions of the transition state theory are not valid any more, and κ B < κ M , see Fig. 9 c).
The situation of small barriers compared to thermal energy can be understood from the behavior at vanishing potential strength ω → 0, when D τ (t) = 2D and the finite time transition rates κ M (t) and κ B (t) can be calculated analytically. Namely, the PDFsP a and P 1 in the definitions (55) and (59) [63] leading to the formulas Since the denominator in the expression for the rate κ M (t) is smaller than that for κ B (t), we conclude that, for weak energy barriers, the inequality κ B (t) ≤ κ M (t) holds. On the other hand, for very large energy barriers, the PDFsP a and P 1 are (almost) identical because the absorbing barrier at x = −R is effectively inaccessible. In such a case, the probability current j(−R, t) is (almost) zero, S(t) ≈ 1, and j * ≈ j − 2D∂ x P 1 | x=−R < 0 leads to the inequality κ B (t) ≥ κ M (t).
To gain a deeper insight into the behavior of the transition rates, let us consider the stationary regime t → ∞. In this regime,P a (−R, ∞) = 0 due to the absorbing boundary condition at x = −R, and where ∂ t ≡ ∂/∂t and ∂ x ≡ ∂/∂x. Then the transition rates κ M (∞) = −j Da and κ B (∞) = −j D1 are determined solely by the diffusive fluxes The smaller the frequency ω, the wider the PDFs P 1 andP a . For small ω, the boundary at x = −R is in the region where the PDF P 1 has its maximum and it is also close to the maximum ofP a . In such a case, the PDFP a , that must vanish at x = −R, changes near the boundary faster thanP 1 , leading to |j D1 | < |j Da | and κ B (∞) ≤ κ M (∞), in accord with the previous paragraph. With increasing ω, the boundary shifts away from maximums ofP 1 andP a . Due to the trajectories trapped in the absorbing state [71,73], the maximum ofP a is slightly farther away from the absorbing boundary than the maximum ofP 1 , and thus, with increasing ω, the tail ofP a with small derivative (small j Da ) reaches the boundary at x = −R before the tail ofP 1 . Hence, for large enough barrier height with respect to the noise strength, the inequality between the rates turns over to κ B (∞) ≥ κ M (∞). Finally, for very stiff traps (ω → ∞), both j Da and j D1 vanish and κ M (∞) = κ B (∞) = 0.
As shown in the figure 8 d), the transition rates converge with time to constant values in the regime (i), where the limits lim t→∞ ω τ (t) and lim t→∞ D τ (t) exist. However, also in the regime (ii), the PDF P 1 assumes, after long times, the time-independent stationary form P 1 = exp (−x 2 /2ν ss ) / √ 2πν ss with the variance ν ss given by Eq. (23). This suggests that also in this regime the transition rate should saturate at long times. However, both the frequency ω τ (t) and the diffusion coefficient 2D τ (t) actually exhibit diverging oscillations caused by the oscillations in the Green's function λ(t) (17), in the regime (ii). These divergences induce problems both in the FPE and in the approximate calculation of the rate using Bullerjahn's method and imply that the Markov description can not be valid (at least) in the dynamical regime (ii). Nevertheless, let us now investigate to what extent the long-time transition rate κ = κ(∞) obtained from BD simulations of the Langevin equation (14) is captured by the predictions (55) and (59) above.

Long-Time Behavior and Kramers' Method
Assuming that at long times the PDFP is time-independent and the limits lim t→∞ ω τ (t) and lim t→∞ D τ (t) exist, we can rewrite the formula (54) as the eigenvalue problem for the long time Markovian transition rate κ M = κ M (∞). We solve this formula numerically using the method described in Refs. [71,72]. The steady state value of the transition rate predicted using Bullerjahn's method reads with the survival probability S 1 (∞) = 1 + Erf R/ √ 2ν ss /2 and the probability current j * (∞) = ω τ (∞)R exp (−R 2 /2ν ss ) / √ 2πν ss . For large barriers, where S 1 (∞) ≈ 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 b given the inverse thermal energy is β, In the penultimate equality, we used the appropriate inverse thermal energy β = 1/2γD τ (∞) and barrier height E b = γω τ (∞)R 2 /2, and also the asymptotic form of the diffusion coefficient D τ (∞) = ω τ (∞)ν ss , which follows from the condition lim t→∞ λ(t) = 0 valid in the dynamical regimes (i) and (ii).
Renormalized Transition Rates Interestingly, the term ω τ (t), which causes divergences of the diffusion coefficient and the frequency of the potential in the FPE (53), does not enter the argument of the exponential in the rates κ B and κ K . This means that it just determines the kinetic prefactor, as can be also observed directly from the long-time form ∂ t P 1 = ω τ (t)∂ x (x + 2ν ss ∂ x ) P 1 of the FPE (53). In the dynamical regime (ii), the kinetic prefactor in Eq. (64) cannot be correct, due to the diverging oscillations in the time-dependent frequency ω τ (t). Nevertheless, the exponential term seems to be reasonable, and thus it is tempting to use in the prefactor of the transition rates simply ω, instead of the problematic ω τ (∞). This substitution gives the correct pre-exponential factor of the rate for vanishing delay τ = 0, where ω τ (∞) = ω and 2D τ (∞) = 2D. We denote the rates with the renormalized prefactor as with x = M, B or K indicating Markov, Bullerjahn, or Kramers, respectively. The necessity to change the kinetic prefactor in the rates stems from the fact that, although the absorbing boundary condition we used in Eq. (62) is correct for Markov dynamics (τ = 0), it can not be precisely valid for the time-delayed dynamics (τ > 0). To see this, it is enough to realize that the delayed system arriving at the boundary at a time t does not feel the energy barrier E b = γωR 2 /2, but the barrier with energy γω[r(t − τ )] 2 /2.
In Figs. 9 a) -c), we compare the various analytical predictions κ x andκ x , x = M, B or K, for the long-time transition rate (lines) with the asymptotic rate κ (symbols) calculated from BD simulations of Eq. (14) using the inverse first passage time for reaching the absorbing boundary at r = 0. For a broad range of parameters fulfilling k B T V (−R, 0), we have found that the rate κ can be predicted reasonably well only for values of τ in the dynamical regime (i) (ωτ < 1/e ≈ 0.37) and in the first part of the dynamical regime (ii) (ωτ < π/2 ≈ 1.57). The rate κ is best approximated by the expressionκ M obtained numerically from Eq. (62) with ω substituted for ω τ (∞) in the operator L. From the analytical expressions works best the re-scaled Bullerjahn expressionκ B [see Eqs. (63) and (65)]. However, in the figure we have used parameters leading to a large barrier E b and thus Kramer's and Bullerjahn's predictions, κ K and κ B , almost coincide. As a consequence, the line forκ K in the figure overlaps withκ B , similarly as the line κ K (that was suppressed in the figure for better readability) overlaps with κ B .
Delay-Dependence of κ In Fig. 9 d), we also show the behavior of the rate κ in the parameter regime (iii) inaccessible for the analytical and numerical formulas due to the diverging oscillation in the system dynamics. The simulated transition rate in Figs. 9 a) -d) is first approximately exponential and thus convex in τ , then its curvature changes to concave and it runs through a maximum and, finally, the rate starts to decrease. The value of τ where the curvature changes sign coincides with the boundary π/2ωτ ≈ 1.57 between the dynamical regimes (ii) and (iii). Interestingly, no qualitative change of κ(τ ) is observed at the boundary ωτ = 1/e ≈ 0.37 between the regimes (i) and (ii). It is tempting to attribute, the (approximate) exponential increase of the rate with τ in regimes (i) and (ii) to the increase of the steady-state variance ν ss , which is given as the ratio of the effective energy barrier V (−R, ∞) and the effective thermal energy γ/2D τ (∞). Although this explanation may work well for small delays, it breaks down for values of τ in the second half of the regime (ii), where the actual rate κ is no longer well approximated by our analytical and numerical predictions. This means that the identification of the parameters γω τ (t) and 2D τ (t) with the effective potential stiffness Figure 10. Transition from the clockwise to the counterclockwise isomer of a trimer molecule. To calculate the transition rate from the clockwise to the counterclockwise isomer, we choose the coordinate frame such that x-axis points from particle 1 to particle 3. By this choice of frame, all transitions are mapped to the crossing of the x-axis by particle 2. The grey region of width ∆x depicts a transition zone or "neutral region" that helps to avoid that small fluctuations around the axis are counted as isomer transitions. and the effective diffusion coefficient, respectively, suggested by the effective Markov equation (53), is reasonable for relatively small values of τ only. In the dynamical regime (iii), the particle undergoes oscillations with an amplitude that increases both with τ and t. The corresponding transition rate, obtained form the BD simulations, thus decreases with the delay τ as a result of oscillations 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 away from the boundary increases with t (so that the transition rate decreases with t), where t is the duration of the simulation.

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, in two dimensions, the trimer molecule can form two different isomers. 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 change to 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 κ = κ(∞) for the isomer transition in the steady state by means of Kramers' theory and compare it to BD simulations. In the BD simulations, we have evaluated the rate κ using the angle ϕ between the abscissas |12| and |13| and a neutral region of width ∆x as exemplarily shown in Fig. 10. A neutral state |0 is introduced to avoid over counting due to fluctuations of ϕ around 0 and π . It is occupied if the smallest height of the triangle formed by the three particles, is smaller ∆x/2. i.e. either if |ϕ| < φ > ≡ max [arcsin(∆x/2r 12 ), arcsin(∆x/2r 13 )] or |ϕ| > φ < ≡ arccos(∆x/2r 12 ) + arccos(∆x/2r 13 ). For ϕ ∈ [φ > , φ < ] the system is said to be in the clockwise state |1 , while for ϕ ∈ [−φ > , −φ < ] it is in the counterclockwise state |−1 . To calculate the transition rate, we have counted the number of transitions between the states |±1 during a specific simulation time window, where the transition occurred if the system underwent the sequence of states |1 → |0 → |−1 or |−1 → |0 → |1 .
The resulting transition rate κ as a function of the time delay τ 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 κ(τ ) 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. Concretely, we only focus on the distance r b of the particle 2 from the x-axis and construct an appropriate effective energy barrier E b,eff and diffusion coefficient D eff . In the steady state, the three particles are most likely found in vertices of an equilateral triangle. Using the notation of Sec. 2.3, we express the vector r b from the particle 2 to the center of the abscissa |13| as r b = r 12 + r 31 /2. Hence the Gaussian white noise η b (t) corresponding to the coordinate r b is given by where we have used the relations (4). Based on these considerations, we approximate the effective diffusion coefficient by D eff = 3D.
For the corresponding effective energy barrier, we make the ansatz E b,eff = C(τ )γωR 2 /12. In order to compose this ansatz, we first considered the minimum value E b = kR 2 /6 = γω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). Then, we introduced into the resulting barrier energy E b a, possibly delay-dependent, unknown dimensionless factor C(τ ) in order to account for the time delay and all additional ways the particle 2 may take in the multidimensional energy profile in order to pass from |1 to |−1 . The resulting Kramers' rate for the transition from |1 to |−1 thus reads where a(τ ) is a further unknown prefactor that may depend on all parameters of the model [see, for example, Eq. (64) for Kramers' rate in a cusped potential].
In order to test the formula (67), we have fitted it to the transition rate κ obtained from BD simulations as a function of the diffusion coefficient D. The results for different values of the time delay are shown in Fig. 11 b) and the corresponding values of the coefficients C(τ ) and a(τ ), obtained from the fits, are given in Tab. 5.2. The presented results prove that the transition rate exhibits exponential increase with the diffusion coefficient for a relatively broad range of values of the time delay, and thus the D -dependence of κ can be relatively well described by the Kramers-type ansatz (64).
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 (67) to the BD data shown in Fig. 11 b).
the replacements ω → ω τ,ss and D → D τ,ss , did not lead to a significant change of the value for C(τ ). We thus conclude that the deviations of C(τ ) from unity are mainly caused by the assumption that the particle will almost in all cases cross the axis at the minimum of the potential energy. An analysis that would take this issue into account would have to employ a multidimensional transition rate theory, which is beyond the scope of the present paper.

Extensions to other memory kernels
So far, we have considered only the interactions involving a given positive time delay. Nevertheless, most of the presented results can be directly applied also to other models with delay or memory. To see this, note that the central equation (B.1) is equivalent to the generalized Langevin equation (GLE) with positive frequency ω > 0 and the memory kernel given by The simplest generalization of systems with the memory kernel φ(t) = δ(t − τ ) are systems with multiple different time delays with the memory kernel Properties of these systems are studied in Ref. [35].
An interesting problem in the context of active matter represents a negative time delay as the individual active agents may react to a future state of their neighborhood which they predict using its present state. For idealized systems capable of a perfect prediction, the GLE (or, equivalently, the linear SDDE) contains the memory kernel φ(t) = δ(t + τ ), τ > 0 and it can be solved using the strategy described in Appendix A. The resulting Green's function is the so called Bruwier series [74] that is convergent for |ω| < |eτ |. Alternatively, the series (70) can be written in the form [75] where s is the absolute value of the smallest root of the equation ρ = −ωe τ ρ . More realistic predictive systems might instead only have an imperfect prediction of the future position x pre (t + τ ) = x(t + τ ). Therefore, we reformulate the deterministic part of the GLE (68) asẋ (t) = −ωx pre (t + τ ).
One of the reasonable strategies for predicting x pre is to use the linear extrapolation which is identical to a small delay expansion of x(t + τ ). The equation of motion (72) then assumes the time local form which is solved by x(t) = x 0 λ(t) with the exponentially decaying Green's function The rescaled frequency ω r = ω/(1 + ωτ ) 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 positive time delays [τ < 0 in Eq. (72) ], the presented first order approximation predicts dynamics with rescaled frequency ω r increasing with increasing time delay (for ωτ < −1). However, such dynamics would lead to smaller variance ν ss for non-zero delays than for a vanishing delay, which contradicts our exact result (23), highlighting the limited practical significance of the naive small-delay expansion, as already pointed out in Ref. [30]. The most frequently used is the exponential memory kernel 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 where Ω = ω 2 − b 2 /4, which is reminiscent of the Green's function (A.7) for the hereby studied system with positive time delay [φ(t) = δ(t − τ )]. The difference is that, while the Green's function (77) always decays exponentially with time, the Green's function (A.7) allows also for negative relaxation times and thus exponential increase with time.
Properties of the Green's function (77) as well as those corresponding to a power-law memory are discussed in detail in Ref. [76].

Conclusion and outlook
Inspired by the ongoing debate on self-organized active matter and the experiments of Khadka et al. [11], we considered a 2d system of N Brownian particles interacting via time-delayed harmonic interactions, as depicted in Fig. 1. 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), which we find to 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. 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 these results with Brownian dynamics simulations of the full models (3). We have found good quantitative agreement in the parameter regimes where the systems evolve relatively close to the minimum energy configurations, 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τ , 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 exhibit exponentially diverging oscillations. Our analysis of systems with N ≥ 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τ than predicted from the linearized equations, the critical value of the product kτ decaying approximately as 1/N . Therefore, larger systems are more unstable than smaller ones. 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 considered in Ref. [11] with constant forces 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 from the system by the feedback per unit time 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. 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 τ . Concerning the analytical results, the best agreement with the true rates was obtained by the Bullerjahn's formula (63) with effective barrier height and diffusion coefficient taken from the time-local FPE (53) and the prefactor rescaled according to Eq. (65). 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.
In the last Sec. 6, we have demonstrated that most of the presented equations can be used also for systems with memory kernels different from that for discrete time delays, i.e. φ(t) = δ(t − τ ). It is enough to substitute the Green's function λ(t) (A.5) corresponding the delayed Langevin equation by the Green's function corresponding to the memory kernel of interest. We have reviewed prototypical memory kernels, discussed the differences and similarities of the corresponding Green's functions, and leave a more detailed study for future work.
As a further extension of our work, it would be interesting to consider physically more realistic interactions that would decay to zero for 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 local a perturbation through the system in case of a large number 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.

Appendix A. Solution of the Noiseless Problem
In this appendix we solve the multi-dimensional linear delay differential equation where ω is a positive semi-definite matrix with real entries and x(t) is a column vector. Laplace transformation of this equation leads to the formula . The solution of this equation for the Laplace transformed variable reads where 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τ /s k+1 is given by (t − τ ) k θ(t − τ )/k! [77] and thus the formula (A.3) can be inverted. Finally, the solution of Eq. (A.1) is given by is a matrix-valued function which solves Eq. (A.1) with the initial condition x(t) = 0 for all t < 0 and x(0) = I. In the present paper, we always assume that the system is initialized at time t at position x 0 witha special history, namely x(t) = 0 for all t < 0, allowing us to simplify Eq. (A.4) to x(t) = λ(t)x 0 . The only fixed point of the delay differential equation (A.1) is x(t) = 0. In order to investigate its stability, it is useful to present an alternative solution of the LDDE (A.1) using the exponential ansatz x(t) ∝ exp (−αt). Inserting this ansatz into Eq. (A.1) leads to the equation α = ω exp (ατ ) for the matrix α. Except for some notable exceptions [78], the solution of this equation is given by where W denotes the matrix valued Lambert W function [79]. The Lambert W function is a multivalued complex function. The long-time behavior of solutions to Eq. (A.1) and thus also the stability of its fixed point are determined by the branch of W yielding the largest real parts of the eigenvalues of the matrix α. The corresponding values of the Lambert W function strongly depend on the dimensionless delay τ ω.
For example, in one dimension, where ω is a positive real number, the branch of the Lambert W function with the largest real part exhibits three qualitatively different regimes as a function of τ ω leading to three different dynamical regimes of solutions 1/t R = (α), Ω = | (α)|, to Eq. (A.1). The boundaries between these regimes can be determined analytically [25,27,30]: (i) 0 < τ ω ≤ 1/e, where α is real and positive and x(t) decays exponentially to 0 with lifetime t R ; (ii) 1/e < τ ω < π/2, where α is complex with a positive real part, producing exponentially damped oscillations of x(t) with frequency Ω and lifetime t R ; and (iii) π/2 < τ ω, where α is complex with a negative real part, corresponding to exponentially diverging oscillations of x(t) with frequency Ω. For τ ω = π/2, 1/t R = 0 and λ(t) oscillates with frequency Ω without any decay.
In the panel a) of Fig. A1, we show that for long times the Green's function λ(t) for Eq. (A.1) is well approximated by the exponential solution (A.7). The above described dynamical regimes are reflected in the behavior of the decay rate 1/t R and the frequency Ω of oscillations, plotted in the panels c) and d), respectively. The panel b) of the figure shows the steady auto correlation of x(t) calculated in Appendix C, which is also well described by the formula (A.7).

Appendix B. Solution with Noise
Appendix B.1. One Dimension Let us now solve Eq. (A.1) with an additional noise term. For simplicity, we present the detailed derivation first in one dimension. We thus want to solve the equatioṅ where η(t) is a white noise fulfilling η(t) = 0 and η(t)η(t ) = δ(t − t ). 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. (B.1) as where λ(t) is given by Eq. (A.5) and ψ(t) = −ω 0 −τ ds λ(t − τ − s)x(s) is determined by the initial condition x(t) for t < 0. As we show in the next section, the formula (B.2) can already be used for the calculation of the time-correlation function for x(t). Here, we differentiate the solution (B.2) which yields the time-local Langevin equatioṅ . The dashed lines delineate the overall exponential decay exp (−t/t R ) of λ(t) and C(t). The dotted line in panel b) depicts the stationary value C(0) of the variance given by Eq. (C.6). The lifetime t −1 R = (α) and the frequency Ω = | (α)| of the oscillations in λ(t) and C(t) are shown in panels c) and d), respectively. The dotted lines in these panels serve just as an eye-guide depicting the linear behavior τ /t R = Ωτ = ωτ . For the purposes of the appendix, we assume that the individual model parameters are dimensionless. Parameters used: ω = 0.9 and τ = σ 2 /2 = 1.
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 (B.3) is completely determined by the joint probability distribution P 2 (x, t; x , t |x 0 , 0) [35].

Appendix B.2. Higher Dimensions
Let us now consider the probleṁ with general matrices ω and σ and the vector η(t) of white noises fulfilling η(t) = 0 and η i (t)η j (t ) = δ ij δ(t − t ). Since this system of Langevin equations is linear, the one-and two-time probability distributions P 1 (x, t|x 0 , 0) and P 2 (x, t; x , t |x 0 , 0) for 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 based solely on the formal solution of the Langevin system (B.9) x(t) = λ(t)x 0 + t 0 ds λ(t − s)ση(s) (B.11) with the initial condition x(t) = 0 for all t < 0 and x(0) = x 0 and with the Green's function λ(t) of the Langevin system given by Eq. (A.5). The mean value µ(t) = x(t) and the elements K ij (t) of the covariance matrix K(t) = x(t)x (t) − x(t) x (t) defining the PDF (B.10) can be obtained by inserting x i (t) from Eq. (B.11) into the definitions and averaging over the noise η(t). These formulas can be generalized straightforwardly to arbitrary initial conditions, where x(t) for t ≤ 0 is drawn from some probability distribution P [x(t), t ≤ 0]. Then the formal solution of the system (B.1) reads x(t) = y(t) + The mean value µ(t) is given by µ(t) = y(t) x(t≤0) , and K(t) = y(t) [y(t)] x(t≤0) + t 0 ds λ(t − s)σσ λ (t − s). The averages • x(t≤0) above are taken with respect to the PDF P [x(t), t ≤ 0]. The long-time behavior of the covariance matrix (B.13) 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 can be evaluated numerically. It is possible to rewrite it in a simpler form. Taking the derivative of 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 λ(t), i.e.
The restriction t > 0 for validity of this equation comes from discontinuity (and thus non-differentiability) of λ(t) at t = 0. The solution to Eq. (C.2) is given by Eq. (A.4) and hence the stationary space-time correlation matrix can be written as where C 0 = C(0) lim t→∞ 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 C(t), t ∈ [−τ, 0]. Fortunately, this initial condition can be calculated using the approach of Frank et al. [31] who calculated the time-correlation function C(t) for t ∈ [0, τ ] in one-dimension (see also Ref. [32]). They utilized the symmetry C(t) = C(−t) following from Eq. (C.1) to rewrite the delay differential equation (C.2) asĊ(t) = −ωC(τ − t) for t ∈ (0, τ ). Taking the derivative of this equation and using Eq. (C.2) yields the second order ordinary differential equation for the initial condition C(t) = C(−t), t ∈ [−τ, 0]. The solution of this equation reads C(t) = C 0 cos(ωt) +Ċ 0 ω −1 sin(ωt), where we still need to determine the unknown coefficients C 0 = C(0) andĊ 0 = lim t→0+Ċ (t). To this end, we need to evaluate independently C(t) and/orĊ(t) for two times t ∈ (0, τ ). Specifically, we show at the end of this appendix that C(τ ) = 0.5ω −1 σσ andĊ 0 = −0.5σσ . Using these results, the final expression for the correlation matrix for t ∈ [−τ, τ ] reads C(t) = C 0 cos(ωt) − 0.5σσ ω −1 sin(ω|t|) (C.5) with the initial value C(0) = C 0 = lim s→∞ K(s) = 1 2 ω −1 σσ + σσ ω −1 sin(ωτ ) cos −1 (ωτ ) (C.6) given by the stationary value of the covariance matrix. The whole time-dependence of C(t) for t ≥ 0 is thus described by the formulas (C.3), (C.5), and (C.6). The correlation matrix for negative times then follows from the symmetry C(t) = C(−t). Finally, let us note than in the one-dimensional case, where ω and σ stand for real numbers and, more generally, if the matrices ω −1 and σσ commute, we can rewrite Eq. (C.6) as where I denotes the identity matrix. An example of the time-correlation function for a one dimensional system is depicted in Fig. A1 b).
The expression for C(τ ) can be obtained by multiplying Eq. (B.9) by x (s), averaging the result over the noise, using the assumed stationarity of the process implying the formulaĊ(0) = d [lim s→∞ x(s)x (s) ] /ds = 0, and applying the symmetry C(t) = C(−t). The result is C(τ ) = ω −1 σ η(s)x (s) = 0.5ω −1 σσ , where the last equality comes after inserting the formal solution (B.11) for x(t) into the average, using the covariance of the noise, and noticing that in the resulting integral we integrate over half of the emerging δ-function only. The expression forĊ 0 then comes simply from Eq. (C.2), which is invalid for t = 0, but can be used for t arbitrarily close to 0 from the right, and using the result for C(τ ).