# Equilibrium Stochastic Delay Processes

Viktor Holubec1,2,⋆, Artem Ryabov1, Sarah A.M. Loos2,3, Klaus Kroy2 1 Charles
University, Faculty of Mathematics and Physics, Department of Macromolecular
Physics, V Holešovičkách 2, CZ-180 00 Praha, Czech Republic  2 Institut für
Theoretische Physik, Universität Leipzig, Postfach 100 920, D-04009 Leipzig,
Germany  3 International Centre for Theoretical Physics, Str. Costiera 11,
34151 Trieste, Italy  $⋆$<EMAIL_ADDRESS>
###### Abstract

Stochastic processes with temporal delay play an important role in science and
engineering whenever finite speeds of signal transmission and processing
occur. However, an exact mathematical analysis of their dynamics and
thermodynamics is available for linear models only. We introduce a class of
stochastic delay processes with nonlinear time-local forces and linear time-
delayed forces that obey fluctuation theorems and converge to a Boltzmann
equilibrium at long times. From the point of view of control theory, such
“equilibrium stochastic delay processes” are stable and energetically passive,
by construction. Computationally, they provide diverse exact constraints on
general nonlinear stochastic delay problems and can, in various situations,
serve as a starting point for their perturbative analysis. Physically, they
admit an interpretation in terms of an underdamped Brownian particle that is
either subjected to a time-local force in a non-Markovian thermal bath or to a
delayed feedback force in a Markovian thermal bath. We illustrate these
properties numerically for a setup familiar from feedback cooling and point
out experimental implications.

††: New J. Phys.

  *

Keywords: Stochastic delay differential equations, Exact solution,
Equilibrium, Fluctuation theorems, Linear response theory, Feedback cooling

##  1 Introduction

Consider the stochastic delay differential equations (SDDEs)

$\displaystyle\dot{x}(t)$ $\displaystyle=$ $\displaystyle v(t),$ (1)
$\displaystyle m\dot{v}(t)$ $\displaystyle=$ $\displaystyle F(t)+F_{\rm
D}(t-\tau)+\eta(t)$ (2)

with a nonlinear time-local force $F(t)=F(x,v,t)$ and a linear delay force
($\tau>0$)

$F_{\rm D}(t-\tau)=-\kappa_{\tau}x(t-\tau)-\gamma_{\tau}v(t-\tau),$ (3)

with constant coefficients $\kappa_{\tau}$ and $\gamma_{\tau}$. The dynamics
is randomly driven by a possibly non-Markovian, zero-mean Gaussian stochastic
noise $\eta(t)$. Intuitively, one can think of Eqs. (1) and (2) as describing
the time evolution of the position $x(t)$ and velocity $v(t)$ of a Brownian
particle with mass $m$ and driven by the combined forces $\eta$, $F$, and
$F_{D}$. These forces can arise from various origins, e.g., from the
environment and the experimental apparatus, including some specifically
tailored feedback mechanisms. Further specifications and various
interpretations will be provided below. Due to finite speeds of information
transfer and processing and elements with slow response, such equations are
ubiquitous in engineering [1], biology [2, 3, 4] and even economics [5, 6, 7,
8]. Most frequently, they are applied in modelling of feedback loops [9, 10,
11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25], neural networks
[26, 27, 28, 29], population dynamics [30, 31], and epidemiology [32, 33].

Rising interest in SDDEs among physicists [34] is driven by recent
experiments. In the so-called feedback cooling experiments with Brownian
particles, one employs a feedback of the particle’s past velocity to achieve a
more localised state [35, 36, 37, 38]. In the surging field of active matter
[39, 40, 41], inevitable time delays in the control of robotic swarms [18] led
to investigations of the stability and localization of many-body systems with
delayed interactions [22, 18, 21, 24, 42]. In agreement with engineering
practice [1, 10, 11, 43], it was found that delay generally introduces
instabilities and oscillations into the dynamics [22, 42] and increases
stability and localization only in special cases [40].

Similarly, inevitable instrumental and feedback delays in micro-manipulation
experiments [44, 22] used to test stochastic thermodynamics [45, 46] has
triggered investigation of the thermodynamic aspects of SDDEs [47, 48, 49, 50,
13]. There are interesting consequences of the acausality of time-reversed
processes in delay systems due to the tracking (future) history for the time-
reversal. If interpreted as feedback-driven systems with information inflow,
their total entropy production rate, $\dot{S}_{\rm tot}$, evaluated as a ratio
of forward to backward path probabilities, is not just the sum, $\dot{S}_{\rm
S}+\dot{S}_{\rm NM}$, of entropy fluxes into the system (S) and into the bath
(B) [51, 47, 48, 50]. This means that the second law $\dot{S}_{\rm tot}\geq 0$
does not imply positivity of $\dot{S}_{\rm S}+\dot{S}_{\rm NM}$. These results
are generic for the system (1)–(2) with a Gaussian white noise
$\eta(t)\propto\xi(t)$, $\left<\xi(t)\right>=0$,
$\left<\xi(t)\xi(t^{\prime})\right>=\delta(t-t^{\prime})$. However, explicit
expressions are currently only available for linear systems [51, 47, 48, 52,
50], which fail to describe a broad range of interesting effects observed in
presence of nonlinear forces [47, 48, 50]. The same can be said about the
probability densities for SDDEs. They are available only for simple linear
setups [53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 42], and nonlinear systems
have been treated by various approximate techniques [59, 60, 63, 18, 64, 42,
65].

Even without time delay, an exact treatment of nonlinear systems is indeed
difficult. However, their stationary and relaxation properties are known
exactly in thermodynamic equilibrium. In this work, we extend this property to
a certain class of SDDEs. Our results can be of interest not only to the
theory of delay processes but also in applied contexts, like in control
theory.

##  2 Main results

As our main result, we identify a class of nonlinear delay processes that
admit a standard thermodynamic description, including the second law
inequality $\dot{S}_{\rm tot}=\dot{S}_{\rm S}+\dot{S}_{\rm NM}\geq 0$. If not
driven, they obey Boltzmann statistics in the steady state. We therefore
characterize these processes as “equilibrium delay processes”. The key idea is
to accompany the time-delayed feedback force applied to the system with a
suitable colored noise $\eta_{\mathrm{FB}}$ and interpret the resulting
overall system as a particle immersed in an equilibrium reservoir and
controlled by time-local external forces. Noteworthy, such feedback noise can
already be realized in state-of-the-art experimental setups [66, 67, 68, 38].
We further point out how to interpret Eqs. (1) and (2) as a feedback-driven
system and how to apply our results therein. Altogether, we provide three
complementary interpretations for the same stochastic process: a special type
of system with time-delayed forces (Sec. 1), system with time-local forces and
a heat bath with memory (Sec. 2.1), and a feedback-driven system (Sec. 2.3).
They differ just in the interpretation of the individual forces on the right
hand side of Eq. (2). In Tab. 1, we summarize relations between the three
interpretations and the definitions of the corresponding forces. In the
following, we take Boltzmann’s constant $k_{\rm B}$ as our unit of entropy.

A: Delay system
---
$m\dot{v}=F+F_{\rm D}+\eta$
_time-local systematic force_ : $F=F_{\rm E}+\kappa_{\tau}x(t)-\gamma_{0}v(t)$
ordinary (non-feedback) external force: $F_{\rm E}=-\partial_{x}U(x,t)+F_{\rm
N}(x,v,t)$
potential component of $F_{\rm E}$: $-\partial_{x}U(x,t)$
non-potential component of $F_{\rm E}$: $F_{\rm N}(x,v,t)$
_time-delayed force_ : $F_{\rm
D}=-\kappa_{\tau}x(t-\tau)-\gamma_{\tau}v(t-\tau)$
_total coloured noise from the environment and experimental apparatus_ :
$\eta(t)$

B: System with non-Markovian heat bath and time-local control
---
$m\dot{v}=F_{\rm E}+F_{\rm F}+\eta$
_time-local external force_ : $F_{\rm E}$
_time-delayed non-Markovian bath friction_:
$F_{\rm
F}=\kappa_{\tau}x(t)-\gamma_{0}v(t)-\kappa_{\tau}x(t-\tau)-\gamma_{\tau}v(t-\tau)=\kappa_{\tau}x(t)-\gamma_{0}v(t)+F_{\rm
D}$
total force from the non-Markovian heat bath at temperature $T$: $F_{\rm
F}+\eta(t)$
_coloured noise from the non-Markovian bath: $\eta(t)$_
heat flux into the system from the non-Markovian bath: $\dot{Q}_{\rm
NM}=\left<(F_{\rm F}+\eta)\dot{x}\right>$

C: Feedback-driven system with Markovian heat bath
---
$m\dot{v}=F_{\rm E}+F_{\rm FB}-\gamma_{0}v(t)+\sqrt{2T_{0}\gamma_{0}}\xi(t)$
_time-local external force_ : $F_{\rm E}$
_feedback force_ :
$F_{\rm
FB}=\kappa_{\tau}x(t)-\kappa_{\tau}x(t-\tau)-\gamma_{\tau}v(t-\tau)=\kappa_{\tau}x(t)+F_{\rm
D}+\eta_{\rm FB}(t)=F_{\rm F}+\gamma_{0}v(t)$
non-Markovian noise exerted by the feedback loop: $\eta_{\rm
FB}(t)=\eta(t)-\sqrt{2T_{0}\gamma_{0}}\xi(t)$
total force from the Markovian bath at temperature $T_{0}$:
$-\gamma_{0}v(t)+\sqrt{2T_{0}\gamma_{0}}\xi(t)$
_time-local friction from the Markovian bath_: $-\gamma_{0}v(t)$
_white noise from the Markovian bath_: $\sqrt{2T_{0}\gamma_{0}}\xi(t)$
heat flux into the system from the Markovian bath: $\dot{Q}_{\rm
M}=\left<(-\gamma_{0}v+\eta)\dot{x}\right>$

Table 1:  Three interpretations of the delay Langevin equation (2) employed in
this paper and the corresponding forces. The forces in the three
interpretations yield the same change in the momentum $m\dot{v}$ and thus the
same stochastic process. The term $\kappa_{\tau}x(t)-\gamma_{0}v(t)$ in the
force $F$ in the interpretation A is introduced to facilitate the
reinterpretation of the delayed force $F_{\rm D}$ as part of the friction
force $F_{\rm F}$ in B. By “ordinary external forces” in A we mean time-local
forces arising from physical interactions and thus not applied via a feedback
loop.

###  2.1 Mapping to time-local control and non-Markovian heat bath (Tab. 1B)

In this section, we describe the reinterpretation of the delay system (Tab.
1A) as an equilibrium system with memory (Tab. 1B). Consider a delay system
described by Eqs. (1)–(2) with the time-local force

$F(t)=F_{\rm E}(x,v,t)+\kappa_{\tau}x(t)-\gamma_{0}v(t).$ (4)

The specific form of the terms proportional to the constants $\gamma_{0}>0$
and $\kappa_{\tau}$ facilitates the reinterpretation of the delay force
$F_{D}$ in Eq. (2) as part of a friction force, below. The remaining force in
Eq. (4),

$F_{\rm E}(x,v,t)=-\partial_{x}U(x,t)+F_{\rm N}(x,v,t)$ (5)

is an arbitrary time-local force applied by external agents. It is composed of
potential and non-potential components $-\partial_{x}U(x,t)$ and $F_{\rm
N}(x,v,t)$. To distinguish the force $F_{\rm E}$ from the force applied via
the feedback loop in the feedback interpretation of Eqs. (1)–(2) (Tab. 1C), we
call it the “ordinary” external force.

Equation (2) now assumes the form

$m\dot{v}(t)=F_{\rm E}(t)+F_{\rm F}(t)+\eta(t)$ (6)

with $F_{\rm
F}(t)\equiv\kappa_{\tau}[x(t)-x(t-\tau)]-\gamma_{0}v(t)-\gamma_{\tau}v(t-\tau).$
It resembles the dynamical equation for the velocity of a particle subjected
to an external force $F_{\rm E}$ and immersed in a viscoelastic solvent
exerting on the particle the overall force $F_{\rm F}+\eta$ with systematic
component (friction) $F_{\rm F}$, and stochastic component (noise) $\eta$.

Such noise and friction can be interpreted to arise from an ordinary
equilibrium heat bath, i.e., a many-body system with infinite heat capacity in
thermal equilibrium, with a somewhat peculiar memory that gives rise to an
“echo” in the noise and friction (Tab. 1B). Notably, for an equilibrium heat
bath with a friction force linear in the variables $x$ and $v$, such as
$F_{\rm F}$, the time-reversal symmetry of the underlying microscopic dynamics
implies that the friction and noise are interrelated by the so-called second
fluctuation-dissipation theorem or fluctuation-dissipation relation (FDR) [69,
70, 71, 72]

$\left<\eta(t)\eta(t^{\prime})\right>={T}\,\Gamma(|t-t^{\prime}|).$ (7)

Here $T$ denotes the temperature and $\Gamma(t)$ is the so called friction
kernel defined by the integral

$F_{\rm
F}(t)=-\int_{-\infty}^{t}dt^{\prime}\,\Gamma(t-t^{\prime})v(t^{\prime}).$ (8)

For a given friction $F_{\rm F}$, the FDR (7) might imply that the noise must
be complex valued. However, in order to admit its ordinary physical
interpretation and realisability in a lab, $\eta(t)$ is required to be a real-
valued function. This condition implies that its power spectrum must be non-
negative,

$S(\omega)=\int_{-\infty}^{\infty}dt\left<\eta(t)\eta(0)\right>\exp(-i\omega
t)\geq 0.$ (9)

For the system of Eqs. (1)–(6), the conditions (7) and (9) can be satisfied
for a certain range of model parameters only, see Secs. 3 and 4. In this
range, Eqs. (1) and (6) can be interpreted as describing a system with
internal Hamiltonian $H=U(x,t)+mv^{2}/2$ acted upon by a non-potential force
$F_{\rm N}$ and coupled to a non-Markovian “equilibrium bath” at temperature
$T$. Let us now review some general properties of this system.

###  2.2 Properties of the mapping (Tab. 1B)

_Average thermodynamics._ If the above equilibrium mapping holds, the system’s
thermodynamics obeys standard relations from classical [73] and stochastic
[45, 46] thermodynamics. Namely, the average entropy flux into the non-
Markovian heat bath at temperature $T$ is given by the Claussius equality

$\dot{S}_{\rm NM}=-\dot{Q}_{\rm NM}/{T}$ (10)

where $\dot{Q}_{\rm NM}=\left<(F_{\rm F}+\eta)\dot{x}\right>$ is the average
heat flux from the heat bath into the system. It can also be interpreted as
the work done by the bath on the system per unit time. Here and below we
employ Stratonovich calculus. The averages $\left<\bullet\right>$ should be
performed over many realizations of the stochastic process.

The average heat flux is related via the first law, $d\langle
H\rangle/dt=\dot{Q}_{\rm NM}+\dot{W}_{\rm E}$, to the average power input,
$\dot{W}_{\rm E}=\langle\partial U/\partial t+F_{\rm N}\dot{x}\rangle$, of the
system, due to external manipulations of the potential $U$ and the non-
potential force $F_{\rm N}$. The sum of the rate of change of the system
entropy, $\dot{S}_{\rm S}$, and the entropy influx $\dot{S}_{\rm NM}$ in Eq.
(10) is the total entropy production, which obeys the second law of
thermodynamics [73]:

$\dot{S}_{\rm tot}=\dot{S}_{\rm S}+\dot{S}_{\rm NM}\geq 0.$ (11)

_Dynamics._ Unlike a general delay system, which can exhibit over-damped,
damped oscillatory, but also diverging behavior [60, 61, 42, 13, 10, 11],
systems obeying the mapping of Sec. 2.1 always eventually relax into a time-
independent steady state for time independent parameters, confining potential
$U$, and stationary non-potential forces $F_{\rm N}$. If the latter vanishes
in Eq. (5), the stationary probability density function (PDF) for position and
velocity is given by the Gibbs canonical distribution,
$p(x,v;T)=p_{x}(x;{T})p_{v}(v;{T})$, with

$\displaystyle p_{x}(x;{T})$ $\displaystyle=$
$\displaystyle\exp[-U(x)/{T}]/Z_{x}({T}),$ (12) $\displaystyle p_{v}(v;T)$
$\displaystyle=$ $\displaystyle\exp(-mv^{2}/(2{T}))/Z_{v}({T}),$ (13)

normalized by $Z_{x}({T})=\int_{-\infty}^{\infty}dx\,\exp(-U(x)/{T})$ and
$Z_{v}({T})=\int_{-\infty}^{\infty}dv\,\exp(-mv^{2}/(2{T}))$. This is an
equilibrium steady state and thus the corresponding entropy production rates
$\dot{S}_{\rm S}$, $\dot{S}_{\rm NM}$, and $\dot{S}_{\rm tot}$ vanish. For
quasi-static variations of the potential, when the system PDF evolves through
a set of such states, the total entropy change $\Delta S_{\rm
tot}=\int_{0}^{t}\dot{S}_{\rm tot}(t^{\prime})dt^{\prime}$ vanishes and the
entropy change in the system, $\Delta S_{\rm S}$, is exactly balanced by the
entropy change in the bath, $\Delta S_{\rm NM}$. The relaxation process to
equilibrium is always accompanied by a decrease in the free energy of the
system. This functional thus represents the Lyapunov function for the
relaxation process that can be easily evaluated from stochastic trajectories
of the system. Even stronger restrictions on the relaxation dynamics are
imposed by the Evans–Searles fluctuation theorem [74, 75]. In contrast,
besides a limited success [76], it is currently unknown if similar general
restrictions also apply to relaxation towards non-equilibrium steady states.

The validity of these results for an arbitrary potential $U(x)$ follows from
general considerations of equilibrium statistical physics [73] and the FDR
[69, 70, 71, 72]. However, a closed dynamical equation, e.g., of Fokker-Planck
type [77, 78], for the PDF of a nonlinear delay process is not known [13]
making a general direct verification difficult. In Sec. 5, we provide an
explicit test for the specific potential $U(x,t)=k_{6}x^{6}/6+k_{3}x^{3}/3$
using Brownian dynamics (BD) simulations of Eqs. (1) and (2). Besides, we
tested the described results for various other polynomial potentials.

We stress that the described equilibrium-like properties of equilibrium delay
processes do not trivialize their dynamics. As an example, consider a
situation when the force $F$ in Eq. (2) is linear in $x$ and $v$ and thus the
system (1)-(2) is exactly solvable. For fixed initial conditions, one finds
that the average position $\langle x(t)\rangle$ and velocity $\langle
v(t)\rangle$ are identical for equilibrium ($\eta_{\rm FB}$ in Tab. 1
determined by the conditions (7) and (9)) and standard ($\eta_{\rm FB}=0$)
delay processes. The four correlation functions $\langle A(t)B(0)\rangle$ for
$A,B=x,v$ may then merely differ in the stationary distribution of the initial
conditions.

_Fluctuation theorems._ From a stochastic-thermodynamics perspective, it is
interesting to also consider a finite-speed protocol rendering the potential
time-dependent. Specifically, in Sec. 5, we test two fluctuation theorems for
the stochastic work $w=\int_{0}^{t}\partial U(x,t^{\prime})/\partial
t^{\prime}\,dt^{\prime}$ done on the system, if $k_{6}=k_{6}(t^{\prime})$,
$t^{\prime}\in(0,t)$ is varied non-quasi-statically, namely the Jarzynski
equality [79]

$\left<\exp(-w/{T})\right>=\exp(-\Delta F/{T})$ (14)

and the Crooks’ fluctuation theorem [80]

$\rho_{\rm F}(w)/\rho_{\rm R}(-w)=\exp[(w-\Delta F)/{T}].$ (15)

Here, $\Delta F$ is the free energy difference between equilibrium states
corresponding to the final and initial values of the potential, $\rho_{\rm F}$
is the probability distribution for work measured along the process when the
potential changes from $U(x,0)$ to $U(x,t)$, and $\rho_{\rm R}$ is the
probability distribution for work measured along the time-reversed process.
For the both fluctuation theorems, the forward process departs from
equilibrium. The validity of Jarzynski’s equality requires the existence of
initial and final Gibbs stationary states and Crooks’ fluctuation theorem
additionally requires the FDR and Gaussianity of the noise [81]. The described
processes fulfill all these requirements and, indeed, our simulations confirm
Eqs. (14) and (15).

_Perturbative expansions._ Even though based on an ad hoc choice of the noise
our results represent first exact analytical solutions to stationary PDFs for
a nonlinear SDDE. As such, they might pave the way for studying steady states
and thermodynamic properties of systems controlled by more general nonlinear
SDDEs. We show in Sec. 6 that linear-response theory [71] can be used to
calculate time-dependent averages in perturbed (nonlinear) equilibrium delay
systems. Besides such classical linear response, one can derive some explicit
approximate formulas for specific perturbations on the level of moments
calculated directly from the nonlinear system of SDDEs (1) and (2).

###  2.3 Equilibrium feedback (Tab. 1C)

The formal interpretation of dynamical equations (1)–(2) as a model for a
system immersed in a non-Markovian equilibrium bath and driven by a time-local
force $F_{\rm E}$, in Sec. 2.1, allowed us to utilize a wealth of known
results. However, in practice, these equations usually describe feedback-
driven systems in contact with a Markovian heat bath exerting a memoryless
friction $-\gamma_{0}v$ and Gaussian white noise
$\sqrt{2T_{0}\gamma_{0}}\xi(t)$ with
$\left<\xi(t)\xi(t^{\prime})\right>=\delta(t-t^{\prime})$. Usually the
system’s environment provides such a bath. To investigate this “more natural”
interpretation, we rewrite the dynamical equation for the velocity as

$m\dot{v}(t)=F_{\rm E}(x,v,t)+F_{\rm
FB}-\gamma_{0}v(t)+\sqrt{2T_{0}\gamma_{0}}\xi(t)$ (16)

and interpret it as describing a system immersed in a standard, i.e., Gaussian
and Markovian, heat bath at temperature $T_{0}$. This system is controlled by
the time-local force $F_{\rm E}$ and the feedback force

$F_{\rm FB}=\kappa_{\tau}[x(t)-x(t-\tau)]-\gamma_{\tau}v(t-\tau)+\eta_{\rm
FB}(t)$ (17)

composed of the systematic delayed component $F_{\rm F}$ and the “feedback”
noise $\eta_{\rm FB}(t)\equiv\eta(t)-\sqrt{2T_{0}\gamma_{0}}\xi(t)$, see Tab.
1C. Given that the conditions (7) and (9) are fulfilled, we call this process
an equilibrium feedback (EFB) process.

Importantly, the formal results concerning the system dynamics, i.e. the
stationary PDFs (12) and (13), are valid regardless of the interpretation, and
thus they apply also for EFB. This means that the EFB is ideal from the point
of view of passivity-based control [82], which is a branch of control theory
that aims to balance the power delivered into the system with its dissipation.
Generic feedback can lead to divergences and instabilities when the energy
influx by the feedback gradually increases the internal energy of the system.
However, EFB processes are always stable and passive in the sense that the
resulting steady states are robust against perturbations and all the energy
injected into the system is dissipated. In A, we moreover show that, under
realistic conditions, the temperature $T$ corresponding to the Boltzmann PDF
reached by the EFB is always larger than the ambient temperature $T_{0}$.

The thermodynamics of EFB has to be treated with care. In particular, the
total entropy production is interpretation-dependent. But the stochastic work
done on the system by varying the potential remains the same, and the
fluctuation theorems (14) and (15) are still valid. Differences arise in the
definitions of the remaining thermodynamic fluxes. With the present definition
of the heat bath, the heat flux reads $\dot{Q}_{\rm
M}=\langle(-\gamma_{0}v+\sqrt{2T_{0}\gamma_{0}}\xi)\dot{x}\rangle$. And, in
addition to the average power $\dot{W}_{\rm E}$ delivered to the system by the
potential and non-potential forces, one has to consider also the power
$\dot{W}_{\rm FB}=\langle F_{\rm FB}\dot{x}\rangle$ associated with the
feedback force $F_{\rm FB}$.

In a conventional feedback process, this power is accompanied by an
information influx [51, 47, 48, 50] that, for example, allows the feedback to
cool the system [35, 37]. The resulting (effective) temperature of the system
is then smaller than the temperature of the ambient bath, implying a positive
heat flux from the bath into the system, $\dot{Q}_{\rm M}>0$. In a steady
state, the conventional feedback is thus able to cool the ambient bath by
extracting the power $-\dot{W}_{\rm FB}=\dot{Q}_{\rm M}>0$ from it. However,
for an arbitrary force $F_{\rm E}$, the second law (11) together with the
relation $\dot{Q}_{\rm NM}=\dot{Q}_{\rm M}+\dot{W}_{\rm FB}=-T\dot{S}_{\rm
NM}$ imposes an upper bound $\dot{Q}_{\rm M}\leq T\dot{S}_{\rm S}-\dot{W}_{\rm
FB}$ on the heat delivered from the bath to the system via the EFB. And, in A,
we show that under equilibrium conditions, $\partial U/\partial t=F_{\rm
N}=\dot{S}_{\rm S}=0$, the EFB brings the system to an effective temperature,
$T$, larger than the ambient temperature, $T_{0}$. Hence, the heat flux
$\dot{Q}_{\rm M}$ is always negative, the EFB performs net work on the system,
$\dot{W}_{\rm FB}=-\dot{Q}_{\rm M}>0$, and it eventually heats the ambient
bath. This means that the EFB cannot be used for standard (zero non-potential
force and time-independent potential) feedback cooling of the system [35, 37].

Sections 3 and 4 clarify when EFB can be realized with time-delayed forces
depending on either the earlier position or velocity, i.e., when the
corresponding feedback noise $\eta_{\rm FB}$ in Tab. 1C can be constrained to
be real valued. The technical details are given in B. The resulting parameter
regimes where the EFB can be realized in these two situations are depicted in
phase diagrams (Figs. 1 and 2). Equilibrium feedback with time-delayed forces
depending on both delayed position and velocity can be investigated in a
similar manner, but the corresponding phase diagram becomes three-dimensional.
In Sec. 5, we verify the validity of our theoretical results by a Brownian
dynamics (BD) simulation of the equilibrium velocity feedback. In Sec. 6, we
discuss several perturbative expansions pushing the theory beyond the
parameter regime of the equilibrium delay processes. We conclude in Sec. 7.

##  3 Equilibrium position feedback

Let us now consider the situation of the position-dependent feedback force
($\gamma_{\tau}=0$ in Eq. (17))

$F_{\rm FB}=\kappa_{\tau}[x(t)-x(t-\tau)]+\eta_{\rm FB}(t)=F_{\rm
F}+\gamma_{0}v(t)+\eta_{\rm FB}(t).$ (18)

The generalized friction force $F_{\rm
F}=\kappa_{\tau}[x(t)-x(t-\tau)]-\gamma_{0}v(t)$ can be written using the
friction kernel

$\Gamma(t)=\left[2\gamma_{0}\delta(t)-\kappa_{\tau}\Theta(\tau-t)\right]\Theta(t),$
(19)

where $\Theta(.)$ denotes the Heaviside step function. This result can be
verified by direct substitution into Eq. (8) and integrating the term
including velocity $v(t)=\dot{x}(t)$ by parts, cf. Eq. (9.14) in Ref. [13].

The conditions (7) and (9) on the EFB imply that the friction kernel (19) and
the total noise $\eta(t)=\eta_{\rm FB}(t)+\sqrt{2T_{0}\gamma_{0}}\xi(t)$ (see
Tab. 1A) must obey the FDR,

$\left<\eta(t)\eta(t^{\prime})\right>/{T}=2\gamma_{0}\delta(t-t^{\prime})-{\kappa}_{\tau}\Theta(\tau-|t-t^{\prime}|),$
(20)

and that the corresponding power spectrum must be non-negative,

$S(\omega)=2\left[{\gamma}_{0}-{\kappa}_{\tau}\tau\frac{\sin(\omega\tau)}{\omega\tau}\right]\geq
0.$ (21)

Using $\max[\sin(x)/x]=1$ and $\min[\sin(x)/x]\approx-0.22$, this implies the
inequalities

$0\leq\max({\kappa}_{\tau},-0.22{\kappa}_{\tau})\tau\leq{\gamma}_{0},$ (22)

which specify the parameter regime where the noise $\eta(t)$ satisfying the
FDR (20) can actually be realized in the lab (for details of noise
realization, see B). The inequalities require non-negative $\gamma_{0}$ which
is always fulfilled in the EFB interpretation, where $\gamma_{0}$ measures the
strength of the background friction. For $\gamma_{0}\geq 0$, the inequalities
(22) bounds the feedback strength ${\kappa}_{\tau}$ as
$-\gamma_{0}/0.22\leq{\kappa}_{\tau}\leq\gamma_{0}$.

Stable,
standard FDR +Stable,
no standard FDR Unstable Figure 1: Phase diagram of the position feedback in
the reduced variables $\tilde{\kappa}_{\tau}=\kappa_{\tau}/m$ and
$\tilde{\gamma}_{0}=\gamma_{0}/m$. In the FDR region,
$\tilde{\gamma}_{0}\geq\max(\tilde{\kappa}_{\tau},-0.22\tilde{\kappa}_{\tau})\tau$
and the system has a positive relaxation time $t_{R}$. Then the system is
stable for arbitrary delay and it is possible to drive it by an equilibrium
position feedback. In the no-FDR region,
$\tilde{\gamma}_{0}<\max(\tilde{\kappa}_{\tau},-0.22\tilde{\kappa}_{\tau})\tau$
and $t_{\rm R}>0$, the system reaches a stable steady state, but the
equilibrium position feedback cannot be realized in practice. In the unstable
region ($t_{\rm R}<0$), the velocity exhibits exponentially diverging
oscillations due to large time delays and thus no steady state exists. For
$\tau=0$, the process is stable.

Under these conditions, the equilibrium position feedback fulfills all the
properties described in Sec. 2. In particular it eventually yields the stable
equilibrium distribution (12) and (13) whenever $\partial U/\partial t=F_{\rm
N}=0$. However, the time delay in a general feedback may yield diverging
trajectories for certain parameter values. As an independent check that the
parameter regime (22) allowing for equilibrium position feedback always leads
to stable stationary solutions, we investigate the overall stability of
position feedback described by Eqs. (16) and (18) for the case $F_{\rm E}=0$,
which can be inspected analytically.

Specifically, the process (16) eventually reaches a stable steady state if all
the corresponding relaxation times, $t_{\rm R}$, are positive. To calculate
them, we substitute the feedback force (18), $F_{\rm E}=0$ , and $v=\dot{x}$
in Eq. (16), set $\eta(t)=0$, and solve the resulting equation using the
exponential ansatz $x=\exp(-\lambda t/\tau)$ 111One can analogously treat
systems with a potential $U$, by linearising it around a (local) minimum and
absorbing the resulting linear time-local force into $\kappa_{\tau}$.. The
obtained transcendental equation

$m\lambda^{2}=\gamma_{0}\tau\lambda-\kappa_{\tau}\tau^{2}\left[\exp(\lambda)-1\right]$
(23)

can in general only be solved numerically and has infinitely many solutions.
As the relaxation time of the system, $t_{\rm R}$, we identify the smallest
$\tau/\Re[\lambda]$ solving Eq. (23), where $\Re[\bullet]$ denotes the real
part. The system eventually relaxes into a stable steady state with $\langle
v(t)\rangle=0$ if $t_{\rm R}>0$. An approximate explicit solution to Eq. (23)
can be obtained in the limit of small delay. Expanding the friction $F_{\rm
F}$ in Eq. (18) up to the first order in $\tau$, we get

$m\dot{v}(t)\approx-(\gamma_{0}-\kappa_{\tau}\tau)v(t)+\sqrt{2T_{0}\gamma_{0}}\xi(t).$
(24)

The last two terms can be interpreted as a noise and friction from an
equilibrium bath with friction coefficient $\gamma_{0}-\kappa_{\tau}\tau$,
which yields stable dynamics where all the energy injected into the system by
the feedback is dissipated (passive dynamics) if

$\kappa_{\tau}\tau\leq\gamma_{0}.$ (25)

The conditions for the general case are depicted in Fig. 1. Indeed, the whole
parameter regime where the EFB can be defined according to the FDR (20)
(green) is found to be stable. This shows that, as expected, EFB is passive
and stable. Nevertheless, the regime of stability, $t_{\rm R}>0$, is broader
(orange). Noteworthy, the system can be stable even for $\gamma_{0}<0$ if the
feedback strength $\kappa_{\tau}$ is also sufficiently negative. This could
have been anticipated from the approximate condition for stability (25), which
predicts the boundary between stable and unstable regimes for
$\tilde{\gamma}_{0}\tau\gtrapprox-2$ remarkably well. The approximate dynamics
allows to define an effective FDR with an effective temperature $T_{\rm
eff}=T_{0}/(1-\kappa_{\tau}\tau/\gamma_{0})$. However, beyond the small delay
approximation, the existence of such effective FDR is not guaranteed. For
smaller values of $\tilde{\gamma}_{0}$, higher order terms in the delay make
the system more unstable than expected from the linear analysis. In the
unstable regime (blue), the mean velocity exhibits exponentially increasing
oscillations [61, 42].

##  4 Equilibrium velocity feedback

Next, we perform the same analysis as in the previous section for the
velocity-dependent feedback force ($\kappa_{\tau}=0$ in Eq. (17))

$F_{\rm FB}=-\gamma_{\tau}v(t-\tau)+\eta_{\rm FB}=F_{\rm
F}+\gamma_{0}v(t)+\eta_{\rm FB}.$ (26)

The friction $F_{\rm F}=-\gamma_{\tau}v(t-\tau)-\gamma_{0}v(t)$ now
corresponds to the friction kernel

$\Gamma(t)=\left[2\gamma_{0}\delta(t)+\gamma_{\tau}\delta(t-\tau)\right]\Theta(t).$
(27)

in Eq. (8). The FDR relation (7) for the total noise $\eta(t)=\eta_{\rm
FB}(t)+\sqrt{2T_{0}\gamma_{0}}\xi(t)$ (see Tab. 1A) now reads

$\left<\eta(t)\eta(t^{\prime})\right>/{T}=2\gamma_{0}\delta(t-t^{\prime})+\gamma_{\tau}\delta(|t-t^{\prime}|-\tau),$
(28)

and thus the condition following from the positivity of the power spectrum (9)
reads

$S(\omega)=2\left[\gamma_{0}+\gamma_{\tau}\cos(\omega\tau)\right]\geq 0.$ (29)

The equilibrium velocity feedback thus can be realized if the inequality

$0\leq|\gamma_{\tau}|\leq\gamma_{0},$ (30)

holds (for detail of the realization, see B). As we have seen for the
equilibrium position feedback, $\gamma_{0}$ must be non-negative, which is
fulfilled in the EFB interpretation. For ${\gamma}_{0}\geq 0$, the
inequalities (30) impose that the amplitude $\gamma_{\tau}$ of the delayed
component of the friction can not exceed that of the Markov component.
Different from the corresponding inequality for $\kappa_{\tau}$ in the
position feedback, this condition is now symmetric with respect to
$\gamma_{\tau}=0$.

Similarly as in the case of the position feedback, we inspect the region of
stability of the general linear velocity feedback for $F_{\rm E}=0$ and
compare it to the region (30) allowing to realize the stable equilibrium
feedback. To this end, we insert the feedback force (26) and $F_{\rm E}=0$ in
Eq. (16), set $\eta(t)=0$, and solve the resulting equation using the
exponential ansatz $v(t)=\exp(-t/t_{\rm R}+i\omega t)$, with real parameters
$t_{\rm R}$ and $\omega$. Solving the resulting algebraic equation for the
relaxation time $t_{\rm R}$, we find

$t_{\rm R}=\frac{\tau}{\Re\left(\tilde{\gamma}_{0}\tau-{\rm
W}\left[-\tilde{\gamma}_{\tau}\tau\exp(\tilde{\gamma}_{0}\tau)\right]\right)},$
(31)

where ${\rm W}[.]$ stands for the Lambert ${\rm W}$ function, $\Re(.)$ denotes
the real part, and $\tilde{\gamma}_{0}=\gamma_{0}/m$ and
$\tilde{\gamma}_{\tau}=\gamma_{\tau}/m$. The Lambert ${\rm W}$ function is a
multivalued function and, in order to assess stability of the system, we
numerically determine the smallest $t_{\rm R}$ resulting from Eq. (31). In
this case, the small-delay expansion of the friction $F_{\rm F}$ in Eq. (26)
yields

$m\dot{v}(t)\approx-(\gamma_{0}+\gamma_{\tau})v(t)+\gamma_{\tau}\tau\dot{v}(t)+\sqrt{2T_{0}\gamma_{0}}\,\xi(t)$
(32)

and thus it suggest that the system will be stable for
$\gamma_{\tau}>-\gamma_{0}$ (positive effective friction coefficient) and
$\gamma_{\tau}\tau/m<1$ (positive effective mass). It also allows to define an
effective FDR with an effective temperature $T_{\rm
eff}=T_{0}/(1+\gamma_{\tau}/\gamma_{0})$ valid for small delays.

These formulas correctly yield the bottom boundary between the unstable and
stable regions in the phase diagram generated using the exact conditions (30)
and (31) depicted in Fig. 2. As for the position feedback, the region where
the FDR (20) and thus the equilibrium velocity feedback can be defined (green)
is indeed stable. And the regime of stability, $t_{\rm R}>0$, is broader than
the FDR regime and still extends to region of negative friction coefficients
$\gamma_{0}<0$ (orange). In the unstable regime (blue), the mean velocity
again exhibits exponentially increasing oscillations [61, 42].

Stable,
standard
FDR Stable,
no standard FDR Unstable Figure 2: Phase diagram of the velocity feedback in
the reduced variables $\tilde{\gamma}_{0}=\gamma_{0}/m$ and
$\tilde{\gamma}_{\tau}=\gamma_{\tau}/m$. In the FDR region,
$0\leq|\gamma_{\tau}|\leq\gamma_{0}$ and the system has a positive relaxation
time $t_{R}$. Then the system is stable for arbitrary delay and it is possible
to drive it by an equilibrium velocity feedback (EFB). In the no-FDR region,
$\gamma_{0}\leq|\gamma_{\tau}|$ and $t_{\rm R}>0$, the system reaches a stable
steady state but the EFB cannot be realized in practice. In the unstable
region ($t_{\rm R}<0$), the velocity exhibits exponentially diverging
oscillations due to large time delays and thus no steady state exists.

##  5 Demonstration of equilibrium velocity feedback

Let us now discuss a specific realization of the equilibrium velocity feedback
and show that it indeed has all the properties described in Sec. 2. As
detailed in B a possible (parsimonious) form of the total noise
$\eta(t)=\eta_{\rm FB}(t)+\sqrt{2T_{0}\gamma_{0}}\xi(t)$, which fulfills the
FDR (28) for equilibrium velocity feedback, is obtained by setting $\eta_{\rm
FB}(t)=\sqrt{\alpha_{\tau}}\xi(t-\tau)$. The parameters of the corresponding
feedback force $F_{\rm FB}(t)$ (26),

$F_{\rm FB}(t)=-\gamma_{\tau}v(t-\tau)+\sqrt{\alpha_{\tau}}\xi(t-\tau),$ (33)

can be tuned to represent various equilibrium velocity delay process. As a
benchmark for the EFB, we consider three processes distinguished by values of
the coefficients $\gamma_{\tau}$ and $\alpha_{\tau}$ above: (i) equilibrium
process (EQ) with $\gamma_{\tau}=\alpha_{\tau}=0$ and thus $F_{\rm FB}(t)=0$;
(ii) non-equilibrium (generic) velocity feedback (NEFB) with $\gamma_{\tau}>0$
and $\alpha_{\tau}=0$ and thus $F_{\rm FB}(t)=-\gamma_{\tau}v(t-\tau)$; and
(iii) equilibrium velocity feedback (EFB) with $\gamma_{0}\geq\gamma_{\tau}>0$
and $\alpha_{\tau}>0$ obeying Eq. (57). The last condition is compatible with
an equilibrium state at arbitrary temperature $T>T_{0}$, if we set

$\displaystyle\alpha_{\tau}$ $\displaystyle=$ $\displaystyle
2\gamma_{0}\left(T-T_{0}\right),$ (34) $\displaystyle\gamma_{\tau}$
$\displaystyle=$ $\displaystyle\pm
2\gamma_{0}\sqrt{\frac{T_{0}}{T}}\sqrt{1-\frac{T_{0}}{T}},$ (35)

where $T_{0}/T\leq 1$. Thus, in agreement with the discussion in A, the
additional noise present in the EFB always agitates or “heats” the system
above the ambient temperature $T_{0}$. Note that the above expressions do not
depend on the delay $\tau$.

###  5.1 Dynamics

To gain intuition concerning the behavior of the equilibrium velocity feedback
process, we now consider the specific system obeying Eqs. (1) and (16) with
the feedback force (33) and the force $F_{\rm E}=-\partial U/\partial x$
induced by the potential

$U(x)=\frac{k_{6}}{6}x^{6}+\frac{k_{3}}{3}x^{3}.$ (36)

We solve the dynamical equations using BD simulations for the NEFB, EFB and EQ
described above. In all our illustrations, we use $1/\tilde{\gamma}_{0}$ as
our time unit and $\sqrt{T_{0}/m}$ as our length unit. Velocity is thus
measured in units of $\tilde{\gamma}_{0}\sqrt{T_{0}/m}$. We show results from
BD simulations for the two parameter sets
$(\tilde{\gamma}_{0}\tau,\tilde{\gamma}_{\tau}\tau)\approx(0.28,0.28)$ and
$(\tilde{\gamma}_{0}\tau,\tilde{\gamma}_{\tau}\tau)\approx(0.42,0.24)$. The
first one yields fast relaxation of both feedback processes for $U=0$. The
second one is optimised to provide small velocity variance for EFB for $U=0$.
For more details, see C.

(a)(b) Figure 3: Stationary PDFs for $x$ (a) and $v$ (b) for the potential
(36) with $k_{6}=k_{3}=1$ in the parameter regime $\tau\approx 0.28$,
$\tilde{\gamma}_{\tau}\tau\approx 0.28$ (optimized for relaxation times of
feedback processes with $U=0$, see C). The data for the equilibrium process
(EQ) and equilibrium velocity feedback (EFB) perfectly agree with the
corresponding Boltzmann PDFs (12) and (13). Concerning the PDFs for generic
velocity feedback (NEFB), no exact analytical formula for the shown PDFs is
known. For all figures, simulation data was obtained from 50000 trajectories
of length 10 with time-step 0.001. (a)(b) Figure 4: The same as in Fig. 3 in
the parameter regime $\tau\approx 0.42$, $\tilde{\gamma}_{\tau}\tau\approx
0.24$ (optimized for velocity variance of EFB, see C). Other parameters are
the same as in Fig. 3. The insets magnify the regions around the global maxima
of the PDFs. The data for the equilibrium process (EQ) and equilibrium
velocity feedback (EFB) again perfectly agree with the corresponding Boltzmann
PDFs (12) and (13).  (a)(b)(c)(d)(e)(f)(g)(h) Figure 5: Relaxation dynamics of
first two central moments of velocity and position for NEFB, EFB and EQ
processes departing (with certainty) from the initial condition $x=0$, $v=1$
for $t\leq 0$. In (a-d) we show results for parameters used in Fig. (3) and in
(e-h) for those used in Fig. (4). The horizontal lines depict stationary
values of the shown moments obtained analytically using the Boltzmann
distributions (12) and (13) in Figs. 3 and 4.

In Figs. 3 and 4 we show the stationary PDFs for $x$ and $v$ obtained for the
first and second parameter set, respectively. In both figures, the simulated
PDFs for EQ and EFB perfectly overlap with the corresponding analytical
Boltzmann PDFs (12) and (13) providing numerical evidence for our claims in
Sec. 2. As expected, the position and velocity fluctuations are always
smallest for the NEFB and largest for the EFB.

To compare the relaxation dynamics of the three processes, we show in Fig. 5
the corresponding mean values $\left<x\right>$ and $\left<v\right>$ and
variances $\sigma_{x}^{2}$ and $\sigma_{v}^{2}$ as functions of time for the
initial condition $(v,x)=(1,0)$ for $t\leq 0$. Interestingly, the first
moments corresponding to the EFB (solid yellow line) relax faster than those
for the NEFB (dot-dashed blue line) and much faster than those for the EQ
process (broken green line). This is clearly a nonlinear effect because, for
$U=0$, the EFB and NEFB share the relaxation time (31). Especially for the EFB
the relaxation is considerably faster for the first parameter set [panels
(a)-(d)] than for the second choice [panels (e)-(h)]. This suggests that at
least some intuition gained from the linear regime $U=0$ also applies to the
nonlinear dynamics. In accord with Figs. 1 and 2, the position and velocity
fluctuations are always smallest for the NEFB and largest for the EFB. Smaller
velocity but also position variance for the EFB is obtained for the second
parameter set.

For the NEFB we were not able to analytically predict both the time evolution
of the depicted variables and their asymptotic values. To solve the full
transient dynamics for the EQ and EFB is also a difficult problem. However,
Fig. 5 shows that the moments in question converge to the values calculated
from the corresponding Boltzmann distributions (12) and (13) with temperatures
${T}_{0}$ (EQ) and ${T}>{T}_{0}$ (NEQ), which provides further numerical
evidence for our claims in Sec. 2.

###  5.2 Heat flux

Figure 6: Heat fluxes $\dot{Q}_{\rm M}$ from the (proper) Markovian bath into
the system for NEFB, EFB, and EQ processes discussed in Fig. 4 (e) – (h). The
horizontal lines depict the stationary values of the corresponding heat
fluxes. For the EQ process, the stationary value is 0. For the EFB process, it
is given by Eq. (37). For the NEFB, we have no universally valid prediction
for the stationary value of $\dot{Q}_{\rm M}$. However, the small-$\tau$
expansion (32) leads to $T_{\rm eff}<T_{0}$ and thus it suggests a positive
value of $\dot{Q}_{\rm M}$.

Let us now investigate the heat flux $\dot{Q}_{\rm
M}=\left<(-\gamma_{0}v+\sqrt{2T_{0}\gamma_{0}}\xi)\dot{x}\right>$ from the
“proper” Markovian bath into the system due to the feedback, see Tab. 1C and
Sec. 2.3. In A, we show that for a general EFB with a positive delay time
$\tau$ the heat flux always reads

$\dot{Q}_{\rm M}^{\rm EFB}=\gamma_{0}\left((\sigma_{v}^{\rm
EQ})^{2}-(\sigma_{v}^{\rm
EFB})^{2}\right)=\frac{2\gamma_{0}}{m}\left({T}_{0}-{T}\right)<0.$ (37)

The EFB thus always performs work on the system, which is eventually
dissipated in the bath. Figure 6 displays how the heat flux in the system
evolves during the relaxation processes for EQ, EFB and NEFB discussed in Fig.
5 (e)–(h). After the initial transient period, the heat flux for EFB converges
to the negative value given by Eq. (37) and thus it heats both the system, as
the corresponding stationary variances are larger than for the EQ process, and
the proper bath. For EQ, the stationary heat flux is zero as imposed by the
second law. For NEFB, the heat flux converges to a positive value. Thus the
NEFB cools the system while absorbing heat from the proper bath.

The result (37) applies for arbitrarily small positive delay $\tau$. The
specific form of the feedback force (33) allows us to also inspect what
happens for vanishing delay. Then the system is still in the Boltzmann
equilibrium state (12)–(13) with temperature ${T}$. However, the corresponding
total noise $\eta(t)=(\sqrt{2\gamma_{0}{T}_{0}}+\sqrt{\alpha_{\tau}})\xi(t)$
and friction $F_{\rm F}=-(\gamma_{0}+\gamma_{\tau})v(t)$ can now be
interpreted as a joint influence of the standard heat bath at temperature
$T_{0}$ and an additional ‘feedback heat bath’ at temperature $T_{\rm
F}=\alpha_{\tau}/2\gamma_{\tau}=\sqrt{{T}/{T}_{0}-1}\,{T}/2$. The laws of
thermodynamics imply that heat flows from hot to cold and thus $\dot{Q}_{\rm
M}$ is positive for $T_{\rm F}/T_{0}>1$ which occurs for $T>2T_{0}$. Further,
the heat flow is zero for $T=2T_{0}$, where $T_{\rm F}=T_{0}$ and thus there
is one global temperature only, and negative otherwise. Evaluating the heat
flux $\dot{Q}_{\rm M}^{\rm EFB}$ from Eqs. (16) and (33) with $\tau=0$ using
the approach of A, we find the expression

$\dot{Q}_{\rm M}^{\rm
EFB}=\frac{\gamma_{0}T_{0}}{m}\sqrt{\frac{T}{T_{0}}-1}\left(1-\sqrt{\frac{T}{T_{0}}-1}\right)$
(38)

which indeed obeys the described properties.

Since $\dot{Q}_{\rm M}^{\rm EFB}$ is strictly negative for $\tau>0$ and can be
both positive and negative for $\tau=0$, it exhibits a discontinuity at
vanishing delay, in accord with the results described in Ref. [50]. Note that
the presented situation with $\tau=0$ is physically weird since it seems
impossible to record the noise and feed it back into the system without any
delay. It also yields a strange behavior as the heat flux vanishes at the
point where temperatures $T_{\rm F}$ and $T_{0}$ are equal but $T=2T_{0}$.
This means that we constructed a bath at temperature $T>T_{0}$ by using two
strictly identical reservoirs at same the temperature $T_{0}$ to which the
system couples via different friction coefficients. The two baths provide the
same realizations of the white noise $\xi(t)$, and thus the total noise
intensity is given by the sum $\sqrt{2\gamma_{0}T_{0}}+\sqrt{{\alpha_{\tau}}}$
of the intensities of the two noises. In contrast, connecting a system to two
standard heat reservoirs always leads to equilibrium (vanishing heat flux)
when the temperatures of the two baths are equal. The mathematical reason is
that different reservoirs necessarily correspond to different noise
realizations, regardless of their temperatures. As an example, consider heat
baths A and B with friction and noise forces given by $-\gamma_{\rm A,B}v$,
$\sqrt{2\gamma_{\rm A,B}T_{0}}\xi_{\rm A,B}(t)$, with independent Gaussian
white noises $\xi_{\rm A,B}(t)$. Then the joint action of these baths is
described by the total friction $-(\gamma_{\rm A}+\gamma_{\rm B})v$ and noise
$\sqrt{2(\gamma_{\rm A}+\gamma_{\rm B})T_{0}}\xi(t)$, where $\xi(t)$ is a unit
variance Gaussian white noise (correlated with $\xi_{\rm A,B}(t)$). To sum up,
the formal identification of the feedback force $F_{\rm FB}$ and noise
$\eta_{\rm FB}$ for $\tau=0$ as effects of a standard heat bath correctly
determines the sign of the heat flux $\dot{Q}_{\rm M}$ in Eq. (38), but it is
physically problematic.

###  5.3 Fluctuation theorems

We conclude the numerical part of the paper by testing the work fluctuation
theorems (14) and (15). To this end, we let the system relax into the steady
state corresponding to the parameter regime of Fig. 4 (e) – (h) and then we
switch on the time-symmetric protocol

$k_{6}=1+0.9\sin(\pi t)$ (39)

with $t\in(0,1)$ for the potential (36). During the time-dependent driving, we
measure the stochastic work

$w=\int_{0}^{1}dt\,\frac{\partial{U}[x(t),t]}{\partial
t}=\int_{0}^{1}dt\,\dot{k}_{6}x^{6}(t)/6$ (40)

and sample its PDF $\rho(w)$. Due to the symmetry of the protocol, the time-
reversed process (R) and the forward process (F) in the fluctuation theorems
(14) and (15) coincide and the free energy difference $\Delta F$ vanishes.
Validity of the Crooks fluctuation theorem (15) for the acquired PDFs thus
implies that

$C(w)\equiv\log\left[\frac{\rho(w)}{\rho(-w)}\exp(w/{T})\right]=0.$ (41)

In Fig. 7 (a) we show the resulting PDFs for work obtained for the NEFB, EFB,
and EQ processes. The panels (b) – (d) then show that from the three processes
only the EFB (c) and EQ (d) yield $C(w)=0$ and thus fulfill the Crooks
fluctuation theorem (41).

The validity of the Jarzynski equality is tested in Fig. 8, where we show
values of averages $\left<\exp(-\beta_{\rm X}w)\right>$ over the sample PDFs
for work as functions of the parameter $\beta_{\rm X}$. For EFB and EQ, we
find that the average equals to one for $\beta_{\rm X}={1/T}$ and $1/{T}_{0}$,
respectively, proving the validity of Jarzynski equality (14) with $\Delta
F=0$ in these cases. For the NEFB, the system starts out of equilibrium so
that it is not clear which (inverse) temperature should be used in Eq. (14).
In the figure, we at least tested that choosing the temperature obtained form
the variance of the velocity, $1/{\beta}_{\rm X}=2\sigma_{v}^{2}$, does not
yield $\left<\exp(-\beta_{\rm X}w)\right>=1$.

(a)(b)(c)(d) Figure 7: Test of the Crooks fluctuation theorem. (a) Probability
densities for work for NEFB, EFB and EQ processes departing from the
stationary state of Fig. 4 and driven with a time symmetric protocol (39) for
the potential (36). The remaining panels show the function (41) for NEFB (b),
EFB (c), and EQ (d). The data suggest that the theorem holds for EFB and EQ,
where $C(w)\approx 0$. The shown results were obtained from $5\times 10^{6}$
runs of BD simulation with time-step $dt=10^{-3}$. Except for the time-
dependent driving, all parameters are the same as in Fig. 4. Figure 8: Test of
the Jarzynski equality for work PDFs from Fig. 7. The vertical lines for the
individual processes correspond to temperatures evaluated from the stationary
variance of the velocity as $1/\beta_{X}=m\left<v^{2}\right>$ assuming
equipartition theorem. For the EFB and EQ this temperature equals to the
temperature measured in any other way (e.g., from position or velocity PDF).
The NEFB in general induces a non-equilibrium steady state and thus
corresponding temperatures measured in different ways are in general different
[83, 84].

##  6 Beyond equilibrium feedback

In this section, we discuss possible analytical extensions of the equilibrium
delay processes that might help to better understand general delay processes.

###  6.1 Classical linear response theory

Any Langevin equation where the friction and noise obey the FDR (7) can be
thought of as a result of coarse-graining the full set of Hamiltonian
equations for the system of interest and the corresponding bath over the bath
degrees of freedom. This means that, the system with equilibrium delay can be
regarded as a standard Hamiltonian system, which implies applicability of the
classical linear response theory [69, 72, 71]. It states that the time-
evolution of the mean value $\left<A(t)\right>$ induced by perturbations of
the equilibrium system with Hamiltonian $H=U+mv^{2}/2$ in the form
$H+\varepsilon f(t)B$ starting at time $t=0$ reads [71]

$\left<A(t)\right>_{1}=\left<A(t)\right>_{0}+\frac{\varepsilon}{T}\int_{0}^{t}ds\,f(s)\left<A(t-s)\dot{B}(0)\right>_{0}.$
(42)

One assumes that the averages in the perturbed system can be expanded as
$\left<\dots\right>=\left<\dots\right>_{0}+\varepsilon\left<\dots\right>_{1}+\dots$,
where the subscript 0 denotes average taken over the unperturbed Boltzmann PDF
corresponding to Hamiltonian $H$ (12) and (13), the subscript 1 denotes
averages taken over the exact PDF up to the order $\varepsilon$, and so on.

We test the linear response theory using the specific equilibrium velocity
feedback system discussed in Sec. 5. We perturb the Hamiltonian by the term
$\varepsilon f(t)x$. This term corresponds to a homogeneous time-dependent
force $-\varepsilon f(t)$ and thus the dynamical equation for velocity reads

$m\dot{v}(t)=-\frac{\partial{U}}{\partial{x}}+F_{\rm
FB}-\gamma_{0}v(t)+\sqrt{2T_{0}\gamma_{0}}\xi(t)-\varepsilon f(t).$ (43)

The potential $U$ is given by Eq. (36) and the feedback force by Eq. (33). We
consider the parameter regime of Fig. 4 and the specific perturbation

$\displaystyle\varepsilon f(t)=\varepsilon\sin(\pi t/10).$ (44)

In Figs. 9 (a)-(d) we show the time correlation functions
$\left<A(t)\dot{x}(0)\right>_{0}=\left<A(t)v(0)\right>_{0}$,
$A=x,v,x^{2},v^{2}$ obtained using BD simulations of this system with
$\varepsilon=0$. The first two central moments of velocity and position
obtained using Eq. (42) with force (44) via these time correlation functions
are depicted in Figs. 9 (e)-(h) together with the corresponding quantities
obtained from BD simulation of the perturbed system. The figures show good
agreement between simulations and linear response theory, which improves for
smaller $\varepsilon$.

The validity of the linear response theory can be rationalized as follows.
Even though we do not have an exact dynamical equation for the PDF for the
equilibrium delay system, we know that the PDF for system and bath obeys a
Liouville equation. The corresponding Liouville operator is composed of the
system Hamiltonian $H$, the bath Hamiltonian, and the system-bath interaction
energy. Even though it is hard to identify the latter two, one can rely on
this Liouville equation as a starting point for perturbation theories around
the parameter regime of the equilibrium delay. For example, one can derive Eq.
(42) using the textbook approach of Ref. [71].

(a)(b)(c)(d)(e)(f)(g)(h) Figure 9: Linear response theory: (a)-(d) response
functions from BD simulations used to evaluate the time evolution of the first
two central moments of velocity and position (e)-(h) in the equilibrium system
of Fig. 4 perturbed by the time-dependent force (44). As expected, agreement
between the approximate theory (dashed lines) and simulation (noisy solid
lines) is better for smaller $\epsilon=0.1$ (black and gray lines) than for
the larger one $\epsilon=0.2$ (pink and orange lines).

###  6.2 Langevin equation

The classical linear response theory (42) applies only for perturbations that
can be subsumed into the Hamiltonian of the system. Other perturbations can be
treated, e.g., on the footing of linear irreversible thermodynamics [73] or
directly on the level of the Langevin equations (1) and (2). In order to
present two simple examples of the latter type, we write these equations in
the form of the equilibrium interpretation of Tab. 1B

$\dot{x}=v,\quad m\dot{v}=\left(-\frac{\partial U}{\partial x}+F_{\rm
F}+\eta\right)+\varepsilon g,$ (45)

where the term proportional to $\varepsilon$ is a perturbation. Perturbations
dependent on time and/or time-delayed variables again require evaluation of
time-correlation functions, which can rarely be obtained analytically.
Perturbations that depend only on position can be absorbed into the potential,
and the stationary PDF can be evaluated exactly. To obtain non-trivial
analytical results, we will investigate two properties of steady states
induced by perturbations of the form $g=g[v(t)]$, i.e. which depend solely on
velocity. However, the obtained general restrictions (47) and (48) on the
system dynamics apply for an arbitrary function $g$. In particular, $g$ can,
for these expressions, be a nonlinear function of position and velocity in the
past. In such a case, the resulting process (45) is a truly nonlinear delay
differential equation.

_Virial theorem:_ The viral theorem states that twice the average kinetic
energy of a system equals to the virial $-\left<Fx\right>$ [85]. For the
system at hand, the total force $F$ is given by the right-hand side (R.H.S.)
of Eq. (45). For the unperturbed system, this implies that

$m\left<v^{2}\right>=-\left<\left(-\frac{\partial U}{\partial x}+F_{\rm
F}+\eta\right)x\right>$ (46)

as can be checked directly from the Langevin equations after multiplying Eq.
(1) by $x$, Eq. (45) with $\varepsilon=0$ by $v$, summing the result, and
taking the stationary ensemble average so that the time derivative of the
cross-correlation $\left<xv\right>$ vanishes. Repeating this procedure for
nonzero perturbation in Eq. (45), we find that if $\left<g(v)\right>_{0}=0$
the unperturbed virtial theorem (46) remains valid up to first order in
$\varepsilon$ because
$\left<xg(v)\right>_{0}=\left<x\right>_{0}\left<g(v)\right>_{0}$. More
generally, we find that

$m\left<v^{2}\right>_{n}+\left<\left(-\frac{\partial U}{\partial x}+F_{\rm
F}+\eta\right)x\right>_{n}=-\varepsilon\left<xg(v)\right>_{n-1}$ (47)

holds for all corrections of order $n\geq 1$. Even though this result cannot
be evaluated explicitly for $n>1$ since
$\varepsilon^{n}\left<xg(v)\right>_{1}$ is unknown, it can provide a stringent
consistency check for simulations.

_Power:_ The power $\left<Fv\right>$ exerted by the total force $F$ on the
R.H.S. of Eq. (45) vanishes in the steady state of the unperturbed system
since the time-derivative of $\left<v^{2}\right>$ vanishes. With the
perturbation switched on, we find

$\left<\left(-\frac{\partial U}{\partial x}+F_{\rm
F}+\eta\right)v\right>_{n}=-\varepsilon\left<vg(v)\right>_{n-1}.$ (48)

Besides providing another set of expressions useful as consistency checks in
simulations, this equation provides an explicit non-trivial result for $n=1$.
Then the R.H.S. $-\varepsilon\left<vg(v)\right>_{0}$ is in general non-zero
and can be evaluated as average over the PDF (13). For example, for $g(v)=v$,
we find

$\dot{\tilde{W}}\equiv\left<\left(-\frac{\partial U}{\partial x}+F_{\rm
F}+\eta\right)v\right>_{n}=-\varepsilon T.$ (49)

Figure 10 verifies this approximate result by BD simulations for $n=1$ and
$n=3$ and the same unperturbed dynamics as in Sec. 6.1 above.

(a)(b) Figure 10: Test of the explicit result (49) (dashed lines) against the
simulation (solid line). (a) $\varepsilon=0.2$ and $n=1$. (b)
$\varepsilon=0.01$ and $n=3$. Other parameters are the same as in Fig. 4.

##  7 Conclusion

Using the second FDR, we have identified situations where one can interpret
time-delayed feedback forces proportional to velocity or position in Langevin
equations as friction forces imposed by an equilibrium bath. Our analysis
reveals a previously unnoticed class of nonlinear stochastic delay
differential equations, which can be treated analytically. They describe
processes that obey standard thermodynamic constraints. In particular, their
long-time distributions are of Gibbs canonical form and they obey standard
fluctuation theorems. From the point of view of control theory, especially
passivity-based control, the corresponding equilibrium feedback processes are
automatically stable and passive. However, their dynamics retains the full
complexity of generic delay processes. One disadvantage is that the
equilibrium feedback always heats up the system and thus, unlike generic
feedback protocols, cannot be used for standard feedback cooling.

As a practical demonstration, we have realized the equilibrium velocity
feedback using Brownian dynamics simulations and shown that it exhibits all
the formally derived properties. For so-called equilibrium velocity feedback,
not only the velocity at an earlier time but also the noise at that time are
used to drive the system at present. Such feedback can nowadays be realized in
practice in feedback experiments [12] with Brownian particles [14, 15, 16, 22,
23, 25, 38] or robots [18, 21, 24]. In particular, the so-called velocity
damping protocols for feedback cooling of trapped microscopic particles record
the particle velocity and later, after an experimental latency, apply to the
particle a force proportional to that velocity [68, 38]. Assuming that the
position and velocity dependence of the systematic force in the dynamical
equation for particle motion is known, measuring both its velocity and
position at time $t$ allows to determine the thermal noise, which can then be
applied to the particle in the future, similarly as the velocity-dependent
force. The potential drawback of this approach is that the measured noise will
be affected by measurement uncertainties and finite measurement time
resolution. The latter means that the obtained noise will effectively be
integrated (low-pass filtered) over one measurement frame, which is however
also the case in our BD simulations. Other promising setups, where the
equilibrium feedback might be realised are state-of-the-art bath engineering
experiments [67]. Finally, equilibrium feedback of the same type as in our
computer simulations can be realised in experimental setups where artificial
noise completely overshadows thermal noise. An example is shaken granular
matter [66], where the noise is realised by shaking the granular system.

From a theoretical perspective, we believe that our results will shed further
light on investigations of the dynamics and thermodynamics of nonlinear
stochastic delay differential equations, which are known to be immensely
resistant to analytical treatments. For example, the known stationary
distributions for the equilibrium delay processes can serve as starting points
of new perturbation theories valid for arbitrarily large delays. And, the
known thermodynamics of the equilibrium delay processes can help to better
understand the individual contributions to the total entropy production
derived for nonlinear stochastic delay differential equations as studied in
Refs. [51, 47, 48].

## Acknowledgments

We acknowledge funding through a DFG-GACR cooperation by the Deutsche
Forschungsgemeinschaft (DFG Project No 432421051) and by the Czech Science
Foundation (GACR Project No 20-02955J). VH gratefully acknowledges support by
Humboldt foundation.

## References

## References

  * [1] Kyrychko Y N and Hogan S J 2010 J. Vib. Control 16 943–960 ISSN 1077-5463 URL https://doi.org/10.1177/1077546309341100
  * [2] Beuter A, Bélair J and Labrie C 1993 Bull. Math. Biol. 55 525–541 ISSN 1522-9602 URL https://doi.org/10.1007/BF02460649
  * [3] Chen Y, Ding M and Kelso J A S 1997 Phys. Rev. Lett. 79(22) 4501–4504 URL https://link.aps.org/doi/10.1103/PhysRevLett.79.4501
  * [4] Novák B and Tyson J J 2008 Nat. Rev. Mol. Cell Biol. 9 981–991 ISSN 1471-0080 URL https://doi.org/10.1038/nrm2530
  * [5] Mackey M C 1989 J. Econ. Theory 48 497 – 509 ISSN 0022-0531 URL http://www.sciencedirect.com/science/article/pii/0022053189900392
  * [6] Voss H U and Kurths J 2002 Analysis of Economic Delayed-Feedback Dynamics (Boston, MA: Springer US) pp 327–349 ISBN 978-1-4615-0931-8 URL https://doi.org/10.1007/978-1-4615-0931-8_16
  * [7] Stoica G 2005 Proc. Am. Math. Soc. 133 1837–1841 URL https://www.jstor.org/stable/4097726
  * [8] Gao Q and Ma J 2009 Nonlinear Dyn. 58 209 ISSN 1573-269X URL https://doi.org/10.1007/s11071-009-9472-5
  * [9] Bechhoefer J 2005 Rev. Mod. Phys. 77(3) 783–836 URL https://link.aps.org/doi/10.1103/RevModPhys.77.783
  * [10] Atay F 2010 Complex Time-Delay Systems: Theory and Applications Understanding Complex Systems (Springer Berlin Heidelberg) ISBN 9783642023286 URL https://books.google.de/books?id=e2SGvHLpursC
  * [11] Lakshmanan M and Senthilkumar D 2011 Dynamics of Nonlinear Time-Delay Systems Springer Series in Synergetics (Springer Berlin Heidelberg) ISBN 9783642149382 URL https://books.google.de/books?id=nV5F0kMsqa4C
  * [12] Gernert R, Loos S A M, Lichtner K and Klapp S H L 2016 Feedback Control of Colloidal Transport (Cham: Springer International Publishing) pp 375–392 ISBN 978-3-319-28028-8 URL https://doi.org/10.1007/978-3-319-28028-8_19
  * [13] Loos S 2020 Stochastic systems with time delay – probabilistic and thermodynamic descriptions of non-Markovian processes far from equilibrium Ph.D. thesis TU Berlin
  * [14] Baraban L, Streubel R, Makarov D, Han L, Karnaushenko D, Schmidt O G and Cuniberti G 2013 ACS Nano 7 1360–1367 ISSN 1936-0851 URL https://doi.org/10.1021/nn305726m
  * [15] Qian B, Montiel D, Bregulla A, Cichos F and Yang H 2013 Chem. Sci. 4(4) 1420–1429 URL http://dx.doi.org/10.1039/C2SC21263C
  * [16] Bregulla A P, Yang H and Cichos F 2014 ACS Nano 8 6542–6550 ISSN 1936-0851 URL https://doi.org/10.1021/nn501568e
  * [17] Vásárhelyi G, Virágh C, Somorjai G, Tarcai N, Szörényi T, Nepusz T and Vicsek T 2014 Outdoor flocking and formation flight with autonomous aerial robots 2014 IEEE/RSJ International Conference on Intelligent Robots and Systems pp 3866–3873
  * [18] Mijalkov M, McDaniel A, Wehr J and Volpe G 2016 Phys. Rev. X 6(1) 011008 URL https://link.aps.org/doi/10.1103/PhysRevX.6.011008
  * [19] Zheng J, Dong J and Xie L 2017 IEEE Trans. Automat. Contr. 62 5866–5872
  * [20] Zhang J, xi Liu Y, Wu R B, Jacobs K and Nori F 2017 Physics Reports 679 1–60 ISSN 0370-1573 quantum feedback: theory, experiments, and applications URL https://www.sciencedirect.com/science/article/pii/S0370157317300479
  * [21] Leyman M, Ogemark F, Wehr J and Volpe G 2018 Phys. Rev. E 98(5) 052606 URL https://link.aps.org/doi/10.1103/PhysRevE.98.052606
  * [22] Khadka U, Holubec V, Yang H and Cichos F 2018 Nat. Commun. 9 3864 ISSN 2041-1723 URL https://doi.org/10.1038/s41467-018-06445-1
  * [23] Lavergne F A, Wendehenne H, Bäuerle T and Bechinger C 2019 Science 364 70–74 ISSN 0036-8075 URL https://science.sciencemag.org/content/364/6435/70
  * [24] Piwowarczyk R, Selin M, Ihle T and Volpe G 2019 Phys. Rev. E 100(1) 012607 URL https://link.aps.org/doi/10.1103/PhysRevE.100.012607
  * [25] Bäuerle T, Löffler R C and Bechinger C 2020 Nat. Commun. 11 2547 ISSN 2041-1723 URL https://doi.org/10.1038/s41467-020-16161-4
  * [26] Foss J, Longtin A, Mensour B and Milton J 1996 Phys. Rev. Lett. 76(4) 708–711 URL https://link.aps.org/doi/10.1103/PhysRevLett.76.708
  * [27] Marcus C M and Westervelt R M 1989 Phys. Rev. A 39(1) 347–359 URL https://link.aps.org/doi/10.1103/PhysRevA.39.347
  * [28] Sompolinsky H, Golomb D and Kleinfeld D 1991 Phys. Rev. A 43(12) 6990–7011 URL https://link.aps.org/doi/10.1103/PhysRevA.43.6990
  * [29] Haken H 2007 Brain Dynamics: An Introduction to Models and Simulations Springer Series in Synergetics (Springer Berlin Heidelberg) ISBN 9783540752387 URL https://books.google.de/books?id=H4amAIqgLGMC
  * [30] Gopalsamy K 2013 Stability and Oscillations in Delay Differential Equations of Population Dynamics Mathematics and Its Applications (Springer Netherlands) ISBN 9789401579209 URL https://books.google.de/books?id=0TfpCAAAQBAJ
  * [31] Mao X, Yuan C and Zou J 2005 J. Math. Anal. Appl. 304 296 – 320 ISSN 0022-247X URL http://www.sciencedirect.com/science/article/pii/S0022247X04007668
  * [32] Beretta E, Hara T, Ma W and Takeuchi Y 2001 Nonlinear Analysis: Theory, Methods & Applications 47 4107–4115 ISSN 0362-546X proceedings of the Third World Congress of Nonlinear Analysts URL https://www.sciencedirect.com/science/article/pii/S0362546X01005284
  * [33] Rihan F A, Alsakaji H J and Rajivganthi C 2020 Advances in Difference Equations 2020 502 ISSN 1687-1847 URL https://doi.org/10.1186/s13662-020-02964-8
  * [34] Otto A, Just W and Radons G 2019 Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences 377 20180389 URL https://royalsocietypublishing.org/doi/abs/10.1098/rsta.2018.0389
  * [35] Bushev P, Rotter D, Wilson A, Dubin F m c, Becher C, Eschner J, Blatt R, Steixner V, Rabl P and Zoller P 2006 Phys. Rev. Lett. 96(4) 043003 URL https://link.aps.org/doi/10.1103/PhysRevLett.96.043003
  * [36] Li T 2013 Millikelvin cooling of an optically trapped microsphere in vacuum Fundamental Tests of Physics with Optically Trapped Microspheres (Springer) pp 81–110
  * [37] Goldwater D, Stickler B A, Martinetz L, Northup T E, Hornberger K and Millen J 2019 Quantum Science and Technology 4 024003 URL https://doi.org/10.1088/2058-9565/aaf5f3
  * [38] Penny T W, Pontin A and Barker P F 2021 Phys. Rev. A 104(2) 023502 URL https://link.aps.org/doi/10.1103/PhysRevA.104.023502
  * [39] Ramaswamy S 2010 Ann. Rev. Cond. Matt. Phys. 1 323–345 URL https://doi.org/10.1146/annurev-conmatphys-070909-104101
  * [40] Bechinger C, Di Leonardo R, Löwen H, Reichhardt C, Volpe G and Volpe G 2016 Rev. Mod. Phys. 88(4) 045006 URL https://link.aps.org/doi/10.1103/RevModPhys.88.045006
  * [41] Gompper G, Winkler R G, Speck T, Solon A, Nardini C, Peruani F, Löwen H, Golestanian R, Kaupp U B, Alvarez L, Kiørboe T, Lauga E, Poon W C K, DeSimone A, Muiños-Landin S, Fischer A, Söker N A, Cichos F, Kapral R, Gaspard P, Ripoll M, Sagues F, Doostmohammadi A, Yeomans J M, Aranson I S, Bechinger C, Stark H, Hemelrijk C K, Nedelec F J, Sarkar T, Aryaksama T, Lacroix M, Duclos G, Yashunsky V, Silberzan P, Arroyo M and Kale S 2020 J. Phys. Condens. Matter 32 193001 URL https://doi.org/10.1088%2F1361-648x%2Fab6348
  * [42] Geiss D, Kroy K and Holubec V 2019 New J. Phys. 21 093014 URL https://iopscience.iop.org/article/10.1088/1367-2630/ab3d76
  * [43] Fridman E and Shaikhet L 2016 Automatica 74 288 – 296 ISSN 0005-1098 URL http://www.sciencedirect.com/science/article/pii/S0005109816303028
  * [44] Jun Y and Bechhoefer J 2012 Phys. Rev. E 86(6) 061106 URL https://link.aps.org/doi/10.1103/PhysRevE.86.061106
  * [45] Sekimoto K 2010 Stochastic Energetics Lecture Notes in Physics (Springer Berlin Heidelberg) ISBN 9783642054112 URL https://books.google.de/books?id=8Fq7BQAAQBAJ
  * [46] Seifert U 2012 Rep. Prog. Phys. 75 126001 URL https://doi.org/10.1088%2F0034-4885%2F75%2F12%2F126001
  * [47] Rosinberg M L, Munakata T and Tarjus G 2015 Phys. Rev. E 91(4) 042114 URL https://link.aps.org/doi/10.1103/PhysRevE.91.042114
  * [48] Rosinberg M L, Tarjus G and Munakata T 2017 Phys. Rev. E 95(2) 022123 URL https://link.aps.org/doi/10.1103/PhysRevE.95.022123
  * [49] Van Vu T and Hasegawa Y 2019 Phys. Rev. E 100(1) 012134 URL https://link.aps.org/doi/10.1103/PhysRevE.100.012134
  * [50] Loos S A M and Klapp S H L 2019 Sci. Rep. 9 2491 ISSN 2045-2322 URL https://doi.org/10.1038/s41598-019-39320-0
  * [51] Munakata T and Rosinberg M L 2014 Phys. Rev. Lett. 112(18) 180601 URL https://link.aps.org/doi/10.1103/PhysRevLett.112.180601
  * [52] Munakata T, Iwama S and Kimizuka M 2009 Phys. Rev. E 79(3) 031104 URL https://link.aps.org/doi/10.1103/PhysRevE.79.031104
  * [53] Adelman S A 1976 J. Chem. Phys. 64 124–130 URL https://doi.org/10.1063/1.431961
  * [54] Fox R F 1977 J. Math. Phys. 18 2331–2335 URL https://doi.org/10.1063/1.523242
  * [55] Hänggi P 1978 Z. Phys. B 31 407–416 ISSN 1431-584X URL https://doi.org/10.1007/BF01351552
  * [56] Sancho J M, Miguel M S, Katz S L and Gunton J D 1982 Phys. Rev. A 26(3) 1589–1609 URL https://link.aps.org/doi/10.1103/PhysRevA.26.1589
  * [57] Hernández-Machado A, Sancho J M, San Miguel M and Pesquera L 1983 Z. Phys. B 52 335–343 ISSN 1431-584X URL https://doi.org/10.1007/BF01307403
  * [58] Küchler U and Mensch B 1992 Stoch. Stoch. Rep. 40 23–42 URL https://doi.org/10.1080/17442509208833780
  * [59] Guillouzic S, L’Heureux I and Longtin A 1999 Phys. Rev. E 59(4) 3970–3982 URL https://link.aps.org/doi/10.1103/PhysRevE.59.3970
  * [60] Frank T D, Beek P J and Friedrich R 2003 Phys. Rev. E 68(2) 021912 URL https://link.aps.org/doi/10.1103/PhysRevE.68.021912
  * [61] McKetterick T J and Giuggioli L 2014 Phys. Rev. E 90(4) 042135 URL https://link.aps.org/doi/10.1103/PhysRevE.90.042135
  * [62] Giuggioli L, McKetterick T J, Kenkre V M and Chase M 2016 J. Phys. A: Math. Theor. 49 384002 URL https://doi.org/10.1088%2F1751-8113%2F49%2F38%2F384002
  * [63] Frank T D 2005 Phys. Rev. E 71(3) 031106 URL https://link.aps.org/doi/10.1103/PhysRevE.71.031106
  * [64] Loos S A M and Klapp S H L 2017 Phys. Rev. E 96(1) 012106 URL https://link.aps.org/doi/10.1103/PhysRevE.96.012106
  * [65] Loos S A M and Klapp S H L 2019 J. Stat. Phys. 177 95–118 ISSN 1572-9613 URL https://doi.org/10.1007/s10955-019-02359-4
  * [66] D’Anna G, Mayor P, Barrat A, Loreto V and Nori F 2003 Nature 424 909–912 ISSN 1476-4687 URL https://doi.org/10.1038/nature01867
  * [67] Murch K W, Vool U, Zhou D, Weber S J, Girvin S M and Siddiqi I 2012 Phys. Rev. Lett. 109(18) 183602 URL https://link.aps.org/doi/10.1103/PhysRevLett.109.183602
  * [68] Ferialdi L, Setter A, Toroš M, Timberlake C and Ulbricht H 2019 21 073019 URL https://doi.org/10.1088/1367-2630/ab2b69
  * [69] Kubo R 1966 Rep. Prog. Phys. 29 255–284 URL https://doi.org/10.1088%2F0034-4885%2F29%2F1%2F306
  * [70] Felderhof B U 1978 J. Phys. A: Math. Gen. 11 921–927 URL https://doi.org/10.1088%2F0305-4470%2F11%2F5%2F021
  * [71] Zwanzig R 2001 Nonequilibrium Statistical Mechanics (Oxford University Press)
  * [72] Kubo R, Toda M and Hashitsume N 2012 Statistical Physics II: Nonequilibrium Statistical Mechanics Springer Series in Solid-State Sciences (Springer Berlin Heidelberg) ISBN 9783642582448 URL https://books.google.de/books?id=cF3wCAAAQBAJ
  * [73] Callen H B 1985 Thermodynamics and an Introduction to Thermostatistics (2nd ed.) (New York: John Wiley & Sons)
  * [74] Evans D J, Searles D J and Williams S R 2009 2009 P07029 URL https://doi.org/10.1088/1742-5468/2009/07/p07029
  * [75] Evans D J, Searles D J and Williams S R 2016 Fundamentals of classical statistical thermodynamics: dissipation, relaxation, and fluctuation theorems (John Wiley & Sons)
  * [76] Maes C, Netočný K and Wynants B 2011 Phys. Rev. Lett. 107(1) 010601 URL https://link.aps.org/doi/10.1103/PhysRevLett.107.010601
  * [77] Risken H 1984 Fokker-Planck Equation (Berlin, Heidelberg: Springer Berlin Heidelberg)
  * [78] Van Kampen N 2011 Stochastic Processes in Physics and Chemistry North-Holland Personal Library (Amsterdam: Elsevier Science) ISBN 9780080475363 URL https://books.google.cz/books?id=N6II-6HlPxEC
  * [79] Jarzynski C 1997 Phys. Rev. Lett. 78(14) 2690–2693 URL https://link.aps.org/doi/10.1103/PhysRevLett.78.2690
  * [80] Crooks G E 1999 Phys. Rev. E 60(3) 2721–2726 URL https://link.aps.org/doi/10.1103/PhysRevE.60.2721
  * [81] Speck T and Seifert U 2007 Journal of Statistical Mechanics: Theory and Experiment 2007 L09002–L09002 URL https://doi.org/10.1088/1742-5468/2007/09/l09002
  * [82] Ortega R, Perez J, Nicklasson P and Sira-Ramirez H 2013 Passivity-based Control of Euler-Lagrange Systems: Mechanical, Electrical and Electromechanical Applications Communications and Control Engineering (Springer London) ISBN 9781447136033 URL https://books.google.cz/books?id=jJzeBwAAQBAJ
  * [83] Cugliandolo L F 2011 Journal of Physics A: Mathematical and Theoretical 44 483001 URL https://doi.org/10.1088/1751-8113/44/48/483001
  * [84] Holubec V, Steffenoni S, Falasco G and Kroy K 2020 Phys. Rev. Research 2(4) 043262 URL https://link.aps.org/doi/10.1103/PhysRevResearch.2.043262
  * [85] Goldstein H, Poole C and Safko J 2001 Classical mechanics (Addison Wesley San Francisco)
  * [86] Shinozuka M and Deodatis G 1991 Appl. Mech. Rev. 44 191–204 ISSN 0003-6900 URL https://doi.org/10.1115/1.3119501
  * [87] Pichot G 2016 Algorithms for stationary Gaussian random field generation Technical Report RT-0484 INRIA Paris URL https://hal.inria.fr/hal-01414707
  * [88] Graham I G, Kuo F Y, Nuyens D, Scheichl R and Sloan I H 2018 SIAM J. Numer. Anal. 56 1871–1895 URL https://doi.org/10.1137/17M1149730

##  Appendix A No feedback cooling with EQ feedback

Consider the heat flux $\dot{Q}_{\rm M}^{\rm EFB}$ from the proper bath into
the system in the steady state created by an equilibrium feedback. The force
applied by the proper bath on the system is
$-\gamma_{0}v(t)+\sqrt{2T_{0}\gamma_{0}}\xi(t)$ and thus

$\displaystyle\dot{Q}_{\rm M}^{\rm EFB}=$ $\displaystyle-\gamma_{0}\langle
v^{2}\rangle+\sqrt{2T_{0}\gamma_{0}}\,\langle\xi(t)v(t)\rangle$ (50)

with $\langle v^{2}\rangle=\lim_{t\to\infty}\langle
v(t)^{2}\rangle=\sigma^{\rm EFB}_{v}={T}/m$. To calculate
$\langle\xi(t)v(t)\rangle$, we use the formal solution

$v(t)=\frac{1}{m}\int_{0}^{t}dt^{\prime}\bigg{[}F[x(t^{\prime}),v(t^{\prime}),t^{\prime}]+F_{\rm
D}[x(t^{\prime}-\tau),v(t^{\prime}-\tau)]+\sqrt{2T_{0}\gamma_{0}}\xi(t^{\prime})+\eta_{\rm
FB}(t^{\prime})\bigg{]}.$ (51)

For a positive delay, $\tau>0$, the causality implies that all terms on the
R.H.S. except for $\xi(t^{\prime})$ are independent of the white noise
$\xi(t)$ at time $t$ (values of velocity and position at time $t^{\prime}\leq
t$ were not affected by the white noise yet, and it is also reasonable to
assume that the feedback noise $\eta_{\rm FB}(t^{\prime})$ can not be
constructed in such a way that it would depend on $\xi(t)$). In symbols we
obtain

$\langle\xi(t)v(t)\rangle=\frac{1}{m}\int_{0}^{t}dt^{\prime}\sqrt{2T_{0}\gamma_{0}}\delta(t-t^{\prime})=\frac{1}{m}\sqrt{\frac{T_{0}\gamma_{0}}{2}}.$
(52)

Using $\sigma^{\rm EQ}_{v}={T}_{0}/m$ we get the stationary flux $\dot{Q}_{\rm
M}$ induced by the EQ feedback in the form

$\dot{Q}_{\rm M}^{\rm EFB}=\gamma_{0}\left((\sigma_{v}^{\rm
EQ})^{2}-(\sigma_{v}^{\rm
EFB})^{2}\right)=\frac{2\gamma_{0}}{m}\left(T_{0}-T\right).$ (53)

As might have been anticipated from the beginning, the EQ feedback can cool
the proper bath ($\dot{Q}_{\rm M}^{\rm EFB}>0$) only if it leads to a smaller
velocity variance (effective temperature) than the proper bath. Let us now
show that this can happen only if the feedback force contains also a time-
local term in velocity.

Inserting the total noise $\eta(t)=\sqrt{2T_{0}\gamma_{0}}\xi(t)+\eta_{\rm
FB}(t)$ (see Tab. 1) into the FDR (7), we find

$\left<\eta(t)\eta(t^{\prime})\right>/{T}=2(\gamma_{0}+\delta)\delta(t-t^{\prime})+\dots=(2T_{0}\gamma_{0}+\epsilon)/{T}\,\delta(t-t^{\prime})+\dots.$
(54)

The first line corresponds to the time-local component of the friction kernel.
Specifically, the term proportional to $\gamma_{0}$ corresponds to the
background friction $\gamma_{0}v(t)$. And the term proportional to $\delta$
stems from the time-local component of the feedback force. The remaining terms
abbreviated by $\dots$ are determined by time non-local components of the
friction kernel. The second line corresponds to the actual noise correlations.
The term proportional to $2T_{0}\gamma_{0}$ is obtained from the background
noise and the term proportional to $\epsilon\delta(t-t^{\prime})$ originates
from the time-local ($t-t^{\prime}=0$) component of $\left<\eta_{\rm
FB}(t)\eta_{\rm FB}(t^{\prime})\right>$. The remaining terms are given by the
cross correlations $\left<\xi(t)\eta_{\rm FB}(t^{\prime})\right>$ and the
$t-t^{\prime}\neq 0$ component of $\left<\eta_{\rm FB}(t)\eta_{\rm
FB}(t^{\prime})\right>$.

Demanding that prefactors in front of the $\delta$-functions in Eq. (54)
equal, we find

$\frac{{T}_{0}}{{T}}=\frac{2(\gamma_{0}+\delta)-\epsilon/T}{2\gamma_{0}}.$
(55)

Since $\epsilon$ is the variance of the feedback noise $\eta_{\rm FB}(t)$, it
must be positive. The feedback can thus cool the system and the ambient bath
($T<T_{0}$ and $\dot{Q}_{\rm M}^{\rm EFB}>0$) only if it contains a strong
enough time-local component of the friction $-\delta v(t)$,
$\delta>\epsilon/2T$. Otherwise ${T}\geq{T}_{0}$ and $\dot{Q}_{\rm M}^{\rm
EFB}\leq 0$ and thus the cooling by equilibrium feedback is not possible.
Realizing a feedback force containing a term proportional to a non-delayed
velocity of the system seems technologically impossible and thus we conclude
that, unlike generic feedback, the equilibrium feedback cannot be used to cool
the system in practical setups [35, 37].

##  Appendix B Noise generation in practice

An arbitrary Gaussian noise with a given positive power spectrum, and thus
also any noise $\eta(t)$ fulfilling the conditions (7) and (9), can be
realized in practice using one of the standard procedures for generating a
Gaussian process with given autocorrelation function [86, 87, 88].

Concerning the position delay of Sec. 3 in the parameter regime (22), our
brute-force implementation of the spectral method using Eq. (35) in Ref. [86]
does not yield satisfactory results. However, we can recommended the method of
Refs. [87, 88] based on a discrete representation of the noise $\eta(t)$ by
$M_{t}{\mathbf{e}}$, where ${\mathbf{e}}$ is a column vector of independent
Gaussian white noises and the matrix $M_{t}$ is given by a square root of a
matrix describing the autocorrelation function (20).

A noise fulfilling the conditions (28) and (30) corresponding to the velocity
delay of Sec. 4 can be constructed analytically by introducing the ansatz

$\eta(t)=\sqrt{\alpha_{0}}\xi(t)+\sqrt{\alpha_{\tau}}\xi(t-\tau).$ (56)

Here, $\alpha_{0}=2{T}_{0}\gamma_{0}$, $\xi(t)$ is the zero-mean, unit-
variance, Gaussian white noise, i.e. $\left<\xi(t)\right>=0$,
$\left<\xi(t)\xi(t^{\prime})\right>=\delta(t-t^{\prime})$. For $\tau>0$, such
$\eta(t)$ obeys Eq. (28) if $(\alpha_{0}+\alpha_{\tau})/{T}=2\gamma_{0}$ and
$\sqrt{\alpha_{0}\alpha_{\tau}}/{T}=\gamma_{\tau}$. Solving these equations
with respect to $\alpha_{0}$ and $\alpha_{\tau}$, we obtain the two roots

$\alpha_{0}/{T}=\gamma_{0}\pm\sqrt{\gamma_{0}^{2}-\gamma_{\tau}^{2}},\quad\alpha_{\tau}/{T}=\gamma_{0}\mp\sqrt{\gamma_{0}^{2}-\gamma_{\tau}^{2}},$
(57)

which are equivalent due to the symmetry
$\alpha_{0}\leftrightarrow\alpha_{\tau}$ of the noise leading to the same
dynamics of $v(t)$. In agreement with our discussion of the power spectrum
(29), this mapping breaks down if $\gamma_{0}<|\gamma_{\tau}|$ when $\eta(t)$
becomes complex.

In standard Brownian dynamics simulations, the stochastic differential
equation $\dot{v}(t)=f(t,t-\tau)+\sigma\xi(t)$ is usually solved using the
Euler method leading to the update rule
$v(t+dt)=v(t)+f(t,t-\tau)dt+\sigma\sqrt{dt}L(t)$, where
$\sqrt{dt}L(t)=\int_{t}^{t+dt}dt^{\prime}\xi(t^{\prime})$ is a zero-mean
Gaussian random variable with variance $dt$. The equilibrium velocity delay
process (4) with the noise (56) has the form
$\dot{v}(t)=f(t,t-\tau)+\sqrt{\alpha_{0}}\xi(t)+\sqrt{\alpha_{\tau}}\xi(t-\tau)$.
And it can be simulated using the update rule
$v(t+dt)=v(t)+f(t,t-\tau)dt+\sqrt{\alpha_{0}dt}L(t)+\sqrt{\alpha_{\tau}dt}L(t-\tau)$.
To simulate the process, one thus has to keep track of the noise sequence
$L(t^{\prime})$ for $t^{\prime}\in[t-\tau,t]$.

##  Appendix C Parameter sets used in Figs. 3 and 4

The two parameter sets considered in Figs. 3 and 4 are chosen as follows. The
first one
$(\tilde{\gamma}_{0}\tau,\tilde{\gamma}_{\tau}\tau)\approx(0.28,0.28)$
minimizes the ratio

$\frac{t_{\rm R}}{t_{\rm R}^{\rm
EQ}}=\frac{\tilde{\gamma}_{0}\tau}{\Re\left[\tilde{\gamma}_{0}\tau-W\left(-e^{\tilde{\gamma}_{0}\tau}\tilde{\gamma}_{\tau}\tau\right)\right]}$
(58)

of the relaxation time $t_{\rm R}^{\rm EQ}\equiv 1/\tilde{\gamma}_{0}$ for the
EQ process and $t_{\rm R}$ (31) for the EFB and NEFB for $U=0$. The second one
$(\tilde{\gamma}_{0}\tau,\tilde{\gamma}_{\tau}\tau)\approx(0.42,0.24)$
minimizes the measure

$\frac{t_{\rm R}}{t_{\rm R}^{\rm EQ}}\left(\frac{\sigma_{v}^{\rm
EFB}}{\sigma^{\rm NEFB}_{v}}\right)^{2}$ (59)

of the trade-off between the relaxation time ratio (58) and the ratio of
stationary velocity variance $(\sigma_{v})^{2}=\left<v^{2}\right>$ for EFB and
NEFB. For NEFB we take the variance for $U=0$, when it can be calculated
analytically as

$\left(\sigma_{v}^{\rm
NEFB}\right)^{2}\equiv\sigma_{v}^{2}=\frac{T_{0}}{m}\frac{\tilde{\gamma}_{0}[\Omega+\tilde{\gamma}_{\tau}\sinh(\Omega\tau)]}{\Omega\left[\tilde{\gamma}_{0}+\tilde{\gamma}_{\tau}\cosh(\Omega\tau)\right]},$
(60)

where $\Omega=\sqrt{\tilde{\gamma}_{0}^{2}-\tilde{\gamma}_{\tau}^{2}}$. For
EFB, the variance $\left(\sigma_{v}^{\rm EFB}\right)^{2}={T}/m$ follows for
any $U$ from equipartition.

In addition, we also tested the parameter set
$(\tilde{\gamma}_{0}\tau,\tilde{\gamma}_{\tau}\tau)\approx(0.34,0.034)$ which
minimizes the ratio $\sigma_{v}^{\rm EFB}/\sigma^{\rm NEFB}_{v}$ but due to
the small magnitude of the feedback force the numerical results for NEFB, EFB,
and EQ are hardly distinguishable and we decided not to show them.

##  Appendix D Derivation of Eq. (60)

Consider the simple velocity process

$\dot{v}(t)=-\tilde{\gamma}_{0}v(t)-\tilde{\gamma}_{\tau}v(t-\tau)+\sqrt{\alpha}\xi(t).$
(61)

In the steady state, $\left<v\right>=0$ and the variance
$\sigma^{2}_{v}=\left<v^{2}\right>$ can be evaluated as follows.

The general solution to Eq. (61) reads [42]

$v(t)=\lambda(t)v_{0}-\tilde{\gamma}_{\tau}\int_{-\tau}^{0}dt^{\prime}\,\lambda(t-t^{\prime}-\tau)v(t^{\prime})+\sqrt{\alpha}\int_{0}^{t}dt^{\prime}\,\lambda(t-t^{\prime})\xi(t^{\prime}),$
(62)

where $v_{0}=v(0)$. The Green’s function $\lambda(t)$ solves Eq. (61) with
vanishing noise term ($\alpha=0$), i.e.

$\dot{\lambda}(t)=-\tilde{\gamma}_{0}\lambda(t)-\tilde{\gamma}_{\tau}\lambda(t-\tau),$
(63)

and the initial condition $\lambda(t<0)=0$ and $\lambda(0)=1$. The most
straightforward way for finding $\lambda(t)$ is to employ a Laplace
transformation in time. The result is [61]

$\lambda(t)=\sum_{l=0}^{\infty}\frac{(-\tilde{\gamma}_{\tau})^{l}}{l!}\left(t-l\tau\right)^{l}{\rm
e}^{-\tilde{\gamma}_{0}\left(t-l\tau\right)}\Theta\left(t-l\tau\right).$ (64)

In the stable regime, $t_{\rm R}>0$ (cf Eq. (31)), where the Green’s function
eventually decays to zero, $\lim_{t\to\infty}\lambda(t)=0$, the general
solution (62) to Eq. (61) can be used for calculation of the time-correlation
function, $C(t)=\lim_{t_{0}\to\infty}\left<v(t_{0}+t)v(t_{0})\right>$. The
stationary variance $\sigma_{v}^{2}=C(0)$, we are actually interested in,
comes as a by-product of this calculation.

We find that for $t>0$ (see also Ref. [42])

$C(t)=\alpha\lim_{t_{0}\to\infty}\int_{0}^{t_{0}}dt^{\prime}\,\lambda(t+t_{0}-t^{\prime})\lambda(t_{0}-t^{\prime}).$
(65)

This expression can already be used for plotting the time correlation
function, however, it is not very suitable for inferring its properties.
Following the approach in Refs. [60, 42], we take the time derivative of Eq.
(65) and use Eq. (63) for the Green’s function to obtain the dynamical
equation

$\dot{C}(t)=-\tilde{\gamma}_{\tau}C(t-\tau)-\tilde{\gamma}_{0}C(t)$ (66)

valid for $t>0$ due to the nonanalyticity of $\lambda(t)$ at $t=0$. The
solution to this equation is given by Eq. (62):

$C(t)=\lambda(t)C_{0}-\tilde{\gamma}_{\tau}\int_{-\tau}^{0}dt^{\prime}\,\lambda(t-t^{\prime}-\tau)C(t^{\prime})$
(67)

and thus the decay time of the time-correlation function is given by the decay
time (31) of the Green’s function $\lambda(t)$. To evaluate the above
expression, we need to find the stationary variance $\alpha_{0}\equiv C(0)$
and the delayed initial condition $C(t)$ for $t\in(-\tau,0)$. This can be done
as follows. Employing the symmetry $C(t)=C(-t)$ of the stationary time-
correlation function, we rewrite Eq. (66) as

$\dot{C}(t)=-\tilde{\gamma}_{\tau}C(\tau-t)-\tilde{\gamma}_{0}C(t).$ (68)

For $t\in(0,\tau)$, we can differentiate this equation once again. The result
is

$\ddot{C}(t)=\Omega^{2}C(t),$ (69)

where we used Eq. (66) and defined the (possibly imaginary) frequency
$\Omega=\sqrt{\tilde{\gamma}_{0}^{2}-\tilde{\gamma}_{\tau}^{2}}$. From Eq.
(69), we find that for $t\in[-\tau,\tau]$

$C(t)=C_{0}\cosh\left(\Omega
t\right)+\dot{C}_{0}\Omega^{-1}\sinh\left(\Omega|t|\right).$ (70)

Here, $\dot{C}_{0}=\lim_{t\to 0+}\dot{C}(t)$ denotes the time-derivative of
the time-correlation function, that is discontinuous at $t=0$ [60], for $t>0$
infinitesimally close to 0. From Eq. (65), we find that

$\dot{C}_{0}=\alpha\lim_{t_{0}\to\infty}\int_{0}^{t_{0}}dt^{\prime}\dot{\lambda}(t_{0}-t^{\prime})\lambda(t_{0}-t^{\prime})=-0.5\alpha\lim_{t_{0}\to\infty}\int_{0}^{t_{0}}dt^{\prime}\frac{d}{dt^{\prime}}\lambda^{2}(t_{0}-t^{\prime})\\\
=-0.5\alpha\lim_{t_{0}\to\infty}\left[\lambda^{2}(0)-\lambda^{2}(t_{0})\right]=-0.5\alpha.$
(71)

In order to evaluate $C_{0}$, we note that $\dot{C}_{0}$ also follows from Eq.
(68) with $t=0$, yielding
$\dot{C}_{0}=-0.5\alpha=-\tilde{\gamma}_{\tau}C(\tau)-\tilde{\gamma}_{0}C_{0}$.
Using
$C(\tau)=C_{0}\cosh\left(\Omega\tau\right)-0.5\alpha\Omega^{-1}\sinh\left(\Omega\tau\right)$
given by Eq. (70) for $t>0$, we finally get the desired result

$C_{0}=\sigma_{v}^{2}=\frac{\alpha}{2}\frac{\Omega+\tilde{\gamma}_{\tau}\sinh(\Omega\tau)}{\Omega\left[\tilde{\gamma}_{0}+\tilde{\gamma}_{\tau}\cosh(\Omega\tau)\right]}.$
(72)

For $\tilde{\gamma}_{0}=0$, we obtain $\Omega=i\sqrt{\tilde{\gamma}_{\tau}}$.
Using the identities $\sinh(ix)=i\sin x$ and $\cosh(ix)=\cos x$, the formula
(72) can be written as

$\sigma_{v}^{2}=\frac{\alpha}{2\tilde{\gamma}_{\tau}}\frac{1+\sin(\tilde{\gamma}_{\tau}\tau)}{\cos(\tilde{\gamma}_{\tau}\tau)},$
(73)

which is the result derived in Refs. [60, 42].