Corrigendum: Brownian molecules formed by delayed harmonic interactions (2019 New J. Phys. 21 093014)

A time-delayed response of individual living organisms to information exchanged within ﬂ ocks 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 toemulate 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 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 interactionstrength. We evaluate the entropy and information ﬂ uxes maintaining these structures and, to quantitatively characterize their stability, develop an approximate time-dependent transition-state theory tocharacterize transitions between different isomers of the molecules. For completeness, we include a comprehensive discussionof the analytical solution procedure for systems of linear stochastic delay differential equations in ﬁ nite dimension, and new results for covariance and time-correlationmatrices.

leads to a decrease in the steady state variance n = á ñ x where we have expanded the result to first order in τ in the final expression. The same result is obtained by expanding the exact expression (2). More details about small-delay expansions in delay Langevin equations can be found in [2]. The expansion (4) up to the first order in the time delay gives intuition why the behavior of equation (1) for small delays is reminiscent of an overdamped harmonic oscillator. Similarly, the expansion of equation (1) up to the second order in time delay gives a hint on why solutions to equation (1) for intermediate delays exhibit damped oscillations. However, one should be aware that especially the higher order Taylor expansions may be problematic [2,3]. Equation which is the same as the result (4), obtained using the first order expansion, and which differs from the second order expansion w wt w t + + ( ) D 1 not guaranteed that the remainder in the Taylor expansion is small and, second, the expansions involve higher derivatives that might not exist. In the case of the linear delay Langevin equation (1), we find that the second derivative of x involves the first derivative of the white noise, which should indeed be treated with care.
Following the mapping to the noisy damped harmonic oscillator further, we can identify the corresponding free eigenfrequency These quantities control the qualitative behavior of average solutions á ñ x to equation (5). For c < 1, thus wt < + » ( ) 1 1 2 0.41, it exhibits overdamped behavior (pure exponential decay), and, for c > á ñ x 1, performs exponentially decaying oscillations . The original model(1) exhibits overdamped behavior for wt < » e 1 0.37and damped oscillations for wt p < < e 1 2. Both approximate models (3) and (5) are unstable for wt > 1 and the original model is unstable for wt p = » 2 1.57. Even though the Taylor expansions above are of limited use, the result that the behavior of overdamped delay systems is reminiscent of the behavior of corresponding underdamped harmonic oscillators seems to be robust. Delay, in general, induces oscillations into the system dynamics, and, for strong potentials, may lead to instabilities. As another example of such behavior, we refer to [4] showing that the position x described the stochastic delay equation x t x t D t 2 9 exhibits qualitatively the same behavior as the velocity =  v x of the constant force oscillator described by the formula In this case, the small-delay Taylor expansion makes no sense at all due to the non-analyticity of the absolute value.

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 Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI. previous frame. The setup allows to create arbitrary time-delayed interactions in the many-body system. In [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 [11], the two-dimensional N-particle system is described by a set of 2N coupled nonlinear stochastic delay differential equations (SDDE). 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 [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 SDDE 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.

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 SDDE [29], which exhibit non-Markovian dynamics. Due to delayinduced 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 nonlinear 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 [30,31,34]. 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 [28,[51][52][53]. 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.

Outline
The rest of the paper is structured as follows. In section 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 section 3. In section 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 section 5. In order to assess the robustness of our findings, section 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 section 7. Most of the technical details are given in appendices A-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 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 harmonic pair interactions given by the potential figure 1 by springs connecting the individual particles. In equation (1), R>0 denotes the equilibrium spring length, k their stiffness, and r t t t r r )| is the distance between the particles i and j located at positions r i and r j at an earlier time t − τ. Clearly, the picture of linear springs can properly represent the time-delayed interactions only for a vanishing time delay τ.
Altogether, the particles diffuse in the composity potential where the summation runs over all pairs (i, j), so that the ith particle is driven by the force ). Because at time t the particle feels the value of the potential corresponding to its position at time t−τ, here x i and y i denote the Cartesian coordinates of the position vector r i in the past. In effect, the N-particle system therefore obeys the set of nonlinear delayed Langevin equations 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 h comprise independent Gaussian white noises satisfying the relations At long times and short enough delays, the particles form molecular-like vibrating structures, while for long delays, they exhibit exponentially diverging oscillations. 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 earlier time t−τ.
The numbers α and β label the components of the vector i h , while i refers to the specific particle. Note that the noise and the friction in equation (3) are related by the fluctuation-dissipation theorem [59] for a vanishing delay τ = 0 only, and that the system is always out of thermodynamic equilibrium for τ>0 [56]. In order to obtain a model that would obey the fluctuation-dissipation theorem, one should consider a different noise correlation function (4).
For small enough values of the time delay, the particles form, after an initial transient period, highly symmetric molecular-like structures, some of which are displayed in figure 4(a) in section 3. For N=2 (dimer) and N=3 (trimer) the steady-state structures occupy the global minimum of the potential V. For N>3 the global minimum becomes inaccessible due the chosen two-dimensional geometry and the resulting structures are thus frustrated in the sense that some of the springs do not reach their equilibrium length in the steady-state. The structures are dynamical, due to the Brownian motion of the particles, which persistently kicks the system out of the minimum of the potential energy (2). The effect of the delay is that the system may exhibit exponentially decaying oscillations on its return to the energy minimum. The decay rate of these oscillations decreases with increasing delay, and, for delays larger than a certain threshold, their amplitude exponentially increases. This is because large delays induce in the system a 'swing effect', when the repulsive force from one side of the potential propels the particle to a 'higher' position at its opposite side, and so on.
Within the equilibrium model that obeys the fluctuation-dissipation theorem, the stationary probability density function (PDF) for positions of the individual particles would simply be the Boltzmann distribution ) with potential V, inverse temperature β, and normalization Z. However, the physical situation at hand, where the delay is interpreted as a result of a feedback control mechanism and thus is independent of the noise, requires the more involved description with equation (4) that leads to non-trivial non-equilibrium steady states. Consequently, the Boltzmann distribution can no longer be assumed. For the simplest case of a dimer with short delay time, we will find an approximately Gaussian distribution, corresponding to a ('deformed') Boltzmann factor at an effective temperature. For larger particle numbers and longer delay times, the situation becomes more complicated.
A similar system with a quasi-constant force between the particles (constant upto a change of sign at distance R), i.e. obeying the set of Langevion equations , K, N, with sign(x) denoting the signum function, was discussed earlier in [11]. The main difference from our setting (3) is that, in equation (5), the absolute value of the force does not depend on the interparticle distances and the particle number N. The main benefit of assuming the harmonic potential in equation (3) is that it allows much more complete analytical treatment. To allow for an easy comparison of the two models, we illustrate our results using parameters inspired by [11]. In the following, we first review some analytical results for stochastic dynamics of dimers and trimers. On this basis, we will return to the discussion of the emerging structures in section 3.

Center of mass
Similarly as for the dynamics considered in [11], the center of mass coordinate N r r i n i c 1 º å = ( ) of the system obeys the Langevin equation denotes Gaussian white noise satisfying equation (4) (with the labels i, j replaced by c) and the diffusion coefficient D c =D/N. Regardless of the interactions, the center of mass performs ordinary Brownian motion and, assuming the center of mass is at time t=0 located at the point r 0 , the PDF for r c reads

Dimer
Let us now consider the simplest case of two interacting particles. For N=2, equation (3) yields the system of four coupled equations of motion: 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 motion describing Gaussian white noise satisfying equation (4) with the vector components i, j replaced by the label r.
Projecting equation (10) on the direction of the bond at time t (by multiplication with t t t e cos , sin j j = ( ) ( ( ) ( ))) and on the direction perpendicular to the bond (by multiplying it with t t e sin , cos j j =j ( ( ) ( ))), we obtain the equations where j(t, t−τ)=j(t)−j(t−τ) denotes the change of orientation of the vector e(t−τ) during time τ. Above, we used the formulas r ≡ re, r r r e e j º + j˙˙a nd e e r r r r h h h º + j j . From symmetry considerations, it follows that the stationary PDF for the orientation must be constant in j.
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 k t t 2 cos , (12) is strong enough such that the terms proportional to [r(t−τ)−R]/R can safely be neglected independently of the noise strength, the formula (13) for the angle assumes the form 5 The corresponding transition PDF (Green's function) for change of the orientation by ].
Using this function, we average equation (12) over j(t, t−τ) obtaining )is the natural relaxation rate. Note that the same formula with ω τ substituted by ω is obtained by simply assuming that the change of the orientation j(t, t−τ) of the bond per delay time τ is small, i.e. for D R 2 1 2 t  . 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 equation (15) in more detail around equation (25) below.
Equation (15) is a linear SDDE which can be solved analytically for r , Î -¥ ¥ ( ) . In our setting, r 0  and thus we should solve equation (15) with a reflecting boundary at the origin. However, since we have assumed that r t R r t 1 t -- | ( ) | ( ) , we already work in the regime where r only seldom significantly deviates from R and thus the solution of equation (15) on the full real axis should approximate well the desired solution on the positive half-line. The solution of equation (15) for r , Î -¥ ¥ ( )and t 0  in terms of deviations of the bond length from its equilibrium length (which we call as shifted bond length), x t r t R, can be derived by several methods. We review two of them (time-convolutionless transform and Gaussian ansatz) for a general multidimensional linear SDDE in appendices A-B. Here, we present just the main formulas. Assuming that the system was initially in state x(0)=x 0 and that x(t)=0 for t<0, the formal solution of equation (15) for r , where the dimensionless Green's function . The symbol θ(x) in equation (18) stands for the Heaviside step function. For an arbitrary initial condition x(t), t<0 and x(0) = x 0 , the expression x 0 λ(t) in equation (17)  delay ω τ τ which is a dimensionless measure for the relevance of the delay relative to the natural relaxation time, the Green's function λ(t) in equation (17) exhibits three different dynamical regimes discussed in detail in appendix A, in figure A1 and also below: (i) monotonic exponential decay to zero for short delays (ii) oscillatory exponential decay to zero for intermediate ('resonant') delays e 1 2   w t p t , and (iii) oscillatory exponential divergence for long delays ω τ τ>π/2.
The stochastic process x(t) in equation (17) arises as a linear combination of white noises and thus the corresponding PDFs must be Gaussian. Indeed, we find that one-and two-time conditional PDFs for x(t) with the initial condition δ(x) for t<0 and δ(x−x 0 ) at t=0 read , denote the mean of the shifted bond-length (16), its variance, and normalized time correlation, respectively. The function P x t x x , , 0 d 1 0 ( | ) 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, ) 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-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 equations (B.6) and (B.7) in appendix B, respectively. Interestingly enough, the diffusion and drift terms in the FPEs are given by the natural scales 2D and ω τ only in the limit τ→0. Moreover, both coefficients acquire a time dependence, determined by the function λ(t).

Specifically, the diffusion and drift coefficients in equation
. 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 τ=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 x x x D , ; , 0 exp 2 ). For a nonzero delay in the regimes (i) and (ii), i.e. when the noiseless solution (18) 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 , see figure 8 in section 5.1. In these cases, our approximate model thus predicts that, in the long run, the delay merely 'deforms' the (approximate) Gaussian equilibrium distribution through a parameter renormalization. This state is characterized by a nonzero entropy production rate [56]. For long time delays, no stationary state exists. In comparison to the analogous setting with a piece-wise constant force discussed previously [11], this destabilization for long delay times τ is a new feature, due to increasingly high systematic forces that may occur for long delays.
In the regimes (i) and (ii), the variance t x t x t 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 x s x s t The formula (24) finally allows us to specify the regime R ss 2 n  where the approximation (12) and (13) is reasonable because the PDF for r is relatively sharply peaked around the mean bond length R. As we already noted, equation (15) is also valid in the case D R 2 1 2 t  when the bond rotates only slightly in each delay interval and thus the angle j(t, t−τ) in equations (12) and (13) is small. However, also in this case we need to additionally assume that R ss 2 n  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 An example of the stochastic evolution of the dimer obtained from BD simulations of the exact system (12) and (13) is depicted in figure 2(a). In figure 2(b), we compare the results obtained from BD simulations with the time evolution of the average shifted bond length (21) for different values of the equilibrium length R. As expected, the approximate analytical formula (21) describes well the exact result for large enough R satisfying the inequality (25). For larger values of ν ss /R 2 , the analytical result underestimates the correct value. This is because the bond length in the BD simulation is obtained from equation (12) with the reflecting boundary at the origin, while we allow negative values of r(t) in the approximate analytical description. Similarly as for the mean value, the analytical formula (22) for the bond length variance ν(t) approximates very well the value obtained from BD simulations for large enough R, as shown in figure 2(c). In figure 2(d), we depict the monotonous rapid If not specified otherwise in the description of the individual panels, we used the experimentally reasonable parameters: ω=2k/γ=10 1 s −1 , τ 1 =0.1 s, D=1 μm 2 s −1 , and R=10 μm. In the BD simulation, we averaged over 10 4 trajectories with time step dt=10 −3 s. divergence of the stationary variance (24) as the time delay τ reaches the unstable regime (iii). This means that the delay tends to delocalize structures. On the other hand, the variance can be optimized as a function of the frequency ω τ . The best localized structure is obtained for the frequency fulfilling the formula ) . For the corresponding value 0.74 of the reduced 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, equation (3) gives the system of six coupled equations of motion: 2 .
and similarly for t r 13 ( ) and t r 32 ( ). To get analytical results for bond lengths r t t r ij ij = ( ) | ( )|, we multiply the formulas for t r ij ( ) by the corresponding unit vectors t t t e r r ij ij ij = ( ) ( ) | ( )| and linearize the resulting equations. To this end, we need to deal with the expressions e α (t−τ)·e β (t), α, β=I, K, III, where we introduced roman numbers as a shorthand indexing I≡12, II≡13, III≡32. For a vanishing delay τ=0, e α ·e α =1 and the scalar products e α ·e β describe the angles of the triangle formed by the three particles (see figure 1 in section 2). We find that up to the leading order in the equilibrium bond length R the triangle is equilateral and thus the internal angles are π/3, leading to the relations e I ·e II =e I ·e III =−e II ·e III ≈1/2+ O[(r α −r β )/R] with a correction that is on the order of (r I −r II )/R for e I ·e II and similarly for the other scalar products. The linearized equation for the relative coordinate x α ≡r α −R thus reads where the lower index mod III a a º is considered as periodic with the period 3, i.e. x IV ≡x I and x V ≡x II . Similarly as in the case of the dimer, equation (27) describes the dynamics of x α (t) well for large equilibrium bond lengths R and for time delays small compared to reorientation times of the unit vectors e α .
For an analytical treatment, it is advantageous to rewrite the system (27) in the matrix form , , where e j a denote jth component of the two-dimensional unit vector e α . The time-correlations between the three components of the noise vector t x ( ) are not mutually independent and read Due to the linearity of equation (28) and Gaussianity of the noise, the Green's function for the one-time PDF P t r r , , 0 1 0 ( | )is Gaussian [37], and determined by the mean value t t x m = á ñ ( ) ( ) and the covariance matrix We review in detail the derivation of these functions in appendices A and B.2.
For the initial condition x(t)=0, t<0 and x(0)=x 0 we get where λ(t) denotes the Green's function for equation (28) given by ] and comparing the result to the one-dimensional Green's function (18), we find that λ(t) undergoes with increasing t (i) monotonous exponential decay to 0 for e 0 3 2 1  wt < , (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  denotes the identity matrix. This formula follows from the results of appendix C after substituting the matrices ω and  ss 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 s s t x x lim s The regime of stability 3ωτ/2<π/2 of the trimer can also be determined from the condition that the matrix M cos wt ( )is not singular, i.e. its determinant cos 3 2 cos 3 2 wt wt ( ) ( )is nonzero. Due to the symmetry of the problem, all diagonal elements of the matrix ss  are identical and the same holds also for all its off-diagonal elements. The approximate analytical, time-dependent solution (32) for the covariance matrix is compared to the exact covariance matrix obtained by BD simulations of the complete model in figure 3. Given the approximations made, we find very good agreement. The analytical results only slightly underestimate the diagonal elements (probably for the same reason as for the dimer) and overestimate the off-diagonal elements. The behavior of the covariance matrix as a function of the frequency ω and delay τ is similar to the behavior of the variance (24) for the dimer. Specifically, the diagonal elements of ss  monotonicly increase (the PDF for the bond lengths become broader) with τ and exhibit a minimum as functions of ω, opening a possibility to optimize the width of the bond length PDF. The off-diagonal elements of ss  monotonously increase (the individual bonds of the trimer become more correlated) both with the delay and with the natural relaxation 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 m ¥ = ( ) , respectively, and (iii) an exponential divergence. The performed BD simulations confirmed that for large R, when the model is well described by the approximate analytical formulas, these regimes can indeed be observed also in the complete model (3). Furthermore, the analytical study predicted that the dimer is in the unstable regime (iii) for ωτ>π/2 and the trimer for ωτ>π/3. Let us now discuss how general the presented findings are.
The stationary average bond length μ can be determined by minimizing the potential energy V V r i j ij with the two-particle potential V 2 given by equation (1). Minimizing the potential in our twodimensional geometry yields the highly symmetric molecule-like structures shown in figure 4(a). Due to the confinement to 2d, the global minimum of the potential corresponding to μ=R is accessible only for the dimer (N=2) and the trimer (N=3). For larger molecules, the average bond length decreases as a result of the infinite range of the potential. The system asymptotically relaxes to the depicted structures if the noise D vanishes and the reduced delay time ωτ is small enough such that the system is in the dynamical regime (i) or (ii). Nonzero D leads to fluctuations around the asymptotic structures and large ωτ causes exponentially diverging oscillations.
We have solved the complete model using BD simulations and depict the behavior of the average bond length r(t) for several values of N in the dynamical regimes (i) and (ii)-(iii) in figures 4(b) and (c), respectively. In the regime (i), we observe that larger systems relax faster than those with smaller N. Furthermore, in figure 4(c), we see that larger systems oscillate with larger amplitudes and that the threshold between the regimes (ii) and (iii) is reached at smaller values of ωτ. More precisely, all the curves in figure 4(c) are plotted using the same parameters except for N and, while the curves for N3 are in the regime (ii), the curves for N>3 correspond to the regime (iii). These observations are in accord with our analytical findings for the dimer and trimer.
By analyzing the mean bond length at late times, we have evaluated the critical reduced delay  wt ( ) determining the threshold between the regimes (ii) and (iii) for N=2, K, 10. In figure 5, we show 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 ]. However, the analytical results would imply that the constant C equals to 2, which is smaller than the value C≈3 obtained from the complete model. The actual dimer and trimer are thus more stable than their linearized versions considered in our analytical study. The found scaling Other parameters used: α = 1 s −1 , R=10 μm, and D=1 μm 2 s −1 . In the simulation, we averaged over 10 4 stochastic trajectories with time step dt=10 −3 s. c crit ≈C/N implies that the stability of a system with rescaled potential stiffness k→k/N (to render the energy extensive [81]) would be (almost) independent of the particle number.
To better understand this behavior, let us first consider the approximate analytical model described in section 2.2. Imagine a particle in a harmonic potential centered around x=0, which is initially located at x(0)=x 0 >0. Assuming that x(t)=0 for t<0 and neglecting the noise, the particle does not feel any force in the time interval tä[0, τ] such that x(t)=x 0 for all tä[0, τ]. In the subsequent time interval tä[τ, 2τ], the particle experiences the force F=−ω τ 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 x 1 0 > | | | |. A similar process then repeats when the particle returns back, with the difference that now it reaches a maximum position 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 doublewell potential depicted figure 7 in section 5.1. We assume that the particle starts in the right well and oscillates with increasing amplitude as discussed above. After some time the amplitude becomes large enough that the particle crosses the barrier to the left well. Due to the presence of the additional well, the potential now contains much wider low-energy region compared to the purely harmonic case. The particle needs longer time to travel from one (unbounded) side of the potential to its other side, and hence also the (reduced) delay ωτ required for inducing diverging oscillations is larger than in the harmonic case. As a consequence, the complete model is seen to be more stable than foreseen by our analytical considerations. Moreover, our discussion reveals the existence of a fourth dynamical regime, preceding the unstable regime (iii), where the particle hops between the individual wells of the potential and the amplitude of the oscillations remains finite.
Compared to the quasi-constant velocity model investigated in [11], our analysis thus reveals two qualitative differences. First, the structures formed in the quasi-constant velocity model oscillate for arbitrary nonzero delay τ, while in the harmonic model these oscillations appear only if ωτ is large enough. Second, the amplitude of the oscillations in the quasi-constant velocity model is always constant in time, while the osculations in the harmonic model either vanish with time, if the system is in the regime (ii), or explode with time, if the system is in the regime (iii). The behavior observed in the harmonic model can be traced back to the increase of the force with the particle distance and thus we expect an analogous behavior also for other systems with time-delayed forces increasing with distance.

Entropy fluxes
The investigated model, much as the model discussed in [11], is inspired by self-organized systems, where a feedback based on the information about the state of the system at a previous time leads to structure formation. Interpreting the delayed interactions in our model as a result of such feedback control, we can investigate the entropy flow out of the system caused by the feedback. Due to the non-analyticity of the model with quasiconstant forces considered in [11], the analysis of entropy flows in the supplementary material therein was performed for vanishing delay only. Using the approximate Gaussian PDFs found in sections 2.2 and 2.3 for the dimer and the trimer, respectively, we can repeat this analysis with nonzero delay.
Without feedback, i.e. without the time-delayed harmonic interactions (2), the particles would spread diffusively and the system entropy would increase accordingly. The feedback control utilizes the information about the particle positions to drive the system into a non-equilibrium steady state with a time independent PDF P(x) and thus with a time constant configurational entropy 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 −1 ) performed by the feedback device against thermal dispersal. Evaluating the stationary entropy production F  due to the feedback control mechanism and the stationary entropy production D  due to the breaking of the fluctuation-dissipation theorem in equations (3) and (4) for τ>0, one can define the feedback efficiency as the ratio F The entropy production D  can be calculated along the lines of [56]. The entropy production q T F H  = 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 [11,64]. The center of mass coordinate of the system is not affected by the feedback and diffuses freely (see section 2.1). The structure formation due to the feedback thus occurs only on the level of the bonds. Let us now consider the time-dependent PDF P(x, t) for the bonds that converges to a time-independent non-equilibrium steady state due to the competition between feedback and diffusion. The rate of change of its Shannon entropy t 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 In this equation, the fist term on the right stands for the diffusive spreading of the PDF. The term t x, [ ] corresponds to the time-delayed force F[x(t−τ)] in the Langevin equation and thus it describes the effect of the feedback. Its concrete form is not relevant for the discussion below and thus we refer to the works [30,65] for more details about its structure.
Inserting equation 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 . It is interesting to adopt the Seifert's idea of trajectory-dependent entropy [57,66] and use equation (39) to define the stochastic (position dependent) entropy flux The average flux (39) then follows as the average t s t x,  = á ñ --( )˙( ) either over the PDF P(x, t) or over the individual stochastic trajectories generated in a BD simulation. To the best of our knowledge, the statistics of the entropy flux (40) has not been investigated yet and thus it is not known whether its PDF fulfills some fluctuation symmetries. Such an investigation would clearly be beyond the scope of the present paper and we leave it for a future work.
Let us now evaluate the three entropy fluxes (36), (38) and (39) for a general d-dimensional Gaussian PDF leading to the rate of change From equation (38) we then find that where λ(t) is given by equation (33)  The formulas (45)- (48) are valid both in the stable regimes (i) and (ii), where the system at long times relaxes to a stationary time-independent structured state, and in the unstable regime (iii). In the unstable regime, the variance ν(t) and the covariance t ( ) diverge in time. As a result, the system entropy t ( ) diverges and the entropy flow 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.
)˙( ). Let us therefore 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 in these two regimes. Using the asymptotic formulas (24) and (34) for the dimer bond-length stationary variance ν ss and trimer covariance ss  , we find from equations (45) for the dimer and for the trimer. In the last two formulas, a 1, 1 ss  = ( )denotes diagonal and b 1, 2 ss  = ( ) off-diagonal elements of the covariance matrix. The system entropies (49) and (51) are determined by the width of the PDFs for the bonds. Therefore, they monotonously increase with temperature T=γ D/k B and with the delay time τ, and exhibit a minimum as functions of the frequencies ω τ (dimer) and ω (trimer), similarly as the variance ν ss and the diagonal matrix elements of the covariance matrix ss  . The quality of the steady-state structures is thus in our model always unfavorably influenced by the delay and, for a given delay, one can tune the frequency in order to minimize this (usually unwanted) effect.
The two entropy fluxes (50) and (52) are negative, highlighting that they correspond to entropy outflows from (or information inflows into) the system. Interestingly enough, the entropy fluxes do not depend on the temperature T (or noise strength D) as already predicted in [11]. This means that the fluxes are discontinuous in the formal limit T → 0 because they must inevitably vanish for zero noise, where the PDF for the system evolution is a δ-function for all times. We plot the stationary entropy fluxes (50) and (52) as functions of the delay τ in figure 6(a) and frequency ω in figure 6(b). Naturally, maintaining a stationary structure in a bigger system (trimer, dashed lines) requires a larger (more negative) entropy flux (or more information) than in the smaller one (dimer, solid lines). The maximum of the fluxes  -|˙|, depicted in figure 6(b), arises as a result of a competition between stronger confinement, corresponding to larger frequencies ω, and gradual destabilization with increasing ωτ, when the system enters the unstable regime (iii). For the figure, we used for simplicity the approximation ω τ ≈ω for the dimer.

Transition rates for isomer transformations
The particles within the molecular structures depicted in figure 4(a) may exchange their positions. Assuming the particles to be distinguishable, different arrangements of the same structure may arise which can be interpreted as different isomers of the same molecule. Their study can provide further insight into the stability properties of our non-equilibrium molecules. In fact, we find that the study for molecules made out of only a few particles is informative also for the phenomenology observed for large particle numbers. While for the purely deterministic motion then isomer transitions only appear for time delays τ in the unstable regime (iii), in a system affected by thermal fluctuations the transitions occur for arbitrary delays. The evaluation of the frequencies of such transitions, which measure the stability of the individual isomers, can thus provide insight into the overall energy landscape responsible for the non-equilibrium structure formation. It is the main topic of transition rate theory [23].
The transition rate t A B k  ( ) 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 . 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.  (panel (b)), respectively. If not specified otherwise we used the parameters D=1 μm 2 s −1 , τ=1 s and ω= 1 s −1 . For the dimer, we used the approximation ω τ =ω, which is accurate for Dτ/R 2 ≈0.
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 t A B  ( ) 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 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 t lim t A B k ¥  ( ) of the transition rate can (approximately) be calculated using Kramers' rate theory [22,23] which was originally developed to describe chemical reaction rates. The approximation works best for a high energy barrier compared to the thermal energy. Kramers' theory was extended to reaction rates for GLEs describing non-Markovian systems. A crucial contribution in this direction came from Grote and Hynes [67] who derived a dynamical correction to Kramers' result. While their analysis was based on a parabolic barrier, Pollak [68] investigated the decay rate of an underdamped particle trapped in a symmetric cusp double well potential obeying the GLE with an arbitrary memory kernel satisfying the fluctuation-dissipation theorem. The timedependent rate t A B k  ( ) for driven overdamped systems can be calculated using the recent theory of Bullerjahn et al [24] for forcible molecular bond breaking.
To the best of our knowledge, the literature on the rate theory of time-delayed systems is scarce. The escape from a cubic metastable well under a time-delayed friction was investigated in [69]. Based on their small-delay approximation, Guillouzic et al [70] calculated the transition rate and the mean first passage time for an overdamped Brownian particle in a delayed quartic potential. From an experimental point of view, Curtin et al [71] studied transitions in a bistable system under time-delayed feedback.
Our model does not belong to any of the previously investigated classes of systems. However, for a vanishing delay one can use Kramers' theory, since the system obeys a Markovian overdamped Langevin equation. Moreover, for nonvanishing delays in the stable regimes (i) and (ii), the one-time PDFs for dimer and trimer can be described by standard (time-local) FPEs with time-dependent coefficients, where Bullerjahn's theory applies and where one can evaluate the transition rate numerically. Furthermore, after long times, the coefficients in these FPEs become time independent suggesting that Kramers' theory may be applied also to obtain the longtime 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.

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 r r 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 r then diffuses in the cusped double-well potential V r r R 2 2 gw = -(˜) (|˜| ) , depicted in figure 7, with the diffusion coefficient 2D.
For a nonzero delay, based on the approximate solution (17) to the Langevin equation (15) assuming r , Î -¥ ¥ ( )and x(t) = 0 for t < 0, we have found that the one-time PDF P P x t P x t x , , ,0

This equation describes diffusion in the harmonic potential
with the time-dependent stiffness γω τ (t) (given by (B.4)) and the time-dependent diffusion coefficient 2D τ (t) (given by (B.8)).
The validity of equation (53) for P 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 κ(t) for switching between the isomers. In the Markovian case, the transition rate for surpassing the (effective) energy barrier γω τ (t)R 2 /2 at r=0 to the other isomer can be calculated from equation (53) with an absorbing boundary at x=−R [63]. We now review several methods suitable for this task, and compare the results to BD simulations of the complete model with energy barrier γωR 2 /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 r R 0 = -> .

Numerical method
We first consider the situation when the system dwells in the state x(t)=0 for t 0  and then starts to diffuse in the time-dependent potential (54). Then, the time-dependent Markovian rate κ M (t) can be determined from the equation t a a M a for the normalized PDF P x t P x t S t , , a a a = ( ) ( ) ( ) for the position of the particle surviving in the right well of the cusped potential [72]. Here, P a (x, t) is the solution to the FPE (53) with absorbing boundary at x=−R and S t xP x t d , We solved it numerically using the method presented in [73].

Bullerjahn's method
Alternatively, one can determine the rate approximately using the analytical theory developed by Bullerjahn et al in [24]. Therein, the rate is constructed from the (Gaussian) solution P 1 (20) of the FPE (53) with natural boundary conditions. Specifically, one approximates the probability current across the absorbing boundary by 7 In the last expression, the symbol Erf denotes the error function and the variance ν(t) is given by equation (22). The approximate Markovian transition rate is then given by In figure 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 w l l = t ( )˙( ) ( ) and D t D t t t 2 2 l w n = + 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, , is controlled by the relaxation time t R for decay of λ(t) to 0, see also figure A1 in appendix A. The effective diffusion coefficient converges to the value D D D 2 1 sin , determined by the stationary variance ss n n ¥ = ( ) , see equation (24). 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 figures 9(a) and (b) below. For long delays, the (effective) barrier height over the (effective) thermal energy decreases so that the assumptions of the transition state theory are not valid any more, and κ B <κ M , see figure 9(c).
The situation of low barriers can be understood from the behavior at vanishing potential strength ω→0, when D τ (t)=2D and the finite time transition rates κ M (t) and κ B (t) can be calculated analytically. Namely, the PDFs P ã and P 1 in the definitions (56) and (60)   Since the denominator in the expression for the rate κ M (t) is smaller than that for κ B (t), we conclude that, for low energy barriers, the inequality t t On the other hand, for very high barriers, the PDFs P ã 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)≈1, and j j D P 2 0 To gain a deeper insight into the behavior of the transition rates, let us consider the stationary regime, t  ¥. In this regime, P R, 0 a -¥ = ( ) , due to the absorbing boundary condition at x=−R, and j x, , due to the conservation of probability P x, The smaller the frequency ω, the wider the PDFs P 1 and P ã . 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 of P ã . In such a case, the PDF P ã , which must vanish at x=−R, changes near the boundary faster than P 1 , leading to j j in accord with the argument put forward in the previous paragraph. With increasing ω, the boundary shifts away from the maxima of P 1 and P ã towards their tails. Due to the trajectories trapped in the absorbing state [72,74], the maximum of P ã is slightly farther away from the absorbing boundary than the maximum of P 1 , and thus, with increasing ω, the tail of P ã , with small derivative (small j Da ), hits the boundary at x=−R before the corresponding tail of P 1 . Hence, for large enough barrier height, the inequality between the rates crosses over to B  suggests that, also in this regime, the transition rate should saturate at long times. However, both the frequency ω τ (t) and the diffusion coefficient 2D τ (t) actually exhibit diverging oscillations caused by the oscillations in the Green's function λ(t) (18), in regime (ii). These divergences cause problems both in the FPE and in the approximate calculation of the rate using Bullerjahn's method. As a consequence, the (effective) Markov description can not be valid in the dynamical regime (ii). Nevertheless, let us now investigate to what extent the long-time transition rate k k = ¥ ( )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 P is time-independent and the limits t lim t w t ¥ ( ) and D t lim t t ¥ ( ) exist, we can rewrite the formula (55) as the eigenvalue problem  [22,23] for the transition rate for leaving one of the wells of a cusped potential with barrier height E b , In the penultimate equality, we used the appropriate inverse thermal energy , and also the asymptotic form of the diffusion coefficient D , 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 ) of the FPE (53). In the dynamical regime (ii), the kinetic prefactor in equation (65) cannot be correct, due to the diverging oscillations in the time-dependent frequency ω τ (t). Nevertheless, the exponential term seems to be reasonable, and thus it is tempting to use in the prefactor of the transition rates simply ω, instead of the problematic w ¥ t ( ). This substitution gives the correct pre-exponential factor of the rate for vanishing delay τ=0, where w w . 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 equation (63) is correct for Markov dynamics (τ=0), it can not be precisely valid for the time-delayed dynamics (τ>0). To see this, it is enough to realize that the delayed system arriving at the boundary at a time t does not feel the energy barrier E b =γ ω R 2 /2, but the barrier with energy ), 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 k obtained numerically from equation (63) with ω substituted for w ¥ t ( )in the operator . From the analytical expressions the re-scaled Bullerjahn expression B k (see equations (64) and (66)) works best. However, in the figure we have used parameters leading to a high barrier E b , and thus Kramer's and Bullerjahn's predictions, κ K and κ B , almost coincide. As a consequence, the line for K k in the figure overlaps with B k , similarly as the lineκ K (suppressed in the figure for better readability) overlaps with κ B .

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

Trimer
Consider the trimer depicted in figure 10 with the three distinguishable particles labeled by the numbers 1-3. We can count the particles either clockwise or anticlockwise and thus two different isomers can form in two dimensions. As in the case of the dimer discussed above, transitions between the two isomers occur with a transition rate depending on the diffusion coefficient, the coupling strength, and the equilibrium spring length.
There are several ways how the clockwise isomer may turn into the anticlockwise one and vice versa. For example, the particles 1 and 2 can switch their positions, or the particle 2 can migrate from above the line connecting particles 1 and 3 to below that line, to name a few. The transition rate for hopping between the two isomers is then given as a sum of the transition rates for all the possible realizations of the transition. In order to make an analytical prediction for the transition rate, we choose the coordinate frame in such a way that the x-axis always points from particle 1 to particle 3 (see figure 10). Then all possible transitions between the two isomers boil down to a single event when particle 2 crosses the x-axis. In particular, this also includes the transitions due to exchanging particles 1 and 3. In this case, the direction of the x-axes changes and thus the particle 2 effectively moves to its other side. In the following, we estimate the long-time transition rate k k = ¥ ( )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 j between the abscissas 12 | | and 13 | | and a neutral region of width x R 3 16 D = as exemplarily shown in figure 10. A neutral state 0ñ | is introduced to avoid over counting due to fluctuations of j around 0 and π. It is occupied if the smallest height of the triangle formed by the three particles, is smaller Δ  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 gray 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.
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 C R 12 b,eff 2 t gw = ( ) . Namely, we first consider 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 equation (2). The (possibly delay-dependent) unknown dimensionless prefactor C(τ) accounts for the time delay and all additional ways the particle 2 may take in the multidimensional energy profile in order to pass from 1ñ | to 1 -ñ | . The resulting Kramers' rate for the transition from 1ñ where a(τ) is a further unknown prefactor that may depend on all parameters of the model (see, for example, equation (65) for Kramers' rate in a cusped potential). In order to test the formula (68), we have fitted it to the transition rate κ obtained from BD simulations as a function of the diffusion coefficient D. The results for different values of the time delay are shown in figure 11(b) and the corresponding values of the coefficients C(τ) and a(τ), obtained from the fits, are given in table 1. The presented results prove that the transition rate exhibits exponential increase with the diffusion coefficient for a relatively broad range of values of the time delay, and thus the D-dependence of κ can be relatively well described by the Kramers-type ansatz (65).
Our attempts to include the effects of the delay in the effective diffusion coefficient and energy barrier, analogously to the approach described in section 5.1, where we made the replacements ,ss w w  t and D D ,ss  t , did not lead to a significant change of the value for C(τ). We thus conclude that the deviations of C(τ) from unity are mainly caused by the assumption that the particle will almost in all cases cross the axis at the minimum of the potential energy. Further improvement would thus require a multidimensional transition rate theory, which is clearly beyond the scope of the present paper. Properties of these systems are studied in [35]. A slightly unusual but interesting variation that makes sense in the context of active matter employs a negative delay. The individual active agents may react to a future state of their neighborhood which they predict in the basis of its present state. For idealized systems capable of a perfect prediction, the GLE (or, equivalently, the linear SDDE) contains the memory kernel f(t)=δ(t+τ), τ>0 and it can be solved using the strategy described in appendix A. The resulting Green's function

Extensions to other memory kernels
is the so called Bruwier series [75] that is convergent for e w t < | | | |. Alternatively, the series (71) can be written in the form [76] t s ). Therefore, we reformulate the deterministic part of the GLE (69) as x t x t . 7 3 One of the reasonable strategies for predicting x pre is to use the linear extrapolation x t x t x t , 7 4 enough, the entropy fluxes do not depend on the noise amplitude D and hence they are discontinuous at D=0, where the steady-state structures are stable without feedback and thus the entropy fluxes vanish.
Assuming the particles to be distinguishable, the steady-state structures (molecules) can form different isomers. Their transition dynamics can provide rich additional insight into the energy landscape underlying the non-equilirbium structure formation. For the dimer and trimer, we have investigated how and when the transitions between the individual isomers can be described by transition state theory. For the dimer, we have applied our analytical results, based on the time-convolutionless transform leading to the time-local FPE (53), to construct several analytical approximations for the transition rate using Kramers' theory [22,23] and Bullerjahn's theory [24]. We have also calculated the transition rate from the FPE numerically. Finally, we have compared the obtained predictions to results of BD simulations of the full problem. While the FPE gives the exact value of the transition rate for vanishing delay, our results show that the obtained rates agree with the true ones for small and moderate values of the delay only. We conjecture that this is caused by the fact that the classical absorbing boundary used in our numerical and analytical evaluation of the transition rate can not be used for larger values of τ. Concerning the analytical results, the best agreement with the true rates was obtained by the Bullerjahn's formula (64) with effective barrier height and diffusion coefficient taken from the time-local FPE (53) and the prefactor rescaled according to equation (66). In the case of the trimer, we have confirmed by BD simulations that the transition rate increases exponentially with the noise strength D even for longer delays and thus Kramers' or Bullerjahn's type predictions can be used also in this case. We plan to further investigate suitable absorbing boundary conditions for delayed systems to predict (at least numerically) transition rates also for large delays.
Finally in section 6, we have considered the robustness of our analytical results with respect to details of the realization of the delay. We demonstrated that most of the presented equations can be used also for systems with memory kernels different from that for discrete time delays, i.e. f(t)=δ(t−τ). It is enough to substitute the Green's function λ(t) (A.5) corresponding the delayed Langevin equation by the Green's function corresponding to the memory kernel of interest. We reviewed some paradigmatic memory kernels and provided an outlook on the differences and similarities of the corresponding Green's functions. A more detailed study is left for future work.
As a further extension of our work, it would be interesting to consider physically more realistic interactions that vanish at large distances. Furthermore, we plan to investigate the reaction of the studied system to an external perturbation. Of particular interest could be the propagation and decay behavior of a local perturbation through the system, especially in case of large numbers of particles. Last but not least, we aim to investigate the behavior of the studied system under the action of an additional deterministic time-dependent driving and study the corresponding stochastic dynamics and thermodynamics. . In the present paper, we always assume that the system is initialized at time t at position x 0 with a special history, namely x(t)=0 for all t<0, allowing us to simplify equation (A.4) to x(t)=λ(t)x 0 .
The only fixed point of the DDE (A.1) is x(t)=0. In order to investigate its stability, it is useful to present an alternative solution of the LDDE (A.1) using the exponential ansatz t t x exp a µ -( ) ( ). Inserting this ansatz into equation (A.1) leads to the equation exp a w at = ( )for the matrix α. Except for some notable exceptions [79], the solution of this equation is given by where W denotes the matrix valued Lambert W function [80]. The Lambert W function is a multivalued complex function. The long-time behavior of solutions to equation (A.1) and thus also the stability of its fixed point are determined by the branch of W yielding the largest real parts of the eigenvalues of the matrix α. The corresponding values of the Lambert W function strongly depend on the reduced delay τω.
For example, in one dimension, where ω is a positive real number, the branch of the Lambert W function with the largest real part exhibits three qualitatively different regimes as a function of τω leading to three different dynamical regimes of solutions I 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. , where α is real and positive and x(t) decays exponentially to 0 with lifetime t R ; (ii) 1/e<τω<π/2, where α is complex with a positive real part, producing exponentially damped oscillations of x(t) with frequency Ω and lifetime t R ; and (iii) π/2<τ ω, where α is complex with a negative real part, corresponding to exponentially diverging oscillations of x(t) with frequency Ω. For τ ω=π/2, 1/t R =0 and λ(t) oscillates with frequency Ω without any decay.
In the panel (a) of figure A1, we show that for long times the Green's function λ(t) for equation (A.1) is well approximated by the exponential solution (A.7). The above described dynamical regimes are reflected in the behavior of the decay rate 1/t R and the frequency Ω of oscillations, plotted in the panels (c) and (d), respectively. The panel (b) of the figure shows the steady auto correlation of x(t) calculated in appendix C, which is also well described by the formula (A.7). 1 1