Heat flow due to time-delayed feedback

Many stochastic systems in biology, physics and technology involve discrete time delays in the underlying equations of motion, stemming, e. g., from finite signal transmission times, or a time lag between signal detection and adaption of an apparatus. From a mathematical perspective, delayed systems represent a special class of non-Markovian processes with delta-peaked memory kernels. It is well established that delays can induce intriguing behaviour, such as spontaneous oscillations, or resonance phenomena resulting from the interplay between delay and noise. However, the thermodynamics of delayed stochastic systems is still widely unexplored. This is especially true for continuous systems governed by nonlinear forces, which are omnipresent in realistic situations. We here present an analytical approach for the net steady-state heat rate in classical overdamped systems subject to time-delayed feedback. We show that the feedback inevitably leads to a finite heat flow even for vanishingly small delay times, and detect the nontrivial interplay of noise and delay as the underlying reason. To illustrate this point, and to provide an understanding of the heat flow at small delay times below the velocity-relaxation timescale, we compare with the case of underdamped motion where the phenomenon of “entropy pumping” has already been established. Application to an exemplary (overdamped) bistable system reveals that the feedback induces heating as well as cooling regimes and leads to a maximum of the medium entropy production at coherence resonance conditions. These observations are, in principle, measurable in experiments involving colloidal suspensions.

A finite heat flow is a generic feature of systems out of thermal equilibrium. In the last decades, special interest has been devoted to heat exchange and other thermodynamic properties of small (mesoscopic) systems coupled to a bath, which are noisy per se 1,2 . Stochastic thermodynamics (ST) has emerged as an elegant and consistent framework to generalize thermodynamic notions to the level of noisy trajectories and to systems far from equilibrium 3 , with numerous applications to soft matter 4 , biological 5,6 , and quantum systems 7,8 . Many fundamental concepts, however, are based on the Markov assumption, although in real-world systems memory effects are in fact often not negligible. While extensions towards several non-Markovian systems have been carried out in the past [9][10][11][12][13][14][15][16][17][18] , the application to continuous systems with discrete time-delays (i. e., delta-like memory kernels in the equations of motion) is still in its infancy. Delays can be of intrinsic nature as in neural systems 19,20 and laser networks 21,22 , or be generated externally, e. g., by a feedback protocol with time lag between signal detection and action of control [23][24][25][26] . Moreover, since time-delay is known to enhance control strategies, it is often included intentionally, for example in Pyragas control 21 . Despite their major importance, delayed systems are still little understood from a thermodynamic perspective, especially in regard to the interplay of delay and nonlinearities.
Recent research [27][28][29][30] has revealed that ST of delayed continuous systems is indeed quite involved, even in the absence of nonlinearities. For linear cases, it has been explicitly shown that in the long-time limit a non-equilibrium steady state (NESS) with finite entropy production 30 is approached (in the absence of time-dependent forces), as the delay pushes the system out of equilibrium. Nonequilibrium inequalities have been found 28 that generalize the second law and provide bounds to the extractable work. Furthermore, the fluctuations of work, heat and entropy production have been investigated 29 . One crucial issue in the thermodynamic description of delayed systems is the acausality of time-reversed processes appearing in the path integral representation of fluctuating heat q and entropy production [27][28][29] . Contrary to the Markovian case 2 , the total average entropy production ΔS tot of a delayed system in a NESS differs from the medium entropy Δ = 〈 〉 S q m s s  (where  is the heat bath's temperature and 〈⋅ ⋅〉 ss denotes NESS ensemble averages). In particular, the second law does not impose nonnegativity on ΔS m alone [27][28][29] . However, while these statements are generic, explicit expressions for the thermodynamic quantities are, so far, only available for systems governed by linear forces [27][28][29][30] , thus excluding wide classes of physically interesting processes which exclusively arise in nonlinear systems.
As a step in this direction, we here apply ST to investigate the heat rate δ = 〈 〉  Q q t /d ss of a classical nonlinear delayed system. By considering =   Q S m  rather than  S tot , we avoid the above-mentioned problem induced by acausality, and at the same time consider a key thermodynamic quantity and nontrivial part of the total entropy production, which already provides important physical insight into the thermodynamics of delayed systems. This strategy enables us to address several fundamental questions: What is the impact of time-delayed feedback on heat exchange and entropy production? Has the overdamped limit consequences for the thermodynamic description? Do thermodynamic quantities reflect delay-induced dynamical behaviour? An example are spontaneous oscillations occurring in many nonlinear delayed systems due to their infinite dimensionality 31,32 .
To this end, we consider a simple exemplary system composed of an overdamped particle subject to nonlinear static forces and a deterministic (i. e., error-free), continuous linear feedback force with a discrete time delay τ ≥ 0. Such a force is imposed, e. g., by optical tweezers 33,34 . The idealized heat bath is assumed to remain at equilibrium, unaffected by the feedback. Our calculation explicitly predicts that delay alone induces a finite heat flow whose direction is tunable. This means, in particular, that the delayed force can generate a steady heat flow from the bath to the particle, i. e., feedback cooling. We moreover unravel an important universal heat flow caused by the nontrivial interplay between delay and noise, which also occurs for underdamped motion. We discuss its physical origin in detail and propose that this heat flow is linked to the entropy pumping known for underdamped systems with velocity-dependent feedback, i. e., "molecular refrigerators" 27,35,36 . We further detect discontinuous behaviour at τ → 0. Discussing the application to a paradigmatic bistable system 37 , we evaluate the heat rate via several approximations and by numerical simulations. In particular, we consider the medium entropy production near coherence resonance 32,[38][39][40], that is, the appearance of regular positional oscillations at a finite thermal energy caused by the interplay of nonlinearity, noise, and delay [41][42][43][44] . Importantly, in contrast to stochastic resonance, there is no periodic external driving. Combining the discretized (Master) equation approach of 41 with ST, we show that the medium entropy production has a maximum at CR, and provide an analytical explanation for this maximum, which has, so far, only been detected numerically 44 .

Model
We consider stochastic processes described by the overdamped Langevin equation (LE) 45 , where the total deterministic force F(x, x τ ) = F con (x) − bx τ depends on the instantaneous, X(t), and on the delayed particle position X(t − τ). We assume that the conservative force F con can be expressed as a polynomial (which holds indeed for a wide class of nonlinear potentials), while the feedback control F d is chosen to be linear corresponding, e. g., optical tweezers 33 . ξ denotes Gaussian white noise with ξ 〈 〉 = t ( ) 0 and 〈ξ(t)ξ(t′)〉 = δ(t − t′), while γ and D 0 are the friction and diffusion coefficients satisfying 46 γ = D k T 0 B , with k B being the Boltzmann constant. Due to the appearance of the delayed position in (1), the corresponding Fokker-Planck equation (given at the end of this report in the Technical aspects section) for the probability density function ρ(x, t) (PDF) is an infinite hierarchy 45,[47][48][49] . We consider natural boundary conditions, i. e., ρ = ∂ x ρ → 0 for x → ±∞, and focus on NESSs, where ∂ t ρ(x, t) = 0. Following the ST framework of Sekimoto 1,30 , the fluctuating heat δq flowing to the reservoir during the infinitesimal time dt, the increment of internal energy du, and the work δw done by the nonconservative (delayed) forces, are given by The ○-symbol indicates usage of Stratonovich calculus. Plugging the LE (1) into Eqs (3,4) results in the NESS ensemble averageṡ  Heat rate and medium entropy production. For nonlinear systems, the linear response function method used in 30 , which bases on the Laplace-transformation, cannot be applied. In the following, we present an alternative approach for the heat rate. In particular, we derive exact expressions for (5, 6) which only involve positional moments at one time. To this end, we utilize two relations that can be derived by projecting the Fokker-Planck equation onto the positional moments and subsequently inserting the LE 50 .
First, the C i (τ)-terms can be substituted by the relation 50 A derivation of Eqs (7,8) can be found at the end of this report in the Technical aspects section. Finally, we combine Eq. (8) with a causality argument to evaluate ) ( ) at τ ≠ 0. Because no physical quantity can be influenced by future noise, statistical independence follows and Substituting Eqs (7) and (8) at n = 1, 2, .., m, respectively, into Eq.
. This is expected, since the net impact of conservative forces should vanish in a NESS 2 .
. Substituting (7) at n = 1, 2, .., m into Eq. (6) (please note that this simplification step is possible due to the polynomial form of force F con ) and using (9) further yields  (10) is an exact expression only involving one-time ensemble averages over X n .  Q and  S m can therefore be computed directly from the steady-state one-time PDF, and hence, on the basis of several approximations known from the literature. The δ τ -term suggests discontinuous behaviour of  Q at τ → 0. But in order to study this limit properly, one also has to investigate the behaviour of the PDF to clarify whether the moments behave continuously. For nonlinear systems, this is a nontrivial task on its own, as, in fact, no exact solutions for the one-time PDF in the presence of delay are known. However, we can nevertheless address this question analytically, since the approximative PDFs become exact in the Markovian limits 45,47 , rendering exact results from Eq. (10).
A detailed description of the approximation schemes is beyond the scope of this report. We refer the interested reader to 45 and references therein, and here only review some main aspects in brief. The small τ approximation 47 is obtained by a first-order Taylor expansion around τ = 0 in the LE, effectively rendering a Markovian system (with exponentially decaying correlations). Two other schemes start from the Fokker-Planck equation for the one-time PDF, which involves the two-time PDF 47 ρ 2 (x, t; x′, t − τ) (see Technical aspects section). The force-linearisation closure 45 (FLC) approximates ρ 2 by the one from the corresponding linear delayed system. The perturbation theory 48,49 (PT) approximates it with ρ 2 of the corresponding Markovian process without delay force. For a doublewell potential, the PT requires an additional approximation, since the corresponding two-time PDF is not known. As in 45 , we use the small-time propagator (PT-st), or the Ornstein-Uhlenbeck (PT-OU) approximation. In the following, we will investigate the heat rate in the two (Markovian) limits, where the approximations become exact.

Limit of Vanishing Feedback Strength.
In the limits where Markovianity is recovered, one expects that the system equilibrates (because of the absence of a driving force), and therefore = =   Q S 0 m . The first limit corresponds to b → 0, i. e., vanishing feedback strength. By construction, the one-time steady-state PDFs from the FLC and from the PT become exact at b → 0. They converge in a continuous manner to the equilibrium (Boltzmann) distribution ρ = b ss 0 , as can be easily seen from the PT-st result 45 ρ ss PT : and normalization constant Z. After performing a partial integration step and plugging in the natural boundary conditions, the (i − 1) st positional moment at τ = 0 can therewith be expressed as    12), we obtain  Plugging both into Eq. (10) in an iterative manner, yields the discontinuous limit This is our first main result. This ubiquitous jump-discontinuity at τ → 0 (see Fig. 1 for an example) indicates an abrupt qualitative change of the thermodynamics when non-Markovianity sets in. It arises due to the discontinuity of the noise-position cross correlations (9) at the onset of causal relationship. Remarkably, the apparent offset of k B b/γ is independent of the details of the potential landscape. We note that the apparent offset has already been observed and discussed in the context of linear systems [27][28][29][30] . In 30 it has been considered as a consequence of inconsistent usage of Ito and Stratonovich calculus. However, as shown here, it is not a mathematical error but also arises within consistently applied Stratonovich calculus. The underlying reason is the interplay of white noise and delay below the short (ballistic) relaxation timescale 51 , which becomes relevant as τ → 0. To further elucidate the behaviour at τ → 0, we also briefly consider the limit in the underdamped case, where the LE con 0γ additionally involves an inertial term with mass m, yielding ballistic motion below the velocity-relaxation timescale, γ/m. For this system, the heat can, in principle, be calculated using the same expression as in the overdamped case, see Eqs (2, 4). As has been shown previously 27 for linear systems,  Q smoothly decays to zero with τ. However, www.nature.com/scientificreports www.nature.com/scientificreports/ the response function method used in 27 is only applicable to linear systems. For a better comparison, we here consider the correlation ) ( ) appearing in Eq. (6), which, in our framework, "causes" the discontinuity from a mathematical point of view [compare Eqs (9,17)]. As opposed to the overdamped case, this correlation is indeed continuous at τ = 0, (not only for linear systems), indicating that  Q behaves continuously. The derivation of (19) is given in the Technical aspects section.
Turning back to the overdamped system, apart from the discontinuous limit, we find from Eq. (17) that for τ > 0,  S m and  Q are nonzero, proving the true non-equilibrium nature of this steady state. According to (17), the rates have negative values at small τ, if b < 0, which is always the case in an optical tweezers setup (see below). The negative signs indicate a steady heat flow from the bath to the particle (feedback cooling). This is a delay-induced phenomenon, which would be impossible in the Markovian counterpart of this system due to the second law 36 . The fact that here <  S 0 m (the bath constantly loses entropy) underlines that further entropic terms must contribute to the (nonnegative) total entropy production (as discussed for underdamped dynamics in [27][28][29] ). In the cooling regime, the amount of medium entropy loss thus provides a lower bound to the entropic cost 27-29 of the feedback. Analytical evaluation of the heat rate for (linear and nonlinear) example systems, which we will provide below, reveals that the (negative) heat flow is actually quite pronounced in the regime of small τ (compare with Fig. 1). Notably, this is also true for underdamped linear systems 27 (where  Q only decays for τ below the velocity-relaxation timescale, compare with Fig. 1 in 27 ). Before proceeding with the concrete examples, we here propose an explanation of this phenomenon.

Discussion of the behaviour for small delay times. When τ is so small that
is a friction-like force. As is well-known for underdamped systems with velocity-dependent (delayed 27 or non-delayed 35,36 ) feedback, such a control amounts to medium entropy reduction due to entropy pumping. (Please note that a non-delayed velocity-dependent feedback ∼  F bX t ( ) d drives a system out of equilibrium 35,36 , contrary to a non-delayed position-dependent control). The reason is that the additional friction-like force reduces the thermal fluctuations of the particle, inducing a net energy transfer from the bath to the particle, hence, a net heat flow. (Reversely, the particle fluctuations are enhanced for b > 0, yielding a positive heat flow.) This means, an additional "entropy pumping" contribution to the heat flow arises, explaining the enlarged heat flow at small delay times. Please note that this effect is independent of the question whether the overdamped limit is used, or not, and should be measurable in experimental setups.
When τ gets below the velocity-relaxation time (m/γ), both system behave differently. For underdamped dynamics, the velocities are correlated, implying that the delay-induced heat flow ) ( ) m / eventually decays smoothly to zero (resulting in a maximum of | |  Q around τ = m/γ, see Fig. 1 in 27 ). For overdamped dynamics, the velocities are, in contrast, uncorrelated (because the white noise directly acts on it ξ ∼  X ), giving rise to the nontrivial limit τδ(τ) of type "0 × ∞". This yields a finite value as τ → 0, as we know from Eq. (17). Hence, in an experimental setup, a position-dependent feedback is expected to induce heat flow, unless the (typically unavoidable) delay is well-below the velocity-relaxation timescale.
Linear delayed system. For linear systems, where the force has the form F = −a 1 x − bx τ , the exact NESS PDFs are known 52 , and (10) simplifies to the exact and closed expression A plot of  Q as a function of τ is given in Fig. 1. As expected,  Q has an apparent offset at τ = 0 and | |  Q grows while τ decays to zero (until a saturation value is reached). Equation (20) is consistent with 30 (where a different approach is used that only applies to linear systems), apart from the additional δ τ -term. We stress that, consequentially, only Eq. (20) correctly predicts Application to a delayed bistable system. We now consider a bistable system composed of a doublewell  4 . This is a prototypical nonlinear noisy system which exhibits both: a nontrivial intrawell dynamics within the asymmetric potential wells, and a noise-and delay-induced dynamical state, i. e., positional oscillations between the wells [41][42][43][44] . For this system, Eq. (10) involves the even moments up to 6 th order. Performing Brownian dynamics (BD) simulations (see Technical aspects section for details) to compute the NESS heat rate, both, from Eqs (2) and (10), we indeed find perfect agreement. In the following, we compare the BD results to those obtained from our analytical approach, i. e., Eq. (10) combined with established approximations. Since the approximations are known to perform best when the particle is likely to stay around a potential minimum 45 , this approach seems (2019) 9:2491 | https://doi.org/10.1038/s41598-019-39320-0 www.nature.com/scientificreports www.nature.com/scientificreports/ appropriate in the low thermal energy regime (γ  D V 0 0 ). At the end, we will introduce a complementary approach for larger noise levels.
Low thermal energy -intrawell dynamics. Figure 1 shows =   Q W and  S m , as functions of the delay time τ, for an exemplary parameter setting in the low noise regime. The simulation results confirm the discontinuous τ-limit and apparent offset as given by Eq. (17), and the predicted feedback cooling. The PT-OU predictions are quantitatively very similar to the ones obtained by the FLC, while the PT-st yields very similar results to the small τ expansion. We will therefore mainly discuss the FLC and the small τ expansion in the following.
While the small τ expansion fails outside the Markovian limit, our analytical approach with FLC makes quantitatively correct predictions for all delay times considered in Fig. 1. The thermodynamic quantities grow with τ, until they approach constant values. Interestingly, the saturation occurs at about τ ≈ x D /(2 ) 0 2 0 , i. e., on the timescale where the mean-squared displacement 2D 0 t of a freely diffusing particle is in the range of the extent of a potential well ≈x 0 (i. e., the particle has explored the whole well within τ). x D /(2 ) 0 2 0 is thus an estimate of the intrawell relaxation time t rel in the delayed system. Figure 1 also displays the results from the corresponding linear system, F = −a 1 x − bx τ [with a 1 = (8V 0 + k)/x 0 and b = −k/x 0 , as given by a second order Taylor expansion of V s + V d around a (deterministically) stable fix point, i. e., (x, x τ ) = ±(x 0 , x 0 )]. Interestingly,  Q is not only equivalent to the nonlinear case for very small τ [as expected from (17)], but also saturates on the same timescale.
The agreement between BD and FLC persists at larger barrier heights V 0 or smaller k. However, for larger k/ (γD 0 ), also the FLC approach breaks down, as can be seen in Fig. 2. The inset shows that the FLC nevertheless captures the qualitative behaviour. Upon increase of k (i. e., the laser intensity in case of optical tweezers 33 ), one may switch from feedback cooling to heating in the nonlinear case.
High thermal energy -interwell dynamics. When the thermal energy is sufficiently high, such that jump processes between the potential wells dominate the dynamics, the discussed approximations all break down by construction. The reason is that it is then the "memory of a jump" expressed in the non-Markovian ρ 2 , which dominates the process and, e. g., triggers subsequent jumps, and precisely ρ 2 is what is (over)simplified in all approaches 45 . In particular, the FLC effectively assumes an harmonic potential well out of which no jumps can occur, while the PT and small τ expansion render Markovian systems which don't exhibit spontaneous oscillations, as the latter are a delay-induced phenomenon. This includes situations, where the interplay of noise and delay leads to spontaneous oscillations of X. However, we can still treat this regime via an alternative strategy, which bases on the discretised approach for multistable systems proposed by Tsimring and Pikovsky 41 . To this end, the X-dynamics is reduced to switching processes between two discrete states s = ±x 0 (corresponding to the two potential wells). There are two different transition rates corresponding to the two possible scenarios that the delayed and instantaneous state have the same, or opposite sign, respectively. These transition rates are approximated by the Kramers formula 41,45,53,54 within a quasistatic approximation. The latter is justified when the delay time is small compared to the intrawell relaxation time, τ  t rel , and γD 0 < V 0 . Further details including the calculation of the transition rates are given at the end of this report in the Technical aspects section.
The simplification of Eq. (6) due to the 2-state reduction is twofold. First, the spatial correlation function s s can now be calculated analytically 41 . Second, all other spatial autocorrelation functions are significantly simplified, since the intrawell-dynamics are omitted, yielding for τ > 0 ) ( ) 0 at τ > 0 due to causality]. Consequentially, Eq. (6) simplifies to  Fig. 1. www.nature.com/scientificreports www.nature.com/scientificreports/ where we have used that the mean change of internal energy must vanish in the NESS. Using the explicit formula for C 1 (τ) [see Eqs (23,34)], can readily be evaluated. The results are plotted in Fig. 3.
The numerical data in Fig. 3 reveals that, when the thermal energy is small compared to V 0 ,  Q increases linearly with γD 0 , as it also does for a linear delayed system [see Eq. (20)]. In this regime, approximating the doublewell by its linearised version (dashed line) even renders the correct slope. By contrast, in the 2-state model =  Q 0, as expected, since the intrawell dynamics is neglected. At larger γD 0 /V 0 , the slope of  Q abruptly changes and nonlinear behaviour sets in. This occurs at the onset of the delay-induced oscillations, as reflected by a sudden increase of the CR order parameter from the 2-state model 41 and as confirmed by numerical simulations (see Fig. 4). The order parameter measures the height of the main peak at a frequency of about 2π/τ in the power spectrum of C 1 (Δt). Hence, a nonzero order parameter indicates the occurrence of subsequent escape events with period τ.
By further increasing γD 0 /V 0 , one enters the regime where delay-induced oscillations dominate the dynamics. At a certain finite thermal energy (of about 0.36 γD 0 /V 0 in Figs 3 and 4, see cyan vertical lines and cyan plot in Fig. 4a.), the peak in the power spectrum is highest, resulting in a maximum of the order parameter and indicating that the delay-induced oscillations with mean period ≈τ are most pronounced. In this range of thermal energies, the escape times are comparable with the delay time τ, and the pronounced peak in the power spectrum signals the resonant response of the system to the noise. Please note that the escape rates can be estimated by the Kramers rates in the corresponding system without delay (τ = 0) (see Technical aspects section). We find that in this range of γD 0 /V 0 , the 2-state reduction renders a good approximation. In particular, it predicts accurately that a region of steep slope is followed by a lower slope of  Q accompanied by a maximum of  S m , and that the latter lies in the regime of CR. This is our second main result.
Analytical explanation for the maximum of entropy production at CR. The behaviour of  Q and  S m can be understood on the basis of the 2-state model. The delay-induced oscillations set in and pause randomly. They have mean period τ, such that C 1 (τ) = 1 for a perfect oscillation (as in the case of no jumps). However, due to their stochastic nature, C 1 (τ) is lowered with the occurrence of oscillatory events. When the timescales of noise-induced escapes and oscillatory motion become comparable, the particle dynamics responds resonantly to the noise. Slightly increasing γD 0 /V 0 then significantly increases the number of occurring oscillation periods, yielding a strong reduction of C 1 (τ).  Q thus steeply increases and  S m reaches high values [Eq. (23)]. At the CR maximum (vertical cyan lines in Figs 3,4), this effect saturates resulting in a reduced slope of  Q and a maximum of  S m at CR. For even higher γD 0 , the superimposing noise generates irregular, low correlated motion with C 1 (τ) → 0, hence →  S 0 m . In contrast,  S m of the full system approaches a (nonzero) constant (not shown here). The breakdown of the 2-state approximation in this limit is indeed expected since the state discretisation then becomes meaningless. We have tested several parameter settings confirming that the  S m is indeed maximal at CR conditions as predicted by the 2-state model. Quantitatively, we always observe an overestimation of the thermodynamic quantities by a factor of about three. This is somewhat surprising since the phase-space reduction inherent  www.nature.com/scientificreports www.nature.com/scientificreports/ to the discretisation is rather expected to yield an underestimation of the steady state quantities 3 . A possible explanation lies in the quasistatic approximation, but the precise reason is subject of future investigations.

Conclusions
In this work we put forward analytical approaches for the heat rate in a prototypical nonlinear delayed system. We have shown that the heat rate can be evaluated based on positional moments, implying that, despite the inherent memory, no temporal correlations are needed. Our formula can be combined with established approximations for the one-time PDF. We have presented analytical results for the Markovian limits and thereby observed a growing heat flow for vanishing delay time. We have further uncovered discontinuous behaviour at the onset of memory, which can also be found for underdamped motion when the delay appears in the velocity. This finite heat flow at small τ-values above and around the (ballistic) velocity-relaxation timescale, is a consequence of the interplay of noise and delay, while the discontinuity at τ → 0 is a consequence of the white noise assumption and the overdamped limit. For experimental realisations of feedback traps this theoretical result implies that the (unavoidable) delay causes a finite steady-state heat flow unless the delay is significantly smaller than the (ballistic) relaxation timescale.
Although the here presented investigations are restricted to static energy landscapes, we indeed expect the results to also hold for time-dependent external potentials, at least, if their changes are slow compared to the other dynamics. Feedback loops offer a way to precisely control the effective temperature of a particle by suppressing its thermal fluctuations 55 . When this technique is used to change temperature, e. g., during the cyclic process of a heat engine 56 , the here discussed heat flow, which is caused even by tiny delays, might play a non-negligible role and should therefore be taken into account in the theoretical calculations of efficiency and entropy production.
Moreover, we have presented an approach capturing jump processes, thereby predicting that the medium entropy production due to the feedback is maximal, when the delay-induced oscillations are coherence-resonant. This work is an important step towards an understanding of thermodynamic notions in delayed systems. Future work will focus on the (non-Gaussian) heat distributions P(q), which appear to violate fluctuation relations 29 and on developing approaches for the total entropy production, possibly via Markovian embedding techniques 57 . We hope that our current findings will stimulate experimental investigations on passive (or even active) colloidal systems, e. g., a validation of the predicted switching from feedback cooling to heating by adjusting the delay force strength, which is tunable in experimental setups.

Technical Aspects
Fokker-Planck equation. The Fokker-Planck equation of the considered non-Markovian system (1) is an infinite hierarchy of coupled equations 45,47,48 , whose first member reads Apart from the one-time PDF, ρ 1 , Eq. (24) also involves the two-time joint PDF, ρ 2 , and is hence not self-sufficient.

Derivation of relations between correlation functions (7,8).
In order to derive (7), we project the Fokker-Planck equation (24) onto its moments by multiplying it with x n+1 , n ≥ 0, and integrating over the spatial domain x ∈ (−∞,∞). Several partial integrations, (where the boundary terms vanish) yield In a NESS, the left side vanishes and (25) reduces to a relation between positional moments and the spatial autocorrelation function, as given in (7). On the other hand, one can deduce directly by plugging in the LE (1)  (26) and (25) readily provides the surprisingly simple and generic relation (8), which, in fact, generally holds for any (nonlinear) force (as can be shown analogously).

Derivation of noise-position (19) and noise-velocity cross correlations for underdamped dynamics.
We aim to calculate the correlations ) ( ) for underdamped motion with position-dependent feedback. To this end, we first we use again the causality argument 58 that the noise cannot influence past position or past velocity. This implies that both correlations must vanish at τ > 0. Second, we rewrite the correlations at τ = 0 as follows, starting from a formal integration over the LE (18)    by using the causality argument, which implies that the integrands are zero for all s < t The constant  amounts for the initial conditions. Third, using the identity t t s t 0 0 , Fourth, we establish that Eqs (28,29) only allows for finite solutions. This can be shown by contradiction, starting with the assumption ξ 〈 〉→ ∞  X t t ( ) ( ) . Eq. (29) then yields a finite value of ξ 〈 〉 X t t ( ) ( ) , which, on the other hand, yields a finite value of ξ 〈 〉  X t t ( ) ( ) from Eq. (29), thus a contradiction. Using this finiteness, we finally obtain for τ ≥ 0 0 Numerical simulations. We perform Brownian dynamics simulations with Euler-Maruyama integration scheme. The temporal discretisation is varied between Δt = 10 −6 and 10 −4 , such that the typical time scales are all properly resolved, in particular, the delay time τ > 1000Δt and the intrawell relaxation time t rel > 1000Δt. Furthermore, steady-state ensemble averages are obtained after appropriate transient times (at least 5 × 10 6 time steps) have been cut off, and by using a sufficiently high number of realisations (at least 10 4 ). Both is checked by insensitivity against further increase of simulated time or sample size. Numerical integrations are performed in Stratonovich calculus.
Reduction to 2-state model. Within the 2-state model, the particle being in one of the two potential wells is represented by one discrete state s = ±x 0 . The transition rates p j∈{1,2} between the states dependent on the delayed position, with p 1 if s(t)s(t − τ) > 0, and p 2 if s(t)s(t − τ) < 0. Within a quasistatic approximation (appropriate if τ  t rel and γD 0 < V 0 ), we consider the corresponding quasistatic potential V qs (x) = V s (x) + V d (x, x τ = ±x 0 ), obtained by fixing x τ at one of the two deterministically stable solutions ±x 0 . The transition rates can then be approximated by the two Kramers rates associated with V qs , given by (32) with p j = p(x min,j ) for escapes out of the two minima, x min,1 = ± x 0 and 2 , while the potential barrier heights read ΔV qs,1 = V 0 + k/2, and ΔV qs,2 = V 0 − 3k/2. This finally yields the transition rates