#  Derivation of the non-equilibrium generalized Langevin equation from a
generic time-dependent Hamiltonian

Roland R. Netz  Fachbereich Physik, Freie Universität Berlin, 144195 Berlin,
Germany  Centre for Condensed Matter Theory, Department of Physics, Indian
Institute of Science, Bangalore 560012, India

###### Abstract

It has been become standard practice to describe steady-state non-equilibrium
phenomena by Langevin equations with colored noise and time-dependent friction
kernels that do not obey the fluctuation-dissipation theorem, but since these
Langevin equations are typically not derived from first-principle Hamiltonian
dynamics it is not clear whether they correspond to physically realizable
scenarios. By exact Mori projection in phase space we derive the non-
equilibrium generalized Langevin equation (GLE) from a generic many-body
Hamiltonian with a time-dependent force h(t) acting on an arbitrary phase-
space dependent observable $A$. The GLE is obtained in explicit form to all
orders in $h(t)$. For non-equilibrium observables that correspond to a
Gaussian process, the resultant GLE has the same form as the equilibrium Mori
GLE, in particular the memory kernel is proportional to the total force
autocorrelation function. This means that the extraction and simulation
methods developed for equilibrium GLEs can be used also for non-equilibrium
Gaussian variables. This is a non-trivial and very useful result, as many
observables that characterize non-equilibrium systems display Gaussian
statistics. For non-Gaussian non-equilibrium variables correction terms appear
in the GLE and in the relation between the complementary force
autocorrelations and the memory kernels, which are explicitly given in term of
cubic correlation functions of $A$. Interpreting the time-dependent force h(t)
as a stochastic process, we derive non-equilibrium corrections to the
fluctuation-dissipation theorem and methods to extract all GLE parameters from
experimental or simulation data, thus making our non-equilibrium GLE a
practical tool to study and model general non-equilibrium systems.

##  I Introduction

The statistical mechanics foundation of non-equilibrium phenomena has occupied
physicists for many decades [1, 2, 3, 4, 5, 6, 7, 8]. More recently, new
experimental techniques, such as single-molecule and optical methods, applied
to non-equilibrium biological systems have accented the need for theories that
are able to deal with non-equilibrium experiments and data [9, 10, 11, 12, 13,
14, 15]. At the same time, novel theoretical approaches were developed and
applied to non-equilibrium driven lattice models [16, 17, 18], interacting
non-equilibrium particle systems [19, 20, 21], non-equilibrium barrier-
crossing phenomena [22, 23] and used to derive non-equilibrium work and
entropy relations [24, 25, 26, 27, 28], generalized fluctuation-dissipation
relations [29, 30, 31, 32, 33] and non-equilibrium entropy-production extremal
principles [34].

The generalized Langevin equation (GLE) has played a key role in the
development of methods to deal with the dynamics of complex systems, as it is
an exact equation of motion for an observable derived by projection from the
many-body Hamiltonian, the GLE thus constitutes a method for exact coarse-
graining of a Hamiltonian system [35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45,
46, 47]. The GLE was applied to protein folding [48, 49, 50, 51], barrier
crossing dynamics [52, 53, 54, 22, 55], motion of living cells [56],
spectroscopy [57, 58, 59], dynamical neworks [60] and data prediction [61].
While the standard GLE formulations describe the motion of an observable in
phase space and thus allow to quantify the approach of a non-equilibrium state
to equilibrium, they do not apply to driven non-equilibrium system as
described by a time-dependent Hamiltonian. Many works dealt with
generalizations of the projection framework to time-dependent and transient
scenarios [62, 63, 64, 65, 66, 67, 68, 69, 70, 71]. Non of these works dealt
with the non-equilibrium Hamiltonian system considered in this paper and
derived the non-equilibrium GLE in closed form. Since most of the
aforementioned theories concerned with non-equilibrium phenomena break the
fluctuation-dissipation theorem more or less by hand, it is instructive to
derive non-equilibrium equations of motion for observables from time-dependent
Hamiltonians. This enables to check which non-equilibrium effective equations
correspond to an underlying Hamiltonian non-equilibrium dynamics and which do
not. This is the vantage point of this paper.

In the presence of a time-dependent force $h(t)$, a Hamiltonian system is
generally out of equilibrium since the force performs work on the system. In
fact, even for constant force $h(t)=h_{0}\neq 0$ such a system is out of
equilibrium if the observable $A$ is unconfined and thus driven into a steady-
state motion by the constant force, as will be explained below. For an
equilibrium system, i.e. for $h(t)=0$, the Mori GLE for the observable $A(t)$
reads [37]

$\displaystyle\ddot{A}(t)=-K(A(t)-\langle A\rangle)-\int_{t_{0}}^{t}{\rm
d}s\,\Gamma(t-s)\dot{A}(s)+F(t),$ (1)

where the stiffness of the effective harmonic potential is denoted as $K$, the
time-dependent friction memory kernel as $\Gamma(t)$ and the complementary
force as $F(t)$. The time at which the projection is done is denoted as
$t_{0}$ and the friction kernel $\Gamma(t)$ is related to the complementary
force autocorrelation function via [37]

$\displaystyle\Gamma(t-s)=\frac{\langle
F(s)F(t)\rangle}{\langle\dot{A}^{2}(t)\rangle}$ (2)

(all averages are phase-space averages, as will be detailed below). Note that
$F(t)$ is often denoted and treated as a random force, this is an
approximation since $F(t)$ is in fact a phase-space dependent deterministic
function and fulfills well-defined initial conditions at $t=t_{0}$, Eq. (1) is
thus deterministic and fully time reversible. The Mori GLE in Eq. (1) is also
exact even for non-linear systems, unless the complementary force $F(t)$ is
approximated as Gaussian.

In this paper we derive the non-equilibrium generalized Langevin equation
(GLE) for an arbitrary phase-space dependent observable $A$, governed by a
general many-body Hamiltonian that includes a time-dependent external force
$h(t)$ acting on $A$, the derivation is exact to all orders of $h(t)$. This
specific time-dependent Hamiltonian is of high relevance, as it forms the
starting point for the text-book derivation of the standard fluctuation-
dissipation theorem (FDT) [8], which is one of the corner stones of
statistical mechanics. Thus the standard FDT and our non-equilibrium GLE are
intimately connected since they stem from the same Hamiltonian.

A key question we address in this paper is whether in the presence of an
external time-dependent force $h(t)$ acting on $A$, the GLE in Eq. (1) and the
relation between the friction kernel and the complementary force
autocorrelation in Eq. (2) still hold. Indeed, one main result of this paper
is that the simple forms of Eqs. (1) and (2) indeed remain valid if the
observable $A$ corresponds to a Gaussian process and if $F(t)$ is replaced by
the combination of $F(t)$ and $h(t)$. This is a non-trivial and relevant
finding, since many biological non-equilibrium processes, such as the motion
of cells, are Gaussian to high accuracy [56]. Conversely, for a non-Gaussian
observable $A$, correction terms appear in Eqs. (1) and (2) that are
explicitly given in terms of three-point (and higher-order) correlation
functions of $A$. These explicit results allow for explicit prediction of the
non-equilibrium correction terms based on experimental or simulation time-
series data. Our derivation of the GLE is non-perturbative and thus exact to
all orders in the non-equilibrium force $h(t)$. For the special case of a
stochastic non-equilibrium force $h(t)$ that is defined by its second moment,
we derive a generalized non-equilibrium FDT which in the limit
$h(t)\rightarrow 0$ simplifies to the standard FDT. Finally, we give explicit
formulas for extracting the parameters of our non-equilibrium GLE from
simulation or experimental time series data, opening the route to the accurate
and data-based modeling of non-equilibrium systems.

Sect. II contains the full derivation of the non-equilibrium GLE, this section
can be skipped by a reader not interested in technical details. In Sect. III
the non-equilibrium GLE is discussed and the role of non-Gaussian observable
fluctuations is explained. In Sect. IV the non-equilibrium force $h(t)$ is
treated as a stochastic variable, which restores time-homogeneity and
simplifies the analysis of the GLE. Here the non-equilibrium FDT is derived.
Sect. V presents a short discussion and an outlook.

##  II Derivation of the non-equilibrium generalized Langevin equation

###  II.1  Definition of the time-dependent Hamiltonian and solution of the
Liouville equation

We consider a time-dependent Hamiltonian for a system of $N$ interacting
particles or atoms in three-dimensional space with a time-dependent force
$h(t)$ that couples to a generic phase-space dependent observable
$A_{S}(\omega)$,

$\displaystyle H(\omega,t)$ $\displaystyle=H_{0}(\omega)-h(t)A_{S}(\omega),$
(3)

the subscript distinguishes this Schrödinger-like, i.e. time-independent,
observable from the time-dependent Heisenberg observable that will be
introduced shortly. Although not really needed for our derivation, the time-
independent Hamiltonian $H_{0}(\omega)$ can be split into kinetic and
potential contributions according to

$\displaystyle H_{0}(\omega)$
$\displaystyle=\sum_{j=1}^{3N}\frac{P_{j}^{2}}{2m_{j}}+V(\mathbf{R})$ (4)

with coordinate-dependent masses $m_{j}$ and where the potential
$V(\mathbf{R})$ contains all interactions between the particles and includes
possible external potentials. A point in 6$N$-dimensional phase space is
denoted by $\omega=(\mathbf{R},\mathbf{P})$, which is a $6N$-dimensional
vector containing the Cartesian particle positions $\mathbf{R}$ and the
conjugate momenta $\mathbf{P}$ and fully specifies the microstate of the
system.

Using the time-dependent Liouville operator

$\displaystyle{\cal L}(\omega,t)$
$\displaystyle=\sum_{j=1}^{3N}\left(\frac{\partial H(\omega,t)}{\partial
P_{j}}\frac{\partial}{\partial R_{j}}-\frac{\partial H(\omega,t)}{\partial
R_{j}}\frac{\partial}{\partial P_{j}}\right),$ (5)

the $6N$-dimensional Hamilton equation of motion can be compactly written as
$\dot{\omega}(t)={\cal L}(\omega,t)\omega(t)$, where $\omega(t)$ is the phase-
space location of the system at time $t$ and $\dot{\omega}(t)={\rm
d}\omega(t)/{\rm d}t$ is the corresponding phase-space velocity. Instead of
following microstate trajectories in phase space, which is the Lagrangian
description of the system dynamics, it is much more convenient to switch to
the Eulerian description and consider the time dependent probability density
distribution as a function of the time-independent phase-space position,
$\rho(\omega,t)$, which obeys the Liouville equation

$\displaystyle\dot{\rho}(\omega,t)$ $\displaystyle=-{\cal
L}(\omega,t)\rho(\omega,t).$ (6)

In all what follows we suppress the dependence of the Liouville operator on
phase space. We observe that a recursive solution of Eq. (6) can be written as

$\displaystyle{\rho}(\omega,t)$
$\displaystyle=\rho(\omega,t_{0})-\int_{t_{0}}^{t}{\rm d}t_{1}{\cal
L}(t_{1})\rho(\omega,t_{1}).$ (7)

By iteration the following exact solution is obtained

$\displaystyle{\rho}(\omega,t)$
$\displaystyle=\exp_{S}\left(-\int_{t_{0}}^{t}{\rm d}s{\cal
L}(s)\right)\rho(\omega,t_{0}),$ (8)

which depends on the initial density distribution at time $t_{0}$ and where
the time-ordered operator exponential in the Schrödinger picture has been
introduced as [72, 73]

$\displaystyle\exp_{S}\left(-\int_{t_{0}}^{t}{\rm d}s{\cal L}(s)\right)$
$\displaystyle\equiv 1+\sum_{n=1}^{\infty}(-1)^{n}\int_{t_{0}}^{t}{\rm
d}t_{1}{\cal L}(t_{1})\int_{t_{0}}^{t_{1}}{\rm d}t_{2}{\cal
L}(t_{2})\int_{t_{0}}^{t_{2}}{\rm d}t_{3}{\cal
L}(t_{3})\cdots\int_{t_{0}}^{t_{n-1}}{\rm d}t_{n}{\cal L}(t_{n}).$ (9)

For a time-independent Liouville operator ${\cal L}(t)={\cal L}_{0}$, all time
integrals can be done and one obtains the solution in the terms of the
standard operator exponential

$\displaystyle{\rho}(\omega,t)$ $\displaystyle=\exp\left(-(t-t_{0}){\cal
L}_{0}\right)\rho(\omega,t_{0}),$ (10)

where the exponential of an operator is defined by its ordinary series
expansion.

###  II.2 From Schrödinger to Heisenberg observables

A system observable can be generally written as a Schrödinger-type phase-space
function $A_{S}(\omega)$, it can for example represent the position of one
particle, the center-of-mass position of a group of particles or of a
molecule, the reaction coordinate describing a chemical reaction or the
folding of a protein. To simplify the notation, we consider a scalar
observable but note that the formalism can be straightforwardly extended also
to multi-dimensional observables. Using the probability density
${\rho}(\omega,t)$, the time-dependent expectation value (or mean) of the
observable $A_{S}(\omega)$ can be written as

$\displaystyle a(t)$ $\displaystyle\equiv\int{\rm
d}\omega\,A_{S}(\omega)\rho(\omega,t)$ $\displaystyle=\int{\rm
d}\omega\,A_{S}(\omega)\exp_{S}\left(-\int_{t_{0}}^{t}{\rm d}s{\cal
L}(s)\right)\rho(\omega,t_{0}).$ (11)

Since the Liouville operator is anti-self adjoint, it follows that [8]

$\displaystyle a(t)$ $\displaystyle=\int{\rm
d}\omega\,\rho(\omega,t_{0})A(\omega,t),$ (12)

where we have defined the Heisenberg observable as

$\displaystyle A(\omega,t)$
$\displaystyle\equiv\exp_{H}\left(\int_{t_{0}}^{t}{\rm d}s{\cal
L}(s)\right)A_{S}(\omega)$ (13)

using the time-ordered operator exponential in the Heisenberg picture (or
Heisenberg propagator)

$\displaystyle\exp_{H}\left(\int_{t_{0}}^{t}{\rm d}s{\cal L}(s)\right)$
$\displaystyle\equiv 1+\sum_{n=1}^{\infty}\int_{t_{0}}^{t}{\rm
d}t_{1}\int_{t_{0}}^{t_{1}}{\rm d}t_{2}\cdots\int_{t_{0}}^{t_{n-1}}{\rm
d}t_{n}{\cal L}(t_{n})\cdots{\cal L}(t_{2}){\cal L}(t_{1}).$ (14)

Obviously, as follows from Eqs. (13) and (14), the Heisenberg observable
satisfies the equation of motion

$\displaystyle\dot{A}(\omega,t)$ $\displaystyle=\frac{{\rm d}A(\omega,t)}{{\rm
d}t}=\exp_{H}\left(\int_{t_{0}}^{t}{\rm d}s{\cal L}(s)\right){\cal
L}(t)A_{S}(\omega)$ (15)

with the initial condition $A(\omega,t_{0})=A_{S}(\omega)$. As derived in
Appendix A, the Heisenberg observable also satisfies the initial differential
boundary condition

$\displaystyle\frac{{\rm d}A(\omega,t)}{{\rm d}t_{0}}$ $\displaystyle=-{\cal
L}(t_{0})\exp_{H}\left(\int_{t_{0}}^{t}{\rm d}s{\cal
L}(s)\right)A_{S}(\omega),$ (16)

which will be later needed to derive operator expansions.

To understand the meaning of a Heisenberg observable, we for the moment
consider the initial density distribution
$\rho(\omega,t_{0})=\delta(\omega-\omega_{0})$, which describes a system that
at time $t_{0}$ is in the microstate $\omega_{0}$. Inserting this into Eq.
(12), we obtain $a(t)=A(\omega_{0},t)$. In other words, $A(\omega_{0},t)$
describes the time-dependent mean of an observable for a system that at time
$t=t_{0}$ was in the microstate $\omega_{0}$, i.e., it describes the temporal
evolution of the conditional mean of the observable $A_{S}(\omega)$. It
transpires that if we derive an equation of motion for $A(\omega,t)$, we have
an equation for how this conditional mean changes in time. This is the central
idea of projection and of GLEs [35, 36, 37].

Taking another time derivative of Eq. (15), we obtain for the acceleration of
the observable

$\displaystyle\ddot{A}(\omega,t)$
$\displaystyle=\exp_{H}\left(\int_{t_{0}}^{t}{\rm d}s{\cal
L}(t)\right)\left({\cal L}^{2}(t)+\dot{\cal L}(t)\right)A_{S}(\omega).$ (17)

Up to now the discussion applied to a general time-dependent Hamiltonian; for
the specific Hamiltonian Eq. (3), where a time dependent force $h(t)$
multiplies the observable $A_{S}(\omega)$, the Liouville operator splits into
two parts

$\displaystyle{\cal L}(t)$ $\displaystyle={\cal L}_{0}-h(t)\Delta{\cal L},$
(18)

with the unperturbed Liouville operator given by

$\displaystyle{\cal L}_{0}$ $\displaystyle=\sum_{j=1}^{3N}\left(\frac{\partial
H_{0}(\omega)}{\partial P_{j}}\frac{\partial}{\partial R_{j}}-\frac{\partial
H_{0}(\omega)}{\partial R_{j}}\frac{\partial}{\partial P_{j}}\right),$ (19)

and the perturbation Liouville operator given by

$\displaystyle\Delta{\cal L}$
$\displaystyle=\sum_{j=1}^{3N}\left(\frac{\partial A_{S}(\omega)}{\partial
P_{j}}\frac{\partial}{\partial R_{j}}-\frac{\partial A_{S}(\omega)}{\partial
R_{j}}\frac{\partial}{\partial P_{j}}\right).$ (20)

These operators have the important properties ${\cal L}_{0}H_{0}(\omega)=0$,
$\Delta{\cal L}A_{S}(\omega)=0$, and $\Delta{\cal L}H_{0}(\omega)=-{\cal
L}_{0}A_{S}(\omega)$, from which we derive, using Eqs. (15) and (17), the
simplified expressions for the observable velocity and acceleration

$\displaystyle\dot{A}(\omega,t)$
$\displaystyle=\exp_{H}\left(\int_{t_{0}}^{t}{\rm d}s{\cal L}(s)\right){\cal
L}_{0}A_{S}(\omega),$ (21) $\displaystyle\ddot{A}(\omega,t)$
$\displaystyle=\exp_{H}\left(\int_{t_{0}}^{t}{\rm d}s{\cal L}(s)\right){\cal
L}(t){\cal L}_{0}A_{S}(\omega).$ (22)

The fact that the velocity Eq. (21) exhibits no time dependence to the right
of the operator exponential is crucial, as it will later on allow us to use
time-independent projection for the derivation of the non-equilibrium GLE.

###  II.3  Projection

Here we follow standard procedures [35, 36, 37, 8]. We introduce a time-
independent projection operator ${\cal P}$ that acts on a phase space function
and its complementary operator ${\cal Q}$ via the relation $1={\cal Q}+{\cal
P}$. Inserting this unit operator into the acceleration Eq. (22), we obtain

$\displaystyle\ddot{A}(\omega,t)$
$\displaystyle=\exp_{H}\left(\int_{t_{0}}^{t}{\rm d}s{\cal L}(s)\right)({\cal
P}+{\cal Q}){\cal L}(t){\cal L}_{0}A_{S}(\omega)$
$\displaystyle=\exp_{H}\left(\int_{t_{0}}^{t}{\rm d}s{\cal L}(s)\right){\cal
P}{\cal L}(t){\cal L}_{0}A_{S}(\omega)$
$\displaystyle+\exp_{H}\left(\int_{t_{0}}^{t}{\rm d}s{\cal L}(s)\right){\cal
Q}{\cal L}(t){\cal L}_{0}A_{S}(\omega),$ (23)

where we used that the Heisenberg propagator is a linear operator. The
projection is performed at time $t_{0}$ at which the time propagation starts
(the relevance of this will become clear later). By inserting the time-
dependent Dyson operator expansion [72, 73, 35, 36, 37] for the Heisenberg
propagator (see Appendix B for a derivation)

$\displaystyle\exp_{H}\left(\int_{t_{0}}^{t}{\rm d}s{\cal
L}(s)\right)=\exp_{H}\left({\cal Q}\int_{t_{0}}^{t}{\rm d}s{\cal
L}(s)\right)+$ (24) $\displaystyle\int_{t_{0}}^{t}{\rm
d}s\exp_{H}\left(\int_{t_{0}}^{s}{\rm d}s^{\prime}{\cal
L}(s^{\prime})\right){\cal P}{\cal L}(s)\exp_{H}\left({\cal Q}\int_{s}^{t}{\rm
d}s^{\prime}{\cal L}(s^{\prime})\right)$

into the second term on the right hand side in Eq. (II.3), we obtain the GLE
in general form

$\displaystyle\ddot{A}(\omega,t)$
$\displaystyle=\exp_{H}\left(\int_{t_{0}}^{t}{\rm d}s{\cal L}(s)\right){\cal
P}{\cal L}(t){\cal L}_{0}A_{S}(\omega)+F(\omega,t_{0},t)$
$\displaystyle+\int_{t_{0}}^{t}{\rm d}s\exp_{H}\left(\int_{t_{0}}^{s}{\rm
d}s^{\prime}{\cal L}(s^{\prime})\right){\cal P}{\cal L}(s)F(\omega,s,t),$ (25)

where the complementary force is defined as

$\displaystyle F(\omega,t_{0},t)$ $\displaystyle\equiv\exp_{H}\left({\cal
Q}\int_{t_{0}}^{t}{\rm d}s^{\prime}{\cal L}(s^{\prime})\right){\cal Q}{\cal
L}(t){\cal L}_{0}A_{S}(\omega).$ (26)

The first term on the right-hand side in Eq. (II.3) will turn out to represent
the conservative force from a potential, the third term represents friction
and non-Markovian effects and the force $F(\omega,t_{0},t)$ represents all
effects that are not included in the other two terms. $F(\omega,t_{0},t)$ is a
function of phase space and evolves in the complementary space, i.e. it
satisfies ${\cal P}F(\omega,t_{0},t)=0$ (as will be explained further below).
While we suppress the $t_{0}$ dependence of the observable $A(\omega,t)$ and
its derivatives, which can cause no confusion since this argument is invariant
throughout most of the calculation, the complementary force
$F(\omega,t_{0},t)$ needs both time arguments since both arguments are varied
in Eq. (II.3).

Clearly, the explicit form of Eq. (II.3) depends on the specific projection
operator $\cal P$. Here we choose the Mori projection, because it is most
straightforward to implement and our main result concerning the effect of non-
Gaussian observables on the structure of the non-equilibrium GLE is accurately
and transparently produced by Mori projection. We note in passing that the
Mori GLE is exact even for non-Gaussian observables, unless one approximates
the complementary force distribution. The Mori projection applied on a general
Heisenberg observable $B(\omega,t)$ using the Schrödinger observable
$A_{S}(\omega)$ as a projection function is given by [37]

$\displaystyle{\cal P}B(\omega,t)=\langle B(\omega,t)\rangle+\frac{\langle
B(\omega,t){\cal L}_{0} A_{S}(\omega)\rangle}{\langle({\cal
L}_{0}A_{S}(\omega))^{2}\rangle}{\cal L}_{0}A_{S}(\omega)$
$\displaystyle+\frac{\langle B(\omega,t)(A_{S}(\omega)-\langle
A_{S}\rangle)\rangle}{\langle(A_{S}(\omega)-\langle
A_{S}\rangle)^{2}\rangle}(A_{S}(\omega)-\langle A_{S}\rangle).$ (27)

Here we have defined the expectation value of an arbitrary phase-space
function $X(\omega)$ with respect to a time-independent projection
distribution $\rho_{\rm p}(\omega)$ as

$\displaystyle\langle X(\omega)\rangle=\int{\rm d}\omega X(\omega)\rho_{\rm
p}(\omega),$ (28)

which we here take to be the equilibrium canonical distribution of the time-
independent Hamiltonian

$\displaystyle\rho_{\rm p}(\omega)=e^{-\beta H_{0}(\omega)+\beta
h_{p}A_{S}(\omega)}/Z,$ (29)

where $Z$ is the partition function. The factor $\beta$ has units of inverse
energy and can be thought of as the inverse thermal energy characterizing the
projection distribution. Note that for generality we added a linear force
$h_{p}$ in the projection Hamiltonian. Time-dependent projection has been used
to derive generic non-equilibrium GLEs [62, 63, 64, 65, 66, 67, 68, 69, 70,
71], but is not needed here because of the specific form of our time-dependent
Hamiltonian. The time-independent Mori projection in Eq. (II.3) projects onto
a constant, the Schrödinger observable $A(\omega,t_{0})=A_{S}(\omega)$ and its
time derivative $\dot{A}(\omega,t_{0})={\cal L}_{0}A_{S}(\omega)$. Thus the
projection in Eq. (II.3) maps any observable $B(\omega,t)$ onto the subspace
of all functions linear in the observables 1, $A_{S}(\omega)$ and ${\cal
L}_{0}A_{S}(\omega)$, meaning that ${\cal P}1=1$, ${\cal
P}A_{S}(\omega)=A_{S}(\omega)$ and ${\cal P}{\cal L}_{0}A_{S}(\omega)={\cal
L}_{0}A_{S}(\omega)$. From this follows immediately that ${\cal Q}1={\cal
Q}A_{S}(\omega)={\cal Q}{\cal L}_{0}A_{S}(\omega)=0$, which are important
properties. In particular, it follows that several expectation values
involving the complementary force vanish, namely $\langle
F(\omega,t_{0},t)\rangle=\langle F(\omega,t_{0},t)A_{S}(\omega)\rangle=\langle
F(\omega,t_{0},t){\cal L}_{0}A_{S}(\omega,)\rangle=0$, which will be used to
extract GLE parameters from non-equlibrium time-series data,

The Mori projection is linear, i.e., for two arbitrary observables
$B(\omega,t)$ and $C(\omega,t^{\prime})$ it satisfies ${\cal P}
(c_{1}B(\omega,t)+c_{2}C(\omega,t^{\prime}))=c_{1}{\cal
P}B(\omega,t)+c_{2}{\cal P}C(\omega,t^{\prime})$, it is idempotent, i.e.,
${\cal P}^{2}={\cal P}$, and it is self-adjoint, i.e. it satisfies the
relation

$\displaystyle\langle C(\omega,t){\cal P}B(\omega,t^{\prime})\rangle=\langle
B(\omega,t^{\prime}){\cal P}C(\omega,t)\rangle.$ (30)

From these properties it follows that the complementary projection operator
${\cal Q}=1-{\cal P}$ is also linear, idempotent and self-adjoint. Thus,
${\cal P}$ and ${\cal Q}$ are orthogonal to each other, i.e. ${\cal P}{\cal
Q}=0={\cal Q}{\cal P}$, details are shown in Appendix C.

##  III Properties of the non-equilibrium Langevin equation

###  III.1 General properties

Using the projection Eq. (II.3) in the generic GLE Eq. (II.3), we obtain the
explicit GLE

$\displaystyle\ddot{A}(\omega,t)=-K(t)(A(\omega,t)-\langle
A_{S}\rangle)-\int_{t_{0}}^{t}{\rm d}s\,\Gamma(s,t)\dot{A}(\omega,s)$
$\displaystyle+\int_{t_{0}}^{t}{\rm d}s\,\Gamma_{A}(s,t)(A(\omega,s)-\langle
A_{S}\rangle)$ $\displaystyle+F(\omega,t_{0},t)+(h(t)-h_{p})/M,$ (31)

the details of the derivation are shown in Appendix D. Eq. (III.1) is an exact
and explicit equation of motion for the Heisenberg observable $A(\omega,t)$
and is time-reversible, which is a consequence of the time-reversibility of
the underlying Hamilton and Liouville equations. Inspection of the GLE shows
that $F(\omega,t_{0},t)$ is the only term (except $h(t)$) in the GLE that
accounts for possible non-linearities (i.e. non-Gaussian contributions) in
$A(\omega,t)$. Thus, imposing $F(\omega,t_{0},t)$ to be a Gaussian variable
corresponds to a severe approximation for non-linear systems. On the other
hand, keeping the full non-Gaussian contributions of $F(\omega,t_{0},t)$ makes
Eq. (III.1) an exact description of the observable dynamics.

The first term in Eq. (III.1) is a force due to an effective harmonic
potential with a time-dependent potential stiffness $K(t)$ given by

$\displaystyle K(t)$ $\displaystyle=K_{0}+K_{1}(t)$ (32)

with

$\displaystyle K_{0}=\frac{\langle({\cal
L}_{0}A_{S}(\omega))^{2}\rangle}{\langle(A_{S}(\omega)-\langle
A_{S}\rangle)^{2}\rangle},$ (33) $\displaystyle
K_{1}(t)=\frac{-\beta(h(t)-h_{p})\langle(A_{S}(\omega)-\langle
A_{S}\rangle)({\cal
L}_{0}A_{S}(\omega))^{2}\rangle}{\langle(A_{S}(\omega)-\langle
A_{S}\rangle)^{2}\rangle}.$ (34)

The second term in Eq. (III.1) accounts for linear friction and depends on the
memory kernel given by

$\displaystyle\Gamma(s,t)=\Gamma_{0}(s,t)+\Gamma_{1}(s,t)$ (35)

with

$\displaystyle\Gamma_{0}(s,t)=\frac{\langle
F(\omega,s,s)F(\omega,s,t)\rangle}{\langle({\cal
L}_{0}A_{S}(\omega))^{2}\rangle},$ (36)
$\displaystyle\Gamma_{1}(s,t)=\frac{-\beta(h(s)-h_{p})\langle
F(\omega,s,t)({\cal L}_{0}A_{S}(\omega))^{2}\rangle}{\langle({\cal
L}_{0}A_{S}(\omega))^{2}\rangle}.$ (37)

There is also a positional memory term which is not present in the equilibrium
GLE in Eq. (1) and which involves the kernel function

$\displaystyle\Gamma_{A}(s,t)=$ (38)
$\displaystyle\frac{\beta(h(s)-h_{p})\langle
F(\omega,s,t)(A_{S}(\omega)-\langle A_{S}\rangle){\cal
L}_{0}A_{S}(\omega)\rangle}{\langle(A_{S}(\omega)-\langle
A_{S}\rangle)^{2}\rangle}.$

The last two terms in Eq. (III.1) are the complementary force
$F(\omega,t_{0},t)$ defined in Eq. (26) and the time-dependent force $h(t)$,
where the mass is given by

$\displaystyle M=\frac{1}{\beta\langle({\cal
L}_{0}A_{S}(\omega))^{2}\rangle}.$ (39)

For vanishing force $h(t)=0$, which renders the equilibrium scenario, and
choosing $h_{p}=0$, we see that $K_{1}(t)=\Gamma_{1}(s,t)=\Gamma_{A}(s,t)=0$
and we thus recover the standard form of the equilibrium Mori GLE in Eq. (1),
where in particular the friction kernel $\Gamma_{0}(s,t)$ is via Eq. (2)
related to the complementary force autocorrelation function (note that in Eqs.
(1) and (2) we have suppressed the phase-space dependence of the observable
$A$ and of the complementary force $F$). In contrast, if $h(t)-h_{p}\neq 0$,
we see that additional terms are present in the GLE and that for
$\Gamma_{1}(s,t)\neq 0$ the friction kernel $\Gamma(s,t)$ does not equal the
complementary force autocorrelation function. We first want to discuss whether
a non-zero $\Gamma_{1}(s,t)$ necessarily indicates the presence of non-
equilibrium effects.

An insightful scenario to address this question is one where the force
$h(t)=h_{0}$ is constant, in which case the complementary force and all memory
kernels become time homogeneous and can be written as
$F(\omega,s,t)=F(\omega,t-s)$, $\Gamma(s,t)=\Gamma(t-s)$,
$\Gamma_{A}(s,t)=\Gamma_{A}(t-s)$. Let us first discuss an unconfined system,
i.e., a system characterized by a diverging second moment,
$\langle(A_{S}(\omega)-\langle A_{S}\rangle)^{2}\rangle=\infty$. In this case
$K=0=\Gamma_{A}(t-s)$ and we must take $h_{p}=0$ in order to have a bounded
projection distribution $\rho_{p}(\omega)$ in Eq. (29). From Eq. (37) we see
that $\Gamma_{1}(t-s)$ can in general be non zero (as we will discuss in more
detail in the next section), in which case the total friction memory kernel
$\Gamma(t-s)$ does not equal the complementary force autocorrelation,
reflecting that an unconfined system under the influence of a constant force
dissipates energy and thus is a non-equilibrium system. In contrast, a
confined system that is characterized by a finite second moment
$\langle(A_{S}(\omega)-\langle A_{S}\rangle)^{2}\rangle$, is in the presence
of a constant force $h(t)=h_{0}$ an equilibrium system. But we see that for
$h_{0}\neq h_{p}$ the terms $K_{1}(t)$, $\Gamma_{1}(s,t)$ and
$\Gamma_{A}(s,t)$ do not necessarily vanish. In other words, unless we choose
as the projection distribution in Eq. (29) the equilibrium distribution, the
friction kernel in the Mori GLE $\Gamma(s,t)$ not necessarily equals the
complementary force autocorrelation. Thus, a non-zero $\Gamma_{1}(s,t)$ not
necessarily indicates a driven non-equilibrium system, but can also be
produced by choosing a specific projection distribution $\rho_{p}(\omega)$ and
characterize the approach of the system towards equilibrium. As we will
discuss next, $\Gamma_{1}(s,t)$ is predicted to vanish for Gaussian non-
equilibrium systems, so the interpretation of non-equilibrium GLEs has to be
done with care. Having made this important point, we from now on put
$h_{p}=0$.

Note that for equilibrium systems, the GLE parameters and the complementary
force can be extracted from simulation or experimental time series data by
various well-established techniques [39, 74, 75, 76, 77, 78, 47, 46]. Similar
extraction techniques for non-equlibrium time-series data will be discussed in
Sec. IV.

###  III.2 Gaussian versus non-Gaussian observables

In order to highlight the role played by non-Gaussian fluctuations of the
observable $A(\omega,t)$, we slightly rewrite the friction memory kernel in
Eq. (35) as

$\displaystyle\Gamma(s,t)=\Gamma_{0}(s,t)+\beta h(s)h(t)/M+\Gamma_{2}(s,t)$
(40)

with

$\displaystyle\Gamma_{2}(s,t)=\frac{-\beta
h(s)\langle(F(\omega,s,t)+h(t)/M)({\cal
L}_{0}A_{S}(\omega))^{2}\rangle}{\langle({\cal
L}_{0}A_{S}(\omega))^{2}\rangle}.$ (41)

For the potential memory term in Eq. (38) we choose the modified form

$\displaystyle\Gamma_{A}(s,t)=$ (42) $\displaystyle\frac{\beta
h(s)\langle(F(\omega,s,t)+h(t)/M)(A_{S}(\omega)-\langle A_{S}\rangle){\cal
L}_{0}A_{S}(\omega)\rangle}{\langle(A_{S}(\omega)-\langle
A_{S}\rangle)^{2}\rangle}.$

From the GLE (III.1) we see that $F(\omega,s,t)+h(t)/M$ is linear in
$A(\omega,\cdot)-\langle A_{S}\rangle$ with a projection time given by
$t_{0}=s$, we conclude from Eqs. (41) and (42) that $\Gamma_{2}(s,t)$ and
$\Gamma_{A}(s,t)$ are proportional to expectation values that are (at least)
cubic in $A(\omega,\cdot)-\langle A_{S}\rangle$. In other words, for an
observable $A(\omega,t)-\langle A_{S}\rangle$ that is Gaussian, the kernel
functions $\Gamma_{2}(s,t)$ and $\Gamma_{A}(s,t)$ (and also the potential
stiffness correction $K_{1}(t)$) vanish. In this case, we are thus led to the
simplified GLE, valid for Gaussian observables,

$\displaystyle\ddot{A}(\omega,t)=-K_{0}(A(\omega,t)-\langle
A_{S}\rangle)-\int_{t_{0}}^{t}{\rm d}s\,\Gamma_{G}(s,t)\dot{A}(\omega,s)$
$\displaystyle+F(\omega,t_{0},t)+h(t)/M$ (43)

where the Gaussian friction kernel is given by the autocorrelation of the sum
of the complementary and non-equilibrium forces according to

$\displaystyle\Gamma_{G}(s,t)=\frac{\langle(F(\omega,s,s)+h(s)/M)(F(\omega,s,t)\rangle+h(t)/M)}{\langle({\cal
L}_{0}A_{S}(\omega))^{2}\rangle}$ $\displaystyle=\frac{\langle
F(\omega,s,s)F(\omega,s,t)\rangle+h(t)h(s)/M^{2}}{\langle({\cal
L}_{0}A_{S}(\omega))^{2}\rangle}$ $\displaystyle=\frac{\langle
F(\omega,s,s)F(\omega,s,t)\rangle}{\langle({\cal
L}_{0}A_{S}(\omega))^{2}\rangle}+\beta h(t)h(s)/M.$ (44)

In the derivation of Eq. (III.2) we used that the non-equilibrium force $h(t)$
is phase-space independent. Thus, we conclude that a Gaussian non-equilibrium
variable is described by a GLE with a friction memory kernel that via Eq.
(III.2) is related to the autocorrelation of the total force acting on the
observable. One notes that Eqs. (III.2) and (III.2) are equivalent to the
standard Mori GLE formulation, Eqs. (1) and (2), provided the force in Eq. (1)
is interpreted as the sum of the complementary force and the non-equilibrium
force. It should be noted that an observable can be Gaussian while the entire
many-body system is non-Gaussian, meaning that other observables and system
coordinates may very well exhibit non-Gaussian fluctuations. Thus the class of
systems exhibiting Gaussian observables is a rather large one and includes for
example moving cells [56].

Another useful equation is derived by averaging the entire GLE Eq. (III.1)
over phase space, resulting in

$\displaystyle\ddot{a}(t)=\int{\rm
d}\omega\,\rho_{p}(\omega)\ddot{A}(\omega,t)$ (45)
$\displaystyle=-K(t)(a(t)-\langle A_{S}\rangle)-\int_{t_{0}}^{t}{\rm
d}s\,\Gamma(s,t)\dot{a}(s)$ $\displaystyle+\int_{t_{0}}^{t}{\rm
d}s\,\Gamma_{A}(s,t)(a(s)-\langle A_{S}\rangle)+h(t)/M,$

where the phase-space averaged Heisenberg observable is denoted as
$a(t)=\langle A(\omega,t)\rangle$ and $\langle F(\omega,t_{0},t)\rangle=0$ was
used. This equation describes how the mean observable $a(t)$ evolves in time
under the influence of the force $h(t)$, it therefore establishes the relation
between the mean observable $a(t)$ and the external force $h(t)$. It is exact
and valid beyond the linear-response approximation and can therefore be viewed
as a generalization of the linear-response FDT, as will be explored in detail
in Sect. IV.2.

To derive yet another GLE, we define the deviation of the Heisenberg
observable $A(\omega,t)$ around its mean as

$\displaystyle\Delta A(\omega,t)=A(\omega,t)-a(t),$ (46)

by subtracting Eqs. (III.1) and (45) the GLE for $\Delta A(\omega,t)$ follows
as

$\displaystyle\Delta\ddot{A}(\omega,t)=-K(t)\Delta
A(\omega,t)-\int_{t_{0}}^{t}{\rm d}s\,\Gamma(s,t)\Delta\dot{A}(\omega,s)$
$\displaystyle+\int_{t_{0}}^{t}{\rm d}s\,\Gamma_{A}(s,t)\Delta
A(\omega,s)+F(\omega,t_{0},t).$ (47)

From Eq. (III.2) we see that the complementary force $F(\omega,t_{0},t)$ is
linear in $\Delta A(\omega,\cdot)$; this means that if $\Delta A(\omega,t)$ is
a Gaussian variable, the kernel functions $\Gamma_{1}(s,t)$ and
$\Gamma_{A}(s,t)$ in Eqs. (37) and (38) vanish and we are thus led to the GLE
for the Gaussian deviatory non-equilibrium variable $\Delta A(\omega,t)$

$\displaystyle\Delta\ddot{A}(\omega,t)=-K_{0}\Delta
A(\omega,t)-\int_{t_{0}}^{t}{\rm d}s\,\Gamma_{0}(s,t)\Delta\dot{A}(\omega,s)$
$\displaystyle+F(\omega,t_{0},t),$ (48)

where the friction kernel $\Gamma_{0}(s,t)$ is given by the complementary
force autocorrelation via Eq. (36). Thus, the deviations of a Gaussian non-
equilibrium variable from its mean are described by a GLE of the form of Eq.
(1) that satisfies the standard FDT relation Eq. (2).

We have demonstrated in this section that a Gaussian observable is described
by a GLE that has the same form as the equilibrium GLE Eq. (1), in other
words, the stochastic behavior of an observable of a non-equilibrium system
differs from an equilibrium system only if the observable is non-Gaussian.
This is a very important finding, since many experimental observables are
Gaussian to a very good degree, for all such variables the standard
equilibrium Mori GLE in the form of Eq. (1), or, more precisely, Eq. (III.2),
is a valid description of the dynamics.

##  IV Stochastic non-equilibrium force

###  IV.1 Infinitesimal response

The GLE discussed so far is difficult to deal with in practice since it is
inhomogeneous in time; this is utterly expected and reflects the presence of
the time-dependent force $h(t)$ in the Hamiltonian but complicates the further
analysis. In many experimental scenarios the time evolution of the force
$h(t)$ is not known or unimportant, it therefore becomes useful to interpret
$h(t)$ as a stochastic variable that is only characterized by its first
moments

$\displaystyle\overline{1}=1,\,\,\,\,\overline{h(t)}=0,\,\,\,\,\overline{h(t)h(s)}=\sigma(t-s).$
(49)

The assumption of a vanishing first moment does not restrict the generality of
the model since (at least for bounded systems) we can subtract a constant from
$h(t)$ and move it into the equilibrium part of the Hamiltonian in Eq. (3).
Here $\sigma(t-s)$ denotes the force autocorrelation function, by only
defining the first two moments of $h(t)$ we are not necessarily implying that
the force is a Gaussian stochastic variable, as will become clearer later on.
By averaging the non-equilibrium GLE Eq. (III.1) over the non-equilibrium
force $h(t)$ we obtain the GLE

$\displaystyle\ddot{A}(\omega,t)=-K_{0}(A(\omega,t)-\langle
A_{S}\rangle)-\int_{t_{0}}^{t}{\rm
d}s\,\overline{\Gamma}(t-s)\dot{A}(\omega,s)$
$\displaystyle+\int_{t_{0}}^{t}{\rm
d}s\,\overline{\Gamma}_{A}(t-s)(A(\omega,s)-\langle
A_{S}\rangle)+\overline{F}(\omega,t-t_{0})+F_{\epsilon}(t),$ (50)

where we have added an infinitesimal generating force $F_{\epsilon}(t)$ by the
substitution $h(t)\rightarrow h(t)+F_{\epsilon}(t)$ in Eq. (III.1) prior to
averaging, which later will be used to derive the infinitesimal response of a
non-equilibrium system. When performing the force averaging, we use that the
observable $A(\omega,t)$ in Eq. (III.1) has no explicit dependence on $h(t)$,
which reflects that its time evolution is completely determined by its initial
value, its initial velocity and the GLE parameters, which do explicitly depend
on $h(t)$ (note that, in contrast, the solution $A(\omega,t)$ of the GLE Eq.
(III.1) does depend on $h(t)$ and the average of this solution over $h(t)$
obviously is not the same as the solution of Eq. (IV.1)). The force-averaged
memory kernels $\overline{\Gamma}(t-s)$ and $\overline{\Gamma}_{A}(t-s)$ and
the force-averaged complementary force $\overline{F}(\omega,t-t_{0})$ are
homogeneous in time, as is shown by perturbative operator expansion to leading
order in powers of $\sigma$ in Appendix E. The friction kernel follows from
Eq. (40) explicitly as

$\displaystyle\overline{\Gamma}(t-s)=\overline{\Gamma}_{0}(t-s)+\beta\sigma(t-s)/M+\overline{\Gamma}_{2}(t-s),$
(51)

where $\overline{\Gamma}_{0}(t-s)$ and $\overline{\Gamma}_{2}(t-s)$ denote the
force averages over Eqs. (36) and (41).

By additionally averaging over the phase space variable $\omega$ we obtain the
GLE for the phase-space and force-averaged Heisenberg observable $a(t)$,

$\displaystyle\ddot{a}(t)=-K_{0}(a(t)-\langle
A_{S}\rangle)-\int_{t_{0}}^{t}{\rm d}s\,\overline{\Gamma}(t-s)\dot{a}(s)$
$\displaystyle+\int_{t_{0}}^{t}{\rm
d}s\,\overline{\Gamma}_{A}(t-s)(a(s)-\langle A_{S}\rangle)+F_{\epsilon}(t).$
(52)

By Fourier transforming Eq. (IV.1) according to
$\tilde{a}(\nu)=\int_{-\infty}^{\infty} dte^{-\imath t\nu}a(t)$, we obtain the
infinitesimal response relation to first order in $F_{\epsilon}(t)$ as

$\displaystyle\tilde{a}(\nu)-2\pi\delta(\nu)\langle
A_{S}\rangle=\tilde{\chi}(\nu)\tilde{F}_{\epsilon}(\nu),$ (53)

where the Fourier-transformed response function is determined by

$\displaystyle
1/\tilde{\chi}(\nu)=K_{0}-\nu^{2}+\imath\nu\tilde{\overline{\Gamma}}^{+}(\nu)-\tilde{\overline{\Gamma}}^{+}_{A}(\nu)$
$\displaystyle=K_{0}-\nu^{2}+\imath\nu(\tilde{\overline{\Gamma}}_{0}^{+}(\nu)+\beta\tilde{\sigma}^{+}(\nu)/M+\tilde{\overline{\Gamma}}_{2}^{+}(\nu))-\tilde{\overline{\Gamma}}^{+}_{A}(\nu).$
(54)

In the last step we have inserted Eq. (51). Note that the response function
$\tilde{\chi}(\nu)$ accounts for the full non-linear dependence of the
observable on the non-equilibrium force $h(t)$, so it describes the response
to an infinitesimal force $F_{\epsilon}(t)$ in the presence of a finite (not
necessarily small) back-ground force $h(t)$. In deriving Eq. (IV.1) we have
shifted the projection time into the far past, $t_{0}\rightarrow-\infty$, and
have introduced causal or single-sided memory kernels and correlations, i.e.
$\overline{\Gamma}^{+}(t)=0$, $\overline{\Gamma}_{A}^{+}(t)=0$,
$\overline{\Gamma}_{0}^{+}(t)=0$, $\overline{\Gamma}_{2}^{+}(t)=0$,
$\sigma^{+}(t)=0$ for $t<0$. In fact, $\tilde{\chi}(\nu)$ can be obtained from
the Volterra equation for the correlation function that follows from the GLE,
as explained in Sec. IV.3. For a Gaussian observable, in which case
$\overline{\Gamma}_{A}(t)=0=\overline{\Gamma}_{2}(t)$ as discussed before, but
in the presence of a non-equilibrium stochastic force, $h(t)\neq 0$, we obtain
for the response function

$\displaystyle
1/\tilde{\chi}(\nu)=K_{0}-\nu^{2}+\imath\nu\tilde{\overline{\Gamma}}_{0}^{+}(\nu)+\imath\nu\beta\tilde{\sigma}_{+}(\nu)/M.$
(55)

In the absence of a non-equilibrium force, i.e. for $h(t)=0$, Eq. (IV.1)
reduces to the standard equilibrium response function

$\displaystyle
1/\tilde{\chi}(\nu)=K_{0}-\nu^{2}+\imath\nu\tilde{\overline{\Gamma}}_{0}^{+}(\nu).$
(56)

By comparison of Eqs. (IV.1), (55) and (56) we see that the presence of a
stochastic non-equilibrium force modifies the response function significantly
and adds terms that are determined by the force autocorrelation function
$\sigma(t)$.

###  IV.2 Fluctuation-dissipation theorem for finite non-equilibrium force
$h(t)$

In order to derive the non-equilibrium version of the FDT, we need to
calculate the two point correlation function

$\displaystyle C(t-t^{\prime})=\overline{\langle(A(\omega,t)-\langle
A_{S}\rangle)(A(\omega,t^{\prime})-\langle A_{S}\rangle)\rangle}$ (57)

for general times $t,t^{\prime}\geq t_{0}$, which is obtained by simultaneous
averaging over phase space and the non-equilibrium force $h(t)$. In Appendix F
we derive Eq. (57) and show that two-point correlation functions are generally
given by a phase-space average over products of Heisenberg variables, for this
we use the product propagation relation derived in Appendix G. In App. H we
show that the Fourier-transformed correlation function is to first order in an
expansion in powers of the second-order force moments
$\tilde{\overline{\Gamma}}_{0}(\nu)$ and $\tilde{\sigma}(\nu)$ given by

$\displaystyle\tilde{C}(\nu)=\tilde{\chi}(\nu)\tilde{\chi}(-\nu)\left[\langle({\cal
L}_{0}A_{S}(\omega))^{2}\rangle\tilde{\overline{\Gamma}}_{0}(\nu)+\tilde{\sigma}(\nu)/M^{2}\right],$
(58)

where we note that the force moments
$\tilde{\overline{\Gamma}}_{0}(\nu)=\tilde{\overline{\Gamma}}_{0}^{+}(\nu)+\tilde{\overline{\Gamma}}_{0}^{+}(-\nu)$
and $\tilde{\sigma}(\nu)=\tilde{\sigma}^{+}(\nu)+\tilde{\sigma}^{+}(-\nu)$ are
Fourier transforms of time-symmetrized functions. Using the expression for the
response function $\tilde{\chi}(\nu)$ in Eq. (IV.1), we can rewrite the
correlation function as

$\displaystyle\frac{\tilde{C}(\nu)}{ \langle({\cal
L}_{0}A_{S}(\omega))^{2}\rangle}=\frac{\tilde{\chi}(-\nu)}{\imath\nu}
-\frac{\tilde{\chi}(\nu)}{\imath\nu} $ (59)
$\displaystyle-\tilde{\chi}(\nu)\tilde{\chi}(-\nu)\left[\tilde{\overline{\Gamma}}_{2}^{+}(\nu)+\tilde{\overline{\Gamma}}_{2}^{+}(-\nu)-\frac{\tilde{\overline{\Gamma}}_{A}^{+}(\nu)-\tilde{\overline{\Gamma}}_{A}^{+}(-\nu)}{\imath\nu}\right].$

Since all time-domain kernel functions are real, we have
$\Re(\tilde{\Gamma}^{+}(\nu))=\Re(\tilde{\Gamma}^{+}(-\nu))$ and
$\Im(\tilde{\Gamma}^{+}(\nu))=-\Im(\tilde{\Gamma}^{+}(-\nu))$, where $\Re$ and
$\Im$ denote the real and imaginary parts of a complex number $X$ according to
$X=\Re X+\imath\Im X$. With this, we can rewrite the correlation function as

$\displaystyle\frac{\nu\tilde{C}(\nu)}{ 2\langle({\cal
L}_{0}A_{S}(\omega))^{2}\rangle}$ $\displaystyle=
-\Im\left(\tilde{\chi}(\nu)\right)$ (60)
$\displaystyle-\tilde{\chi}(\nu)\tilde{\chi}(-\nu)\left[\nu\Re(\tilde{\overline{\Gamma}}_{2}^{+}(\nu))-\Im(\tilde{\overline{\Gamma}}_{A}^{+}(\nu))\right],$

which is the FDT in the presence of a non-equilibrium stochastic force. As
discussed in Sec. III.2, if the observable $A(\omega,t)$ is Gaussian, the
memory kernel contributions $\tilde{\overline{\Gamma}}_{2}(\nu)$ and
$\tilde{\overline{\Gamma}}_{A}(\nu)$ vanish and we recover an equation that
resembles the standard FDT [8]

$\displaystyle\frac{\nu\tilde{C}(\nu)}{ 2\langle({\cal
L}_{0}A_{S}(\omega))^{2}\rangle}$ $\displaystyle=
-\Im\left(\tilde{\chi}(\nu)\right).$ (61)

However, the correlation function $\tilde{C}(\nu)$ and the response function
$\tilde{\chi}(\nu)$ do depend on the non-equilibrium force $h(t)$, so Eq. (61)
is a non-trivial generalization of the standard FDT (which does not depend on
the non-equilibrium force $h(t)$) to non-equilibrium Gaussian systems. In
Appendix I we derive the standard FDT by leading-order perturbation analysis
of the Heisenberg observable in Eq. (13). The comparison with this derivation
is instructive, since it deviates slightly from the text-book derivation of
the FDT and highlights the non-perturbative character of our non-equilibrium
GLE. Comparison of Eqs. (60) and (61) shows that our generalized FDT for non-
Gaussian non-equilibrium systems in Eq. (60) contains additional terms that
are essentially proportional to the non-equilibrium force moment $\sigma$ and
which result from a systematic first-order expansion in powers of $\sigma$, as
put forward in Appendix H.

A particularly transparent formulation of the non-equilibrium fluctuation
dissipation theorem Eq. (60) is given by

$\displaystyle-\frac{\nu\tilde{C}(\nu)/\langle({\cal
L}_{0}A_{S}(\omega))^{2}\rangle}{ 2\Im\left(\tilde{\chi}(\nu)\right)}= $
$\displaystyle\frac{\tilde{\overline{\Gamma}}_{0}(\nu)+\beta\tilde{\sigma}(\nu)/M}{2\left(\Re(\tilde{\overline{\Gamma}}^{+}(\nu))-\Im(\tilde{\overline{\Gamma}}_{A}^{+}(\nu))/\nu\right)}=1-\Xi(\nu),$
(62)

where in the first equation we used Eqs. (58) and (IV.1). Eq. (IV.2) contains
two alternative forms of our non-equilibrium FDT, the first in terms of the
correlation function and the imaginary part of the response function, which is
preferred when dealing with experimental data and when the non-equilibrium
force $h(t)$ is not known, the second in terms of the force correlations
$\tilde{\overline{\Gamma}}_{0}(\nu)$, $\tilde{\sigma}(\nu)$ and the real part
of the memory function, which is useful when the non-equilibrium force $h(t)$
is known and which will be further explained in Sec. IV.4. In Eq. (IV.2) we
defined the frequency-dependent non-equilibrium correction

$\displaystyle\Xi(\nu)=\frac{\tilde{\chi}(\nu)\tilde{\chi}(-\nu)}{
\Im\left(\tilde{\chi}(\nu)\right)}\left[\Re(\nu\tilde{\overline{\Gamma}}_{2}^{+}(\nu))-\Im(\tilde{\overline{\Gamma}}_{A}^{+}(\nu))\right],$
(63)

which has been previously introduced to quantify the departure from
equilibrium in biological non-equilibrium data [33]. Obviously, for Gaussian
observables, i.e. for
$\tilde{\overline{\Gamma}}_{2}^{+}(\nu)=0=\tilde{\overline{\Gamma}}_{A}^{+}(\nu)$
we have $\Xi(\nu)=0$. The expression for $\Xi(\nu)$ is proportional to the
factor $\tilde{\chi}(\nu)\tilde{\chi}(-\nu)/\Im(\tilde{\chi}(\nu))$, which
depends on the response function $\chi(\nu)$ and can be described in terms of
rather generic models. The non-Gaussian memory contributions
$\tilde{\overline{\Gamma}}_{2}^{+}(\nu)$ and
$\tilde{\overline{\Gamma}}_{A}^{+}(\nu)$ are conversely rather system-specific
and defy a generic approach.

###  IV.3  Response function from the Volterra equation

The standard way of extracting the memory kernel from time-series data is by
turning the stochastic GLE for the phase-space dependent observable
$A(\omega,t)$ into a non-stochastic integro-differential equation for the two-
point correlation function, which can be solved by Fourier transformation or
recursively after discretization in the time domain [39]. Here we show that
the same recipe also works for our non-equilibrium GLE. To proceed, we
multiply the GLE in Eq. (III.1) by $\dot{A}(\omega,t_{0})={\cal
L}_{0}A_{S}(\omega)$ and average over phase space $\omega$ and the non-
equilibrium force $h(t)$, by which we obtain the equation

$\displaystyle-\dddot{C}_{0}(t-t_{0})=K_{0}\dot{C}_{0}(t-t_{0})+\int_{0}^{t-t_{0}}{\rm
d}s\overline{\Gamma}(s)\ddot{C}_{0}(t-t_{0}-s)$
$\displaystyle-\int_{0}^{t-t_{0}}{\rm
d}s\overline{\Gamma}_{A}(s)\dot{C}_{0}(t-t_{0}-s)$ (64)

for the two point correlation function

$\displaystyle C_{0}(t-t_{0})=\overline{\langle(A(\omega,t_{0})-\langle
A_{S}\rangle)(A(\omega,t)-\langle A_{S}\rangle)\rangle}.$ (65)

Note that the two-point correlation function $C_{0}(t)$ defined here differs
from the one defined in Eq. (57) in that one of the times coincides with the
projection time $t_{0}$. This makes a fundamental difference, as will become
clear shortly. When deriving Eq. (IV.3), we used that the phase-space average
over the product of $\dot{A}(\omega,t_{0})$ and the forces $F(\omega,t_{0},t)$
or $h(t)$ in the GLE Eq. (III.1) vanishes and that the resulting equation
becomes homogeneous in time due to the average over the stochastic force
$h(t)$, as explained in Appendix E.

By performing a single-sided Fourier transform of Eq. (IV.3) while assuming
$t\geq t_{0}$ and defining $\tilde{C}_{0}^{+}(\nu)=\int_{0}^{\infty}
dte^{-\imath t\nu}C_{0}(t)$, we obtain the solution

$\displaystyle\tilde{C}_{0}^{+}(\nu)=\frac{\ddot{C}_{0}(0)\tilde{\chi}(\nu)}{\imath\nu}
+\frac{C_{0}(0)}{\imath\nu} $ (66)

in terms of the response function $\tilde{\chi}(\nu)$ defined in Eq. (IV.1).
This expression for $\tilde{C}_{0}^{+}(\nu)$ is equivalent to the expression
Eq. (61) for $\tilde{C}(\nu)$ in the Gaussian limit but for a non-Gaussian
observable differs from the full expression for $\tilde{C}(\nu)$, in other
words, $C_{0}(t)$ and $C(t)$ are completely different correlation function for
non-Gaussian non-equilibrium systems. It transpires that from the time-domain
correlation function $C_{0}(t)$, which can be straightforwardly obtained in
experiments or simulations by turning on the non-equilibrium force at time
$t_{0}$ (using that the distribution at time $t_{0}$ prior to application of
the non-equilibrium force equals the canonical projection distribution
$\rho_{p}(\omega)$ in Eq. (29)), the single-sided Fourier transform
$\tilde{C}_{0}^{+}(\nu)$ and the values $\ddot{C}_{0}(0)$ and $C_{0}(0)$
follow, from which $\tilde{\chi}(\nu)$ is determined via Eq. (66) by direct
inversion. The zero-frequency part of the response function can be furthermore
extracted by noting that $\int_{-\infty}^{\infty}{\rm
d}\nu\tilde{C}_{0}(\nu)=2\pi\ddot{C}_{0}(0)/(K_{0}-\tilde{\overline{\Gamma}}_{A}^{+}(0))$
where $\tilde{C}_{0}(\nu)=\tilde{C}_{0}^{+}(\nu)+\tilde{C}_{0}^{+}(-\nu)$, as
follows by residual calculus. Alternatively, instead of Fourier
transformation, Eq. (IV.3) can be recursively solved by discretization [39,
74, 46].

###  IV.4 Extracting non-equilibrium GLE parameters from time series data

In the previous section we showed that the correlation function $C_{0}(t)$,
defined in Eq. (65), is the solution of the differential equation Eq. (IV.3)
and can be used to calculate the response function $\tilde{\chi}(\nu)$ via
inversion of Eq. (66). This is a very practical method to obtain
$\tilde{\chi}(\nu)$ from an experimental or simulated correlation function. In
fact, from $\tilde{\chi}(\nu)$ all parameters of the non-equilibrium-force-
averaged GLE in Eq. (IV.1) can be obtained. To see this, we introduce the
running integral over the memory function
$\overline{G}_{A}(s)=\int_{s}^{t-t_{0}}{\rm
d}s^{\prime}\,\overline{\Gamma}_{A}(s^{\prime})$, with which Eq. (IV.1) can
after partial integration be rewritten as

$\displaystyle\ddot{A}(\omega,t)=-\left(K_{0}-\overline{G}_{A}(0)\right)(A(\omega,t)-\langle
A_{S}\rangle)$ (67) $\displaystyle-\int_{t_{0}}^{t}{\rm
d}s\,\left(\overline{\Gamma}(t-s)+\overline{G}_{A}(t-s)\right)\dot{A}(\omega,s)+\overline{F}(\omega,t-t_{0}).$

The GLE in Eq. (67) now depends on a single combined kernel function
$\overline{\Gamma}(t)+\overline{G}_{A}(t)$. Shifting the projection time into
the far past, $t_{0}\rightarrow-\infty$, the single-sided Fourier transform of
this kernel function can be related to the Fourier-transformed response
function $\tilde{\chi}(\nu)$ given in Eq. (IV.1) according to

$\displaystyle\tilde{\overline{\Gamma}}^{+}(\nu)+\tilde{\overline{G}}^{+}_{A}(\nu)=\frac{1}{\imath\nu}
\left(\frac{1}{\tilde{\chi}(\nu)} -\frac{1}{\tilde{\chi}(0)}+\nu^{2}\right).$
(68)

The stiffness of the effective harmonic potential that appears in Eq. (67) is
determined by the zero-frequency limit of the response function according to

$\displaystyle K_{0}-\overline{G}_{A}(0)=\frac{1}{\tilde{\chi}(0)},$ (69)

as follows from Eq. (IV.1). It transpires that all parameters of the GLE in
Eq. (67) can be derived from the response function $\tilde{\chi}(\nu)$
according to Eqs. (68) and (69). Furthermore, by inverting Eq. (67), the
complementary force trajectory $F(\omega,t_{0},t)$ can be calculated from a
trajectory of the observable $A(\omega,t)$ (which becomes the force-averaged
complementary force trajectory $\overline{F}(\omega,t-t_{0})$ only after
averaging over different realizations of $h(t)$). From that and using the
definition Eq. (36), we can by averaging over $h(t)$ calculate the
complementary force correlation $\overline{\Gamma}_{0}(t)$, which allows to
check whether time series data is of equilibrium or non-equilibrium nature.
For this we rewrite the FDT in Eq. (IV.2) using the definition of
$\overline{G}_{A}(t)$ as

$\displaystyle\frac{\tilde{\overline{\Gamma}}_{0}(\nu)+\beta\tilde{\sigma}(\nu)/M}{2\Re\left(\tilde{\overline{\Gamma}}^{+}(\nu)+\tilde{\overline{G}}_{A}^{+}(\nu)\right)}=1-\Xi(\nu),$
(70)

which for an equilibrium system simplifies to

$\displaystyle\frac{\tilde{\overline{\Gamma}}_{0}(\nu)}{2\Re\left(\tilde{\overline{\Gamma}}^{+}(\nu)\right)}=1.$
(71)

Knowing $\tilde{\overline{\Gamma}}_{0}(\nu)$ from the determined complementary
force trajectory and
$\tilde{\overline{\Gamma}}^{+}(\nu)+\tilde{\overline{G}}^{+}_{A}(\nu)$ from
Eq. (68), we can check whether Eq. (71) is violated and in that case infer
that the system is of non-equilibrium nature (the inverse conclusion cannot be
easily drawn since Eq. (71) could be satisfied even if the system is out of
equilibrium due to fortunate cancellation of terms in Eq. (70)). In
conclusion, knowledge of the response function $\tilde{\chi}(\nu)$, which can
be obtained from the correlation function $C_{0}(t)$ (determined by Eq.
(IV.3)) via Eq. (66), allows to detect whether a system is in equilibrium or
not. In addition, if the non-equilibrium force trajectory $h(t)$ is known, the
correction terms $\tilde{\overline{\Gamma}}_{2}(\nu)$ and
$\tilde{\overline{\Gamma}}_{A}(\nu)$, that appear in the non-equilibrium FDT
in Eq. (70) using the definition of $\Xi(\nu)$ in Eq. (63), can, according to
Eqs. (41) and (42), be calculated explicitly. This allows for an independent
check of the non-equilibrium FDT in Eq. (70).

More suitable for certain experimental systems where the non-equilibrium force
$h(t)$ can not be turned on and off at will, knowing the correlation function
$C(t)$ defined in Eq. (57) and the response function $\tilde{\chi}(\nu)$
defined via the infinitesimal response relation in Eq. (53), violation of the
equilibrium FDT in Eq. (61) can be checked and the non-Gaussian memory terms
in the non-equilibrium FDT in Eq. (60) can be extracted in an alternative
fashion.

###  IV.5 Joint observable distribution from path integrals

It remains to elucidate under which conditions the Heisenberg variable
$A(\omega,t)$ is described by a Gaussian process and thus the GLE Eq. (III.2)
is valid. For this we consider the two-point joint probability distribution of
$A(\omega,t)$, which is defined as

$\displaystyle\rho(A_{2},t_{2};A_{1},t_{1})=\overline{\langle\delta(A_{2}-A(\omega,t_{2})\delta(A_{1}-A(\omega,t_{1})\rangle}$
(72)

and involves averages over phase space $\omega$ and the non-equilibrium force
$h(t)$ according to Eqs. (28) and (49). For the delta functions we use the
Fourier representation and obtain

$\displaystyle\rho(A_{2},t_{2};A_{1},t_{1})=\int_{-\infty}^{\infty}\frac{{\rm
d}q_{1}}{2\pi} \int_{-\infty}^{\infty}\frac{{\rm d}q_{2}}{2\pi} $
$\displaystyle\times e^{\imath q_{1}A_{1}+\imath q_{2}A_{2}}
\overline{\langle\exp\left(-\imath q_{1}A(\omega,t_{1})-\imath
q_{2}A(\omega,t_{2})\right)\rangle}.$ (73)

Using the Fourier-transformed linear response relation in Eq. (53) and the
response function defined in Eq. (IV.1), we write the time-domain solution of
the GLE Eq. (III.1) as

$\displaystyle A(\omega,t)-\langle A_{S}\rangle=\int_{-\infty}^{\infty}{\rm
d}s\chi(s)(F(\omega,t-s)+h(t-s)/M).$ (74)

Note that in deriving Eq. (74) we have only partially averaged Eq. (III.1)
over $h(t)$ so that the equation becomes time-homogeneous but still has the
$h(t)$ source term. We also neglect the $h(t)$ dependence of the complementary
force so that it becomes homogeneous in time and can be written as
F$(\omega,t-s)$, the calculation is therefore first order in the force
correlation functions. The averages over the forces $F(\omega,t)$ and $h(t)$
in Eq. (IV.5) we express by Gaussian path integrals as

$\displaystyle\langle X(F,h)\rangle=\int_{-\infty}^{\infty}\frac{{\cal
D}F(\omega,\cdot)}{{\cal N}_{F}}X(F,h)$
$\displaystyle\exp\left(-\int_{-\infty}^{\infty}{\rm d}s{\rm
d}s^{\prime}\frac{F(\omega,s)F(\omega,s^{\prime})\Gamma_{0}^{-1}(s-s^{\prime})}{2{\langle({\cal
L}_{0}A_{S}(\omega))^{2}\rangle}}\right),$ (75)
$\displaystyle\overline{X(F,h)}=\int_{-\infty}^{\infty}\frac{{\cal
D}h(\cdot)}{{\cal N}_{h}}X(F,h)$
$\displaystyle\exp\left(-\int_{-\infty}^{\infty}{\rm d}s{\rm
d}s^{\prime}\frac{h(s)h(s^{\prime})\sigma^{-1}(s-s^{\prime})}{2}\right),$ (76)

where ${{\cal N}_{F}}$ and ${{\cal N}_{h}}$ are normalization constants and
$\Gamma_{0}^{-1}(s)$ and $\sigma^{-1}(s)$ are the inverse functions of the
Gaussian kernels defined in Eqs. (36) and (49) according to $\int{\rm
d}s\Gamma_{0}^{-1}(t-s)\Gamma_{0}(s-t^{\prime})=\delta(t^{\prime}-t)$. Non-
Gaussian fluctuations of $F(\omega,t)$ and $h(t)$ need not be explicitly
considered here as by construction they do not change their two-point
correlations. Performing the Gaussian path integrals, the two-point
distribution follows in matrix notation as

$\displaystyle\rho(A_{2},t_{2};A_{1},t_{1})=$
$\displaystyle\frac{\exp\left(-(A_{j}-\langle
A_{S}\rangle)I_{jk}^{-1}(A_{k}-\langle A_{S}\rangle)/2\right)}{\sqrt{\det 2\pi
I}},$ (77)

where the indices $j,k=1,2$ are summed over and the entries of the two-by-two
matrix

$\displaystyle I_{jk}=C(t_{j}-t_{k})$ (78)

are given by the two-point correlation function defined in Eq. (58), details
of the derivation are given in Appendix J. From Eq. (IV.5) we see that if the
forces $F(\omega,t)$ and $h(t)$ are described by general Gaussian processes,
as assumed in Eqs. (IV.5) and (IV.5), then also the observable
$A(\omega,t)-\langle A_{S}\rangle$ is a Gaussian process determined by the
correlation function $C(t)$ defined in Eq. (58). In this case, the GLE Eq.
(III.2) is valid, which has the structure of an equilibrium GLE. Conversely,
if one of the two forces $F(\omega,t)$ or $h(t)$ is non-Gaussian, then also
the observable $A(\omega,t)-\langle A_{S}\rangle$ is non-Gaussian and the GLE
in Eq. (III.1) applies, which does not satisfy the FDT in Eq. (III.2). To
reiterate this point, if the FDT in Eq. (III.2) is violated, this can be due
to non-Gaussian contributions in the complementary force $F(\omega,t)$ or in
the non-equilibrium force $h(t)$. It will in general be difficult to tell in
an experiment which of the forces generates the non-Gaussian behavior of the
observable $A(\omega,t)$, in particular since $F(\omega,t)$ depends on $h(t)$,
unless the non-equilibrium force $h(t)$ is generated externally and thus
explicitly known.

##  V Summary and Discussion

In the paper we have derived the non-equilibrium GLE from a many-body
Hamiltonian which contains a non-equilibrium time-dependent force $h(t)$ that
acts on a general phase-space-dependent observable $A_{S}(\omega)$. This is
the same Hamiltonian one uses for deriving the standard FDT, in other words,
the GLE we derive is conjugate to the standard FDT with one important
distinction: While the standard FDT describes the first-order response of the
time-dependent mean of the observable to the external force $h(t)$ and thus
the relation between the response function $\chi(t)$ and the two-point
correlation function $C(t)$ is independent of the force $h(t)$, we derive the
non-equilibrium GLE Eq. (III.1) non-perturbatively, i.e. exactly to all orders
in $h(t)$. Having pointed out the exact nature of our derived non-equilibrium
GLE, one should add that the kernel functions and the complementary force that
appear in the GLE Eq. (III.1) depend implicitly on the force trajectory via
the propagators. While linear and non-linear response theory describes how the
mean of an observable depends on the external force $h(t)$, the GLE is an
equation of motion for the fluctuating observable, it therefore opens up a
complementary field of applications and is particular relevant for the
description of time series data.

From the exact GLE Eq. (III.1) we infer that Gaussian non-equilibrium
observables are described by a GLE that takes the form of the equilibrium GLE
in Eq. (1). This is an important and rather non-trivial finding, as key
observables of many non-equilibrium systems are in fact Gaussian. For example,
the motion of cancer cells and algae [56] has been shown to be described by a
Gaussian process provided one looks at single-cell data. The correction terms
in the non-equilibrium GLE Eq. (III.1) that account for non-Gaussian effects
turn out to be three-point (and higher-order) correlation functions that
involve the complementary force.

The non-equilibrium GLE in Eq. (III.1) breaks time-homogeneity and thus is
difficult to deal with in practice. We therefore derive a GLE assuming that
the external force $h(t)$ is stochastic and defined by its second moment. The
preaveraging over $h(t)$ reinstalls time homogeneity of the GLE and allows us
to derive the non-equilibrium FDT to first order in a cumulant expansion in
terms of the external force and the complementary force. Similar to the GLE,
the non-equilibrium FDT for Gaussian observable has the same form as the
equilibrium FDT, only for non-Gaussian observables correction terms appear in
the non-equilibrium FDT that again are related to three-point (and higher-
order) correlation functions that involve the complementary force.

We also introduce different methods for extracting the parameters of our non-
equilibrium GLE from time series data. Here we distinguish between methods
that require knowledge of the external force trajectory $h(t)$, as is the case
for simulation data and for experiments where a system is perturbed by an
externally applied time-dependent force, for example by laser-traps or atomic-
force microscopes, and methods where the force trajetcory $h(t)$ is not known,
as is the case in most experiments on biological systems.

It is a widespread misconception that deviations from the equality in Eq. (2)
would signal a break-down of the standard FDT. The FDT relates the first-order
response function of the observable $A$ to the force $h(t)$ with the two-point
correlation function of $A$, thus, by definition, the standard FDT is
independent of $h(t)$, similar to Eq. (2) for an equilibrium system. For
finite force $h(t)\neq 0$ the relation between the friction kernel
$\Gamma(t-s)$ and the complementary force autocorrelation in Eq. (2) is
replaced by a relation that explicitly depends on $h(t)$. Likewise, the
relation between the non-linear response of $A$ to an infinitesimal field
increase $h(t)+MF_{\epsilon}(t)$ and the fluctuations of $A$ becomes dependent
on $h(t)$ and in fact involves three-point correlation functions of $A$, as we
show in this paper.

A second widespread misconception is that the Mori GLE Eq. (1) is approximate
and only holds on the Gaussian level. The opposite is true, Eq. (1) is exact
and any non-linear properties that the observable $A$ might have are
accurately represented by non-Gaussian contributions of $F$. This means that
the non-equilibrium GLE we derive, Eq. (III.1), also applies to non-Gaussian
non-equilibrium observables, provided the non-Gaussian contributions from the
complementary force are correctly included.

One conclusion from our work is that for Gaussian systems, it does not make
sense to consider a GLE that violates the FDT, this nicely explains our
previous finding that for Gaussian systems, there is no way of detecting non-
equilibrium properties from time-series data [56].

We mention in passing that the external force h(t) in general performs work on
the system, meaning that the total energy of the system will in general
increase with time. Since the GLE Eq. (III.1) is exact, it correctly takes
into account these transient energetic effects.

The time-dependent Hamiltonian Eq. (3) includes only the linear coupling
between the time-dependent force $h(t)$ and the observable $A_{S}(\omega)$,
more complicated coupling terms are conceivable. The motivation for studying
this simple time-dependent Hamiltonian is that it leads to the standard FDT,
our non-equilibrium GLE thus is conjugate to the standard FDT. In the future,
it will be interesting to derive GLEs from more complex time-dependent
Hamiltonians.

###### Acknowledgements.

We acknowledge support by Deutsche Forschungsgemeinschaft Grant CRC 1449
”Dynamic Hydrogels at Biointerfaces”, ID 431232613 Project A03, by the ERC
Advanced Grant 835117 NoMaMemo and by the Infosys Foundation.

##  Appendix A Derivative of Heisenberg propagator with respect to initial
time

Using the Heavyside function $\theta(t)$, defined as $\theta(t)=1$ for $t>0$
and $\theta(t)=0$ for $t<0$, the time-ordered operator exponential in the
Heisenberg picture Eq. 14 can be rewritten as

$\displaystyle\exp_{H}\left(\int_{t_{0}}^{t}{\rm d}s{\cal L}(s)\right)$
$\displaystyle=1+\sum_{n=1}^{\infty}\int_{t_{0}}^{t}{\rm
d}t_{1}\int_{t_{0}}^{t}{\rm d}t_{2}\theta(t_{1}-t_{2})\int_{t_{0}}^{t}{\rm
d}t_{3}\theta(t_{2}-t_{3})\cdots\int_{t_{0}}^{t}{\rm
d}t_{n}\theta(t_{n-1}-t_{n}){\cal L}(t_{n})\cdots{\cal L}(t_{3}){\cal
L}(t_{2}){\cal L}(t_{1}).$ (79)

We now reorder the integration variables to obtain

$\displaystyle\exp_{H}\left(\int_{t_{0}}^{t}{\rm d}s{\cal L}(s)\right)=$
$\displaystyle 1+\sum_{n=1}^{\infty}\int_{t_{0}}^{t}{\rm
d}t_{n}\int_{t_{0}}^{t}{\rm d}t_{n-1}\theta(t_{n-1}-t_{n})\int_{t_{0}}^{t}{\rm
d}t_{n-2}\theta(t_{n-2}-t_{n-1})\cdots\int_{t_{0}}^{t}{\rm
d}t_{1}\theta(t_{1}-t_{2}){\cal L}(t_{n}){\cal L}(t_{n-1}){\cal
L}(t_{n-2})\cdots{\cal L}(t_{1}),$ (80)

note that the Liouville operators cannot be reordered since ${\cal L}(t_{1})$
and ${\cal L}(t_{2})$ in general do not commute for $t_{1}\neq t_{2}$. We now
replace the Heavyside functions by the appropriate integration boundaries and
obtain

$\displaystyle\exp_{H}\left(\int_{t_{0}}^{t}{\rm d}s{\cal
L}(s)\right)=1+\sum_{n=1}^{\infty}\int_{t_{0}}^{t}{\rm
d}t_{n}\int_{t_{n}}^{t}{\rm d}t_{n-1}\int_{t_{n-1}}^{t}{\rm
d}t_{n-2}\cdots\int_{t_{2}}^{t}{\rm d}t_{1}{\cal L}(t_{n}){\cal
L}(t_{n-1}){\cal L}(t_{n-2})\cdots{\cal L}(t_{1}).$ (81)

From this expression Eq. (16) follows directly.

##  Appendix B Derivation of the time-dependent Dyson operator expansion

The operator expansion for two general time-dependent operators ${\cal V}(t)$
and ${\cal W}(t)$ reads

$\displaystyle\exp_{H}\left(\int_{t_{0}}^{t}{\cal V}(s)+{\cal W}(s){\rm
d}s\right)=\exp_{H}\left(\int_{t_{0}}^{t}{\rm d}s{\cal
V}(s)\right)+\int_{t_{0}}^{t}{\rm d}s\exp_{H}\left(\int_{t_{0}}^{s}{\rm
d}s^{\prime}{\cal V}(s^{\prime})\right){\cal
W}(s)\exp_{H}\left(\int_{s}^{t}{\cal V}(s^{\prime})+{\cal W}(s^{\prime}){\rm
d}s^{\prime}\right).$ (82)

To prove this relation, we use Eqs. (15) and (16) to obtain from Eq. (82)

$\displaystyle\exp_{H}\left(\int_{t_{0}}^{t}{\cal V}(s)+{\cal W}(s){\rm
d}s\right)=\exp_{H}\left(\int_{t_{0}}^{t}{\rm d}s{\cal
V}(s)\right)-\int_{t_{0}}^{t}{\rm d}s\frac{{\rm d}}{ {\rm
d}s}\exp_{H}\left(\int_{t_{0}}^{s}{\rm d}s^{\prime}{\cal
V}(s^{\prime})\right)\exp_{H}\left(\int_{s}^{t}{\cal V}(s^{\prime})+{\cal
W}(s^{\prime}){\rm d}s^{\prime}\right).$ (83)

Now the integral can be performed and the equality is obtained. An alternative
operator expansion relation reads

$\displaystyle\exp_{H}\left(\int_{t_{0}}^{t}{\cal V}(s)+{\cal W}(s){\rm
d}s\right)=\exp_{H}\left(\int_{t_{0}}^{t}{\rm d}s{\cal
V}(s)\right)+\int_{t_{0}}^{t}{\rm d}s\exp_{H}\left(\int_{t_{0}}^{s}{\cal
V}(s)+{\cal W}(s){\rm d}s^{\prime}\right){\cal
W}(s)\exp_{H}\left(\int_{s}^{t}{\rm d}s^{\prime}{\cal V}(s^{\prime})\right),$
(84)

which can be proven analogously. Choosing ${\cal V}(s)={\cal Q}{\cal L}(s)$
and ${\cal W}(s)={\cal P}{\cal L}(s)$ in Eq. (84) we obtain Eq. (24).

##  Appendix C Derivation of essential Mori projection properties

For the following derivations it is useful to split the Mori projection
operator in Eq. (II.3) into three parts according to

$\displaystyle{\cal P}B(\omega,t)={\cal P}_{1}B(\omega,t)+{\cal
P}_{2}B(\omega,t)+{\cal P}_{3}B(\omega,t)$ (85)

with

$\displaystyle{\cal P}_{1}B(\omega,t)=\langle B(\omega,t)\rangle,$ (86)
$\displaystyle{\cal P}_{2}B(\omega,t)=\frac{\langle B(\omega,t){\cal L}_{0}
A_{S}(\omega)\rangle}{\langle({\cal L}_{0}A_{S}(\omega))^{2}\rangle}{\cal
L}_{0}A_{S}(\omega),$ (87) $\displaystyle{\cal P}_{3}B(\omega,t)=\frac{\langle
B(\omega,t)(A_{S}(\omega)-\langle
A_{S}\rangle)\rangle}{\langle(A_{S}(\omega)-\langle
A_{S}\rangle)^{2}\rangle}(A_{S}(\omega)-\langle A_{S}\rangle).$ (88)

The linearity of the Mori projection, i.e., the fact that for two arbitrary
observables $B(\omega,t)$ and $C(\omega,t^{\prime})$ the property ${\cal P}
(c_{1}B(\omega,t)+c_{2}C(\omega,t^{\prime}))=c_{1}{\cal
P}B(\omega,t)+c_{2}{\cal P}C(\omega,t^{\prime})$ holds, is self-evident,
${\cal Q}$ is also easily seen to be linear.

The idempotency of ${\cal P}$, i.e., the fact that ${\cal P}^{2}={\cal P}$, is
not self-evident and will be proven. We split the proof in three parts. First,

$\displaystyle{\cal P}{\cal P}_{1}B(\omega,t)$ $\displaystyle=\langle\langle
B(\omega,t)\rangle\rangle+\langle B(\omega,t)\rangle\frac{\langle{\cal L}_{0}
A_{S}(\omega)\rangle}{\langle({\cal L}_{0}A_{S}(\omega))^{2}\rangle}{\cal
L}_{0}A_{S}(\omega)+\langle
B(\omega,t)\rangle\frac{\langle(A_{S}(\omega)-\langle
A_{S}\rangle)\rangle}{\langle(A_{S}(\omega)-\langle
A_{S}\rangle)^{2}\rangle}(A_{S}(\omega)-\langle A_{S}\rangle)$
$\displaystyle={\cal P}_{1}B(\omega,t).$ (89)

For the first and third terms we used that $\langle\langle
B(\omega,t)\rangle\rangle=\langle B(\omega,t)\rangle$, which holds since the
probability distribution in Eq. (29) is normalized. For the second term we
assume that $A_{S}$ is a function of position only, i.e.
$A_{S}(\omega)=A_{S}(\bf R)$, such that ${\cal L}_{0} A_{S}(\omega)$ is linear
in the momenta and the average $\langle{\cal L}_{0} A_{S}(\omega)\rangle$
vanishes.

Second,

$\displaystyle{\cal P}{\cal P}_{2}B(\omega,t)$ $\displaystyle=\frac{\langle
B(\omega,t){\cal L}_{0}A_{S}(\omega)\rangle}{\langle({\cal
L}_{0}A_{S}(\omega))^{2}\rangle}\left[\langle{\cal
L}_{0}A_{S}(\omega)\rangle+\frac{\langle({\cal
L}_{0}A_{S}(\omega))^{2}\rangle}{\langle({\cal
L}_{0}A_{S}(\omega))^{2}\rangle}{\cal
L}_{0}A_{S}(\omega)+\frac{\langle(A_{S}(\omega)-\langle A_{S}\rangle){\cal
L}_{0}A_{S}(\omega)\rangle}{\langle(A_{S}(\omega)-\langle
A_{S}\rangle)^{2}\rangle}(A_{S}(\omega)-\langle A_{S}\rangle)\right] $
$\displaystyle={\cal P}_{2}B(\omega,t),$ (90)

where again we used that $A_{S}(\omega)=A_{S}(\bf R)$ such that ${\cal L}_{0}
A_{S}(\omega)$ is linear in the momenta and the averages $\langle{\cal L}_{0}
A_{S}(\omega)\rangle$ and $\langle(A_{S}(\omega)-\langle A_{S}\rangle){\cal
L}_{0}A_{S}(\omega)\rangle$ vanish.

Third,

$\displaystyle{\cal P}{\cal P}_{3}B(\omega,t)$ $\displaystyle=\frac{\langle
B(\omega,t)(A_{S}(\omega)-\langle
A_{S}\rangle)\rangle}{\langle(A_{S}(\omega)-\langle A_{S}\rangle)^{2}\rangle}$
$\displaystyle\times\left[\langle A_{S}(\omega)-\langle
A_{S}\rangle\rangle+\frac{\langle(A_{S}(\omega)-\langle A_{S}\rangle){\cal
L}_{0}A_{S}(\omega)\rangle}{\langle({\cal
L}_{0}A_{S}(\omega))^{2}\rangle}{\cal
L}_{0}A_{S}(\omega)+\frac{\langle(A_{S}(\omega)-\langle
A_{S}\rangle)^{2}\rangle}{\langle(A_{S}(\omega)-\langle
A_{S}\rangle)^{2}\rangle}(A_{S}(\omega)-\langle A_{S}\rangle)\right] $
$\displaystyle={\cal P}_{3}B(\omega,t),$ (91)

where again we used that $A_{S}(\omega)=A_{S}(\bf R)$. Adding Eqs. (C), (C),
(C) we see that ${\cal P}^{2}={\cal P}({\cal P}_{1}+{\cal P}_{2}+{\cal
P}_{3})={\cal P}_{1}+{\cal P}_{2}+{\cal P}_{3}={\cal P}$ and thus ${\cal P}$
is idempotent. From the idempotency of ${\cal P}$ it follows that ${\cal Q}$
is also idempotent, to prove this one writes

$\displaystyle{\cal Q}^{2}B(\omega,t)=(1-{\cal P})^{2}B(\omega,t)=(1-2{\cal
P}+{\cal P}^{2})B(\omega,t)=(1-{\cal P})B(\omega,t)={\cal Q}B(\omega,t).$ (92)

The self-adjointedness of ${\cal P}$, Eq. (30), is straightforwardly proven by
writing

$\displaystyle\langle C(\omega,t){\cal P}B(\omega,t^{\prime})\rangle=$ (93)
$\displaystyle\langle C(\omega,t)\rangle\langle
B(\omega,t^{\prime})\rangle+\langle C(\omega,t){\cal
L}_{0}A_{S}(\omega)\rangle\frac{\langle B(\omega,t^{\prime}){\cal
L}_{0}A_{S}(\omega)\rangle}{\langle({\cal
L}_{0}A_{S}(\omega))^{2}\rangle}+\langle C(\omega,t)(A_{S}(\omega)-\langle
A_{S}\rangle)\rangle\frac{\langle B(\omega,t^{\prime})(A_{S}(\omega)-\langle
A_{S}\rangle)\rangle}{\langle(A_{S}(\omega)-\langle A_{S}\rangle)^{2}\rangle}$
(94) $\displaystyle=\langle B(\omega,t^{\prime}){\cal P}C(\omega,t)\rangle.$
(95)

By using ${\cal Q}=1-{\cal P}$ we see straightforwardly that ${\cal Q}$ is
also self-adjoint.

Using similar arguments as above, one can show that ${\cal P}c$= c, ${\cal
P}(A_{S}(\omega)-\langle A_{S}\rangle)=(A_{S}(\omega)-\langle A_{S}\rangle)$,
${\cal P}{\cal L}_{0}A_{S}(\omega)={\cal L}_{0}A_{S}(\omega)$, from which
follows that also ${\cal P}A_{S}(\omega)=A_{S}(\omega)$. From these relations
we can directly conclude that ${\cal Q}c$= 0, ${\cal Q}(A_{S}(\omega)-\langle
A_{S}\rangle)=0$, ${\cal Q}{\cal L}_{0}A_{S}(\omega)=0$, and also ${\cal
Q}A_{S}(\omega)=0$.

From the idempotency of ${\cal P}$ or ${\cal Q}$ we follow that ${\cal P}{\cal
Q}={\cal P}(1-{\cal P})={\cal P}-{\cal P}^{2}=0$ and, similarly, ${\cal
Q}{\cal P}=0$, thus, the operators ${\cal P}$ and ${\cal Q}$ are orthogonal to
each other.

##  Appendix D Derivation of non-equilibrium GLE

We consider the first term in Eq. (II.3), which can be split into two terms
and reads, apart from the propagator in front,

$\displaystyle{\cal P}{\cal L}(t){\cal L}_{0}A_{S}(\omega)={\cal P}{\cal
L}_{0}^{2}A_{S}(\omega)-h(t){\cal P}\Delta{\cal L}{\cal L}_{0}A_{S}(\omega).$
(96)

We apply the projection operator Eq. (II.3) on the first term in Eq. (96),
which generates three contributions. The first contribution is given by

$\displaystyle\langle{\cal L}_{0}^{2}A_{S}(\omega)\rangle=-\beta
h_{p}\langle({\cal L}_{0}A_{S}(\omega))^{2}\rangle,$ (97)

where we used that ${\cal L}_{0}$ is anti-self-adjoint and that ${\cal
L}_{0}\rho_{p}(\omega)=\beta h_{p}\rho_{p}(\omega){\cal L}_{0}A_{S}(\omega)$.
Using the same properties of ${\cal L}_{0}$, the second contribution can be
written as (apart from the normalization factor)

$\displaystyle\langle({\cal L}_{0}A_{S}(\omega)){\cal
L}_{0}^{2}A_{S}(\omega)\rangle=-\beta h_{p}\langle({\cal
L}_{0}A_{S}(\omega))^{3}\rangle/2.$ (98)

The third contribution (apart from the normalization factor) follows as

$\displaystyle\langle((A_{S}(\omega)-\langle A_{S}\rangle){\cal
L}_{0}^{2}A_{S}(\omega)\rangle$ $\displaystyle=-\langle({\cal
L}_{0}A_{S}(\omega))^{2}\rangle-\beta h_{p}\langle(A_{S}(\omega)-\langle
A_{S}\rangle)({\cal L}_{0}A_{S}(\omega))^{2}\rangle.$ (99)

We now apply the projection operator Eq. (II.3) on the second term in Eq.
(96), which again generates three contributions. The first contribution is
given by

$\displaystyle-h(t)\langle\Delta{\cal L}{\cal L}_{0}A_{S}(\omega)\rangle=\beta
h(t)\langle({\cal L}_{0}A_{S}(\omega))^{2}\rangle,$ (100)

where we used that $\Delta{\cal L}$ is anti-self-adjoint and that $\Delta{\cal
L}\rho_{p}(\omega)=\beta\rho_{p}(\omega){\cal L}_{0}A_{S}(\omega)$. Using the
same properties of $\Delta{\cal L}$, the second contribution can be written as
(apart from the normalization factor)

$\displaystyle-h(t)\langle({\cal L}_{0}A_{S}(\omega))\Delta{\cal L}{\cal
L}_{0}A_{S}(\omega)\rangle=\beta h(t)\langle({\cal
L}_{0}A_{S}(\omega))^{3}\rangle/2.$ (101)

The third contribution (apart from the normalization factor) follows as

$\displaystyle-h(t)\langle((A_{S}(\omega)-\langle A_{S}\rangle)\Delta{\cal
L}{\cal L}_{0}A_{S}(\omega)\rangle$ $\displaystyle=\beta
h(t)\langle(A_{S}(\omega)-\langle A_{S}\rangle)({\cal
L}_{0}A_{S}(\omega))^{2}\rangle.$ (102)

Again assuming that the observable is a function of position only,
$A_{S}(\omega)=A_{S}(\bf R)$, as we did in Appendix. C when we derived the
idempotency of the projection operator ${\cal P}$, we see that Eqs. (98) and
(101) vanish because $\langle({\cal L}_{0}A_{S}(\omega))^{3}\rangle$ is odd in
the momenta.

Combining the results in Eqs. (97), (D), (100), (D), the first term in Eq.
(II.3) reads

$\displaystyle\exp_{H}\left(\int_{t_{0}}^{t}{\rm d}s{\cal L}(s)\right){\cal
P}{\cal L}(t){\cal L}_{0}A_{S}(\omega)$
$\displaystyle=-K(t)(A(\omega,t)-\langle A_{S}\rangle)+(h(t)-h_{p})/M,$ (103)

where $K(t)$ and $M$ are defined in Eqs. (32) and (39).

We now consider the last term in Eq. (II.3), which reads, without the time
integral and the propagator in front, ${\cal P}{\cal L}(s)F(\omega,s,t)$. The
projection operator Eq. (II.3) generates three contributions. The first
contribution is given by

$\displaystyle\langle{\cal L}(s)F(\omega,s,t)\rangle=\beta(h(s)-h_{p})\langle
F(\omega,s,t){\cal L}_{0}A_{S}(\omega)\rangle$
$\displaystyle=\beta(h(s)-h_{p})\langle F(\omega,s,t){\cal Q}{\cal
L}_{0}A_{S}(\omega)\rangle$ $\displaystyle=0,$ (104)

where in the first line we used that ${\cal L}(s)$ is anti-self-adjoint and
that ${\cal L}(s)\rho_{p}(\omega)=\beta(h_{p}-h(s))\rho_{p}(\omega){\cal
L}_{0}A_{S}(\omega)$, in the second line that ${\cal Q}$ is idempotent and
self-adjoint, and in the third line that ${\cal Q}{\cal L}_{0}A_{S}(\omega)=0$
(as derived in Appendix C). The second contribution can be written as (apart
from the normalization factor)

$\displaystyle\langle({\cal L}_{0}A_{S}(\omega)){\cal
L}(s)F(\omega,s,t)\rangle$ $\displaystyle=-\langle
F(\omega,s,t)\rangle\rangle{\cal L}(s){\cal L}_{0}A_{S}(\omega)$
$\displaystyle+\beta(h(s)-h_{p})\langle({\cal
L}_{0}A_{S}(\omega))^{2}F(\omega,s,t)\rangle$ $\displaystyle=-\langle
F(\omega,s,s)F(\omega,s,t)\rangle\rangle$
$\displaystyle+\beta(h(s)-h_{p})\langle({\cal
L}_{0}A_{S}(\omega))^{2}F(\omega,s,t)\rangle,$ (105)

where in the first equation we used the same properties of ${\cal L}(s)$ as
before and in the second equation that ${\cal Q}$ is idempotent and self-
adjoint and the definition of the complementary force in Eq. (26).

The third contribution can be written as (apart from the normalization factor)

$\displaystyle\langle(A_{S}(\omega)-\langle A_{S}\rangle){\cal
L}(s)F(\omega,s,t)\rangle$ $\displaystyle=-\langle F(\omega,s,t){\cal
L}(s)(A_{S}(\omega)-\langle A_{S}\rangle)\rangle$
$\displaystyle+\beta(h(s)-h_{p})\langle F(\omega,s,t)(A_{S}(\omega)-\langle
A_{S}\rangle){\cal L}_{0}A_{S}(\omega)\rangle$
$\displaystyle=\beta(h(s)-h_{p})\langle F(\omega,s,t)(A_{S}(\omega)-\langle
A_{S}\rangle){\cal L}_{0}A_{S}(\omega)\rangle,$ (106)

where in the first equation we used the same properties of ${\cal L}(s)$ as
before and in the second equation that ${\cal Q}$ is idempotent and self-
adjoint and that $\Delta{\cal L}A_{S}(\omega)=0$.

Combining the results in Eqs. (D), (D), (D), (D), the expression for the
general GLE in Eq. (II.3) leads to the explicit GLE in Eq. (III.1).

##  Appendix E Derivation of time-homogeneous memory kernels for stochastic
force

In this section we show that the memory kernels of the GLE in Eq. (III.1)
become time-homogeneous for a stochastic non-equilibrium force as defined by
Eq. 49. The proof is done to quadratic order in powers of the force $h(t)$. We
start with the operator expansion for two general time-dependent operators
${\cal V}(t)$ and ${\cal W}(t)$, Eq. 82, and note that we can construct a
systematic perturbative expansion in the operator ${\cal W}(t)$ by writing

$\displaystyle\exp_{H}\left(\int_{t_{0}}^{t}{\cal V}(s)+{\cal W}(s){\rm
d}s\right)$ $\displaystyle=\exp_{H}\left(\int_{t_{0}}^{t}{\rm d}s{\cal
V}(s)\right)+\int_{t_{0}}^{t}{\rm d}s\exp_{H}\left(\int_{t_{0}}^{s}{\rm
d}s^{\prime}{\cal V}(s^{\prime})\right){\cal
W}(s)\exp_{H}\left(\int_{s}^{t}{\cal V}(s^{\prime})+{\cal W}(s^{\prime}){\rm
d}s^{\prime}\right)$ $\displaystyle=\exp_{H}\left(\int_{t_{0}}^{t}{\rm
d}s{\cal V}(s)\right)+\int_{t_{0}}^{t}{\rm
d}s\exp_{H}\left(\int_{t_{0}}^{s}{\rm d}s^{\prime}{\cal
V}(s^{\prime})\right){\cal W}(s)\exp_{H}\left(\int_{s}^{t}{\cal
V}(s^{\prime}){\rm d}s^{\prime}\right)+{\cal O}({\cal W}^{2}).$ (107)

Expressions valid to higher order in ${\cal W}$ can be constructed by
recursively inserting the second line of Eq. (E) into the last operator
exponential in the first line.

Following this method, the complementary force can be expanded in powers of
the non-equilibrium force as

$\displaystyle F(\omega,t_{0},t)$ $\displaystyle\equiv\exp_{H}\left({\cal
Q}\int_{t_{0}}^{t}{\rm d}s^{\prime}{\cal L}(s^{\prime})\right){\cal Q}{\cal
L}(t){\cal L}_{0}A_{S}(\omega)$ $\displaystyle=\exp_{H}\left({\cal
Q}\int_{t_{0}}^{t}{\cal L}_{0}-h(s^{\prime})\Delta{\cal L}{\rm
d}s^{\prime}\right){\cal Q}({\cal L}_{0}-h(t)\Delta{\cal L}){\cal
L}_{0}A_{S}(\omega)$ $\displaystyle=\left[e^{(t-t_{0}){\cal QL}_{0}}
-\int_{t_{0}}^{t}{\rm d}s^{\prime}e^{(s^{\prime}-t_{0}){\cal
QL}_{0}}h(s^{\prime}){\cal Q}\Delta{\cal L}e^{(t-s^{\prime}){\cal
QL}_{0}}\right] {\cal Q}({\cal L}_{0}-h(t)\Delta{\cal L}){\cal
L}_{0}A_{S}(\omega)+{\cal O}(h^{2}),$ (108)

where higher-order terms can be recursively derived. We calculate here
exemplarily the force average over the non-equilibrium friction memory kernel
defined in Eq. (37)

$\displaystyle\overline{\Gamma}_{1}(s,t)=\frac{-\beta\overline{h(s)\langle
F(\omega,s,t)({\cal L}_{0}A_{S}(\omega))^{2}\rangle}}{\langle({\cal
L}_{0}A_{S}(\omega))^{2}\rangle}=\frac{-\beta(X_{1}(s,t)+X_{2}(s,t))}{\langle({\cal
L}_{0}A_{S}(\omega))^{2}\rangle},$ (109)

which we split into two terms using the expansion in Eq. (E). The first term
reads, after performing the average over $h(t)$,

$\displaystyle X_{1}(s,t)=-\sigma(t-s)\left\langle({\cal
L}_{0}A_{S}(\omega))^{2}e^{(t-s){\cal QL}_{0}}{\cal Q}\Delta{\cal L}{\cal
L}_{0}A_{S}(\omega)\right\rangle,$ (110)

which clearly is a function of $t-s$ only. The second term reads

$\displaystyle X_{2}(s,t)=$ $\displaystyle-\left\langle({\cal
L}_{0}A_{S}(\omega))^{2}\int_{s}^{t}{\rm d}s^{\prime}e^{(s^{\prime}-s){\cal
QL}_{0}}\sigma(s^{\prime}-s){\cal Q}\Delta{\cal L}e^{(t-s^{\prime}){\cal
QL}_{0}}{\cal Q}{\cal L}_{0}^{2}A_{S}(\omega)\right\rangle$
$\displaystyle=-\left\langle({\cal L}_{0}A_{S}(\omega))^{2}\int_{0}^{t-s}{\rm
d}\tilde{s}^{\prime}e^{\tilde{s}^{\prime}{\cal
QL}_{0}}\sigma(\tilde{s}^{\prime}){\cal Q}\Delta{\cal
L}e^{(t-s-\tilde{s}^{\prime}){\cal QL}_{0}}{\cal Q}{\cal
L}_{0}^{2}A_{S}(\omega)\right\rangle,$ (111)

which also clearly is a function of $t-s$ only and where we have used the
change of integration variable according to $\tilde{s}^{\prime}=s^{\prime}-s$.
Thus, $\overline{\Gamma}_{1}(t-s)$ defined in Eq. (37) is a function of $t-s$
only, the same holds for $\overline{\Gamma}_{2}(t-s)$ defined in Eq. (41).
Using similar techniques, it can be also shown that
$\overline{\Gamma}_{0}(t-s)$ defined in Eq. (36), $\overline{\Gamma}_{A}(t-s)$
defined in Eq. (38) and $\overline{F}(\omega,t-s)$ are functions of $t-s$
only, which proves the functional dependencies assumed in Eq. (IV.1).

##  Appendix F Derivation of two-point correlation functions from Heisenberg
observables

We derive the two-point correlation function of the Heisenberg observable, to
reduce the notational complexity we here use a time-independent Liouville
operator ${\cal L}_{0}(\omega)$. Splitting the time propagation of the density
distribution in Eq. (10) into two steps we obtain

$\displaystyle{\rho}(\omega,t)$ $\displaystyle=e^{-(t-t^{\prime}){\cal
L}_{0}(\omega)}e^{-(t^{\prime}-t_{0}){\cal L}_{0}(\omega)}\rho(\omega,t_{0}).$
(112)

Introducing delta functions we obtain

$\displaystyle{\rho}(\omega,t)$ $\displaystyle=\int{\rm
d}\omega^{\prime}\int{\rm d}\omega_{0}e^{-(t-t^{\prime}){\cal
L}_{0}(\omega)}\delta(\omega-\omega^{\prime})e^{-(t^{\prime}-t_{0}){\cal
L}_{0}(\omega^{\prime})}\delta(\omega^{\prime}-\omega_{0})\rho(\omega_{0},t_{0})$
$\displaystyle $ $\displaystyle\equiv\int{\rm d}\omega^{\prime}\int{\rm
d}\omega_{0}{\rho}(\omega,t;\omega^{\prime},t^{\prime};\omega_{0},t_{0}),$
(113)

where in the second line we defined the three-point joint distribution
function ${\rho}(\omega,t;\omega^{\prime},t^{\prime};\omega_{0},t_{0})$. With
this distribution, the two-point correlation function defined in Eq. (57) can
be written as

$\displaystyle C(t,t^{\prime})$ $\displaystyle=\int{\rm d}\omega\int{\rm
d}\omega^{\prime}\int{\rm d}\omega_{0}\left(A_{S}(\omega)-\langle
A_{S}\rangle\right)\left(A_{S}(\omega^{\prime})-\langle
A_{S}\rangle\right){\rho}(\omega,t;\omega^{\prime},t^{\prime};\omega_{0},t_{0}).$
(114)

Inserting the definition of
${\rho}(\omega,t;\omega^{\prime},t^{\prime};\omega_{0},t_{0})$ from Eq. (F)
and using that ${\cal L}_{0}(\omega)$ is anti-self-adjoint, we obtain

$\displaystyle C(t,t^{\prime})$ $\displaystyle=\int{\rm
d}\omega{\rho}(\omega,t_{0})e^{(t^{\prime}-t_{0}){\cal
L}_{0}(\omega)}\left(A_{S}(\omega)-\langle
A_{S}\rangle\right)e^{(t-t^{\prime}){\cal
L}_{0}(\omega)}\left(A_{S}(\omega)-\langle A_{S}\rangle\right).$ (115)

We now use the product propagation relation proved in Appendix G and obtain

$\displaystyle C(t,t^{\prime})$ $\displaystyle=\int{\rm
d}\omega{\rho}(\omega,t_{0})\left(e^{(t^{\prime}-t_{0}){\cal
L}_{0}(\omega)}\left(A_{S}(\omega)-\langle
A_{S}\rangle\right)\right)e^{(t-t_{0}){\cal
L}_{0}(\omega)}\left(A_{S}(\omega)-\langle A_{S}\rangle\right),$ (116)

which can, using the definition of the Heisenberg observable, be rewritten as

$\displaystyle C(t,t^{\prime})$ $\displaystyle=\int{\rm
d}\omega{\rho}(\omega,t_{0})\left(A(\omega,t_{0},t)-\langle
A_{S}\rangle\right)\left(A(\omega,t_{0},t^{\prime})-\langle
A_{S}\rangle\right).$ (117)

In this expression $t_{0}$ denotes the projection time and $t\geq t_{0}$ and
$t^{\prime}\geq t_{0}$ denote arbitrary moments in time, note that we suppress
the dependence of the Heisenberg observable on the projection time in the rest
of the paper.

##  Appendix G Derivation of product propagation relation

To reduce the notational complexity, we here use a time-independent Liouville
operator ${\cal L}_{0}(\omega)$ and suppress the phase-space dependence of
${\cal L}_{0}(\omega)$ and of the two general phase-space function $A(\omega)$
and $B(\omega)$ in the following. The Liouville propagator acting on the
product of two phase-space functions follows from the series definition of the
exponential as

$\displaystyle e^{t{\cal L}_{0}}AB=\sum_{n=0}^{\infty}\frac{t^{n}{\cal
L}_{0}^{n}}{n!}AB.$ (118)

Since ${\cal L}_{0}$ is a linear differential operator the product rule
applies, which reads

$\displaystyle{\cal L}_{0}AB=A{\cal L}_{0}B+B{\cal L}_{0}A.$ (119)

Applying the product rule recursively, we find

$\displaystyle{\cal L}_{0}^{n}AB=\sum_{m=0}^{n}\frac{n!}{m!(n-m)!}\left({\cal
L}_{0}^{m}A\right){\cal L}_{0}^{n-m}B.$ (120)

Combining Eqs. (118) and (120) we find

$\displaystyle e^{t{\cal L}_{0}}AB$
$\displaystyle=\sum_{n=0}^{\infty}\sum_{m=0}^{n}\frac{t^{n}}{m!(n-m)!}\left({\cal
L}_{0}^{m}A\right){\cal L}_{0}^{n-m}B$ $\displaystyle $
$\displaystyle=\sum_{m=0}^{\infty}\sum_{n=m}^{\infty}\frac{t^{n}}{m!(n-m)!}\left({\cal
L}_{0}^{m}A\right){\cal L}_{0}^{n-m}B$ $\displaystyle $
$\displaystyle=\sum_{m=0}^{\infty}\sum_{n=0}^{\infty}\frac{t^{n+m}}{m!n!}\left({\cal
L}_{0}^{m}A\right){\cal L}_{0}^{n}B.$ (121)

Using Eq. (118) we finally arrive at the product popagation relation

$\displaystyle e^{t{\cal L}_{0}}AB=\left(e^{t{\cal L}_{0}}A\right)e^{t{\cal
L}_{0}}B.$ (122)

##  Appendix H Force averaging of two-point correlation function

In this section we derive the two-point correlation function Eq. (57) that
results from an average over phase space and over the stochastic non-
equilibrium force $h(t)$. We start from the non-equilibrium GLE Eq. (III.1).
To reduce the notational burden we consider an unconfined system that is
characterized by a diverging second moment, $\langle(A_{S}(\omega)-\langle
A_{S}\rangle)^{2}\rangle=\infty$, in which case $K(t)=0=\Gamma_{A}(s,t)$. Eq.
(III.1) in this scenario simplifies to

$\displaystyle
F(\omega,t_{0},t)+h(t)/M=\ddot{A}(\omega,t)+\int_{t_{0}}^{t}{\rm
d}s\,\Gamma(s,t)\dot{A}(\omega,s).$ (123)

Averaging the product of Eq. (123) at two different times $t$ and $t^{\prime}$
over phase space and $h(t)$ we obtain the expression

$\displaystyle\overline{\left\langle\left(F(\omega,t_{0},t)+ \frac{h(t)}{M}
\right)\left(F(\omega,t_{0},t^{\prime})+ \frac{h(t^{\prime})}{M}
\right)\right\rangle}=\overline{\left\langle
F(\omega,t_{0},t)F(\omega,t_{0},t^{\prime})\right\rangle}+\overline{
\frac{h(t)h(t^{\prime})}{M^{2}} }$ $\displaystyle $
$\displaystyle=\overline{\left\langle\left(\ddot{A}(\omega,t)+\int_{t_{0}}^{t}{\rm
d}s\,\Gamma(s,t)\dot{A}(\omega,s)\right)\left(\ddot{A}(\omega,t^{\prime})+\int_{t_{0}}^{t^{\prime}}{\rm
d}s^{\prime}\,\Gamma(s^{\prime},t^{\prime})\dot{A}(\omega,s^{\prime})\right)\right\rangle}$
$\displaystyle $
$\displaystyle=\left\langle\left(\overline{\ddot{A}(\omega,t)+\int_{t_{0}}^{t}{\rm
d}s\,\Gamma(s,t)\dot{A}(\omega,s)}
\right)\left(\overline{\ddot{A}(\omega,t^{\prime})+\int_{t_{0}}^{t^{\prime}}{\rm
d}s^{\prime}\,\Gamma(s^{\prime},t^{\prime})\dot{A}(\omega,s^{\prime})}
\right)\right\rangle+\Delta X,$ (124)

where in the first equation we used that $h(t)$ is independent of phase space
and that $\langle F(\omega,t_{0},t)\rangle=0$. We defined a correction term as

$\displaystyle\Delta X=$
$\displaystyle\overline{\left\langle\left(\ddot{A}(\omega,t)+\int_{t_{0}}^{t}{\rm
d}s\,\Gamma(s,t)\dot{A}(\omega,s)\right)\left(\ddot{A}(\omega,t^{\prime})+\int_{t_{0}}^{t^{\prime}}{\rm
d}s^{\prime}\,\Gamma(s^{\prime},t^{\prime})\dot{A}(\omega,s^{\prime})\right)\right\rangle}$
$\displaystyle $
$\displaystyle-\left\langle\left(\overline{\ddot{A}(\omega,t)+\int_{t_{0}}^{t}{\rm
d}s\,\Gamma(s,t)\dot{A}(\omega,s)}
\right)\left(\overline{\ddot{A}(\omega,t^{\prime})+\int_{t_{0}}^{t^{\prime}}{\rm
d}s^{\prime}\,\Gamma(s^{\prime},t^{\prime})\dot{A}(\omega,s^{\prime})}
\right)\right\rangle.$ (125)

Using that the observable $A(\omega,t)$ as the function that is determined by
the differential equation (123) does not explicitly depend on the non-
equilibrium force $h(t)$, the correction term can be rewritten as

$\displaystyle\Delta X=$ $\displaystyle\int_{t_{0}}^{t}{\rm
d}s\,\int_{t_{0}}^{t^{\prime}}{\rm
d}s^{\prime}\,\left\langle\dot{A}(\omega,s)\dot{A}(\omega,s^{\prime})\right\rangle\left(\overline{\Gamma(s,t)\Gamma(s^{\prime},t^{\prime})}-\overline{\Gamma(s,t)}
\overline{\Gamma(s^{\prime},t^{\prime})}\right)$ (126)

in terms of the second non-equilibrium force cumulant of the memory kernel.
Performing the average over $h(t)$ and neglecting the correction term $\Delta
X$ in Eq. (H) we obtain

$\displaystyle\overline{\left\langle
F(\omega,t-t_{0})F(\omega,t^{\prime}-t_{0})\right\rangle}+
\frac{\sigma(t-t^{\prime})}{M^{2}}
=\left\langle\left(\ddot{A}(\omega,t)+\int_{t_{0}}^{t}{\rm
d}s\,\overline{\Gamma}(t-s)\dot{A}(\omega,s)
\right)\left(\ddot{A}(\omega,t^{\prime})+\int_{t_{0}}^{t^{\prime}}{\rm
d}s^{\prime}\,\overline{\Gamma}(t^{\prime}-s^{\prime})\dot{A}(\omega,s^{\prime})
\right)\right\rangle.$ (127)

Performing a double Fourier integral over $t$ and $t^{\prime}$ and using Eqs.
(36) and (IV.1), we obtain

$\displaystyle 2\pi\delta(\nu-\nu^{\prime})\left[\langle({\cal
L}_{0}A_{S}(\omega))^{2}\rangle\tilde{\overline{\Gamma}}_{0}(\nu)+\tilde{\sigma}(\nu)/M^{2}\right]=\tilde{\chi}^{-1}(\nu)\tilde{\chi}^{-1}(\nu^{\prime})\left\langle\tilde{A}(\omega,\nu)\tilde{A}(\omega,\nu^{\prime})\right\rangle.$
(128)

Realizing that the double Fourier integral over $t$ and $t^{\prime}$ of the
two-point correlation function $C(t-t^{\prime})$ is given by

$\displaystyle
2\pi\delta(\nu-\nu^{\prime})\tilde{C}(\nu)=\left\langle\tilde{A}(\omega,\nu)\tilde{A}(\omega,\nu^{\prime})\right\rangle,$
(129)

Eq. (58) follows. Inspection of Eq. (58) shows that it is linear in the
complementary force and non-equilibrium force second moments
$\tilde{\overline{\Gamma}}_{0}(\nu)$ and $\tilde{\sigma}(\nu)$. The correction
term $\Delta X$ contains higher-order powers of these moments and is
neglected.

##  Appendix I Derivation of the standard FDT in the Heisenberg picture

Here we derive the standard (i.e. non-equilibrium) FDT by leading-order
perturbation analysis of the Heisenberg observable in Eq. (13), which deviates
from the text-book derivation of the standard FDT. To proceed, using the
expression for the Heisenberg observable Eq. (13) and the expression for the
time-dependent Liouville operator Eq. (18), we obtain

$\displaystyle A(\omega,t)=\exp_{H}\left(\int_{t_{0}}^{t}{\rm d}s{\cal
L}(s)\right)A_{S}(\omega)=\exp_{H}\left(\int_{t_{0}}^{t}{\rm d}s\left({\cal
L}_{0}-h(s)\Delta{\cal L}\right)\right)A_{S}(\omega),$ (130)

where the Heisenberg propagator is defined in Eq. (14). Expanding the operator
exponential to first order in $h(t)$, as shown in Eq. (E), we obtain

$\displaystyle A(\omega,t)=\left[ e^{(t-t_{0}){\cal L}_{0}}
-\int_{t_{0}}^{t}{\rm d}sh(s)e^{(s-t){\cal L}_{0}}\Delta{\cal L}e^{(t-s){\cal
L}_{0}}\right]A_{S}(\omega).$ (131)

The mean observable follows using the definition Eq. (28) as

$\displaystyle a(t)=\left\langle A(\omega,t)\right\rangle=\left\langle
e^{(t-t_{0}){\cal L}_{0}}A_{S}(\omega)\right\rangle
-\left\langle\int_{t_{0}}^{t}{\rm d}sh(s)e^{(s-t){\cal L}_{0}}\Delta{\cal
L}e^{(t-s){\cal L}_{0}}A_{S}(\omega)\right\rangle,$ (132)

which can be rewritten, using the anti-self-adjointedness of ${\cal L}_{0}$
and $\Delta{\cal L}$, as

$\displaystyle a(t)$ $\displaystyle=\left\langle A_{S}(\omega)\right\rangle
-\int_{t_{0}}^{t}{\rm d}s\beta h(s)\left\langle A_{S}(\omega){\cal
L}_{0}e^{(t-s){\cal L}_{0}}A_{S}(\omega)\right\rangle$
$\displaystyle=\left\langle A_{S}(\omega)\right\rangle -\int_{t_{0}}^{t}{\rm
d}s\beta h(s)\frac{{\rm d}}{{\rm d}t} \left\langle A_{S}(\omega)e^{(t-s){\cal
L}_{0}}A_{S}(\omega)\right\rangle$ $\displaystyle=\left\langle
A_{S}(\omega)\right\rangle +\int_{t_{0}}^{t}{\rm d}sh(s)\chi(t-s)/M.$ (133)

Here we defined the response function as

$\displaystyle\chi(t)=-\beta M\frac{{\rm d}}{{\rm d}t} \left\langle
A_{S}(\omega)e^{t{\cal L}_{0}}A_{S}(\omega)\right\rangle=-\beta M\frac{{\rm
d}}{{\rm d}t} C(t),$ (134)

which is the standard FDT as one finds it in text books, where the derivation
is typically done in the Schrödinger picture as there is no advantage of the
Heisenberg picture for this derivation. We in this paper use the Heisenberg
picture since it is needed when deriving the GLE. Note that Eq. (134) has no
dependence on the non-equilibrium force $h(t)$ whatsoever and thus is
equivalent to Eq. (61), which is the Fourier-transformed version of the FDT,
in the limit $h(t)=0$.

To make the equivalence between Eqs. (134) and (61) obvious, we Fourier
transform Eqs. (I) and (134). For this, we define the response function
$\chi(t)$ as single-sided and choose $t_{0}\rightarrow-\infty$, after which
Eq. (I) reads

$\displaystyle a(t)$ $\displaystyle=\left\langle A_{S}(\omega)\right\rangle
+\int_{-\infty}^{\infty}{\rm d}sh(s)\chi(t-s)/M.$ (135)

Fourier transformation yields

$\displaystyle\tilde{a}(\nu)=2\pi\delta(\nu)\langle
A_{S}\rangle+\tilde{\chi}(\nu)\tilde{h}(\nu)/M,$ (136)

which is equivalent to Eq. (53). Accounting for the single-sidedness of
$\chi(t)$, Eq. (134) can be written as

$\displaystyle\chi(t)=-\beta M\theta(t)\frac{{\rm d}}{{\rm d}t} C(t),$ (137)

where $\theta(t)$ denotes the Heavyside function. After Fourier transformation
Eq.(137) reads

$\displaystyle\tilde{\chi}(\nu)=\beta MC(0)-\imath\nu\beta
M\tilde{C}^{+}(\nu).$ (138)

The odd part of this equation reads

$\displaystyle\tilde{\chi}(\nu)-\tilde{\chi}(-\nu)=-\imath\nu\beta
M\tilde{C}(\nu),$ (139)

where we used that the Fourier transform of the symmetric correlation function
is $\tilde{C}(\nu)=\tilde{C}^{+}(\nu)+\tilde{C}^{+}(-\nu)$. Since the time-
domain response function $\chi(t)$ is a real function, we finally obtain

$\displaystyle 2\Im\left(\tilde{\chi}(\nu)\right)=-\nu\beta
M\tilde{C}(\nu)=-\frac{ \nu\tilde{C}(\nu)}{\langle({\cal
L}_{0}A_{S}(\omega))^{2}\rangle},$ (140)

where we used the definition of the mass M in Eq. (39) and which is equivalent
to Eq. (61), provided we take the equilibrium limit $h(t)=0$ in Eq. (61).

##  Appendix J Derivation of joint distribution function via path integrals

Inserting the Fourier-transform of Eq. (74) into Eq. (IV.5) we obtain

$\displaystyle\rho(A_{2},t_{2};A_{1},t_{1})=\int_{-\infty}^{\infty}\frac{{\rm
d}q_{1}}{2\pi} \int_{-\infty}^{\infty}\frac{{\rm d}q_{2}}{2\pi} e^{\imath
q_{1}(A_{1}-\langle A_{S}\rangle)+\imath q_{2}(A_{2}-\langle A_{S}\rangle)}
\overline{\left\langle\exp\left(-\imath\int_{-\infty}^{\infty}\frac{{\rm
d}\nu}{2\pi}\left(\tilde{F}(\omega,\nu)+\tilde{h}(\nu)/M\right)\tilde{g}(\nu)\right)\right\rangle}$
(141)

where we have defined

$\displaystyle\tilde{g}(\nu)\equiv\tilde{\chi}(\nu)\left(q_{1}e^{\imath\nu
t_{1}}+q_{2}e^{\imath\nu t_{2}}\right).$ (142)

We diagonalize the path integrals in Eqs. (IV.5) and (IV.5) by Fourier
transformation and obtain

$\displaystyle\langle
X(\tilde{F},\tilde{h})\rangle=\int_{-\infty}^{\infty}\frac{{\cal
D}\tilde{F}(\omega,\cdot)}{\tilde{\cal
N}_{F}}X(\tilde{F},\tilde{h})\exp\left(-\int_{-\infty}^{\infty}\frac{{\rm
d}\nu}{2\pi}
\frac{\tilde{F}(\omega,\nu)\tilde{F}(\omega,-\nu)}{2{\langle({\cal
L}_{0}A_{S}(\omega))^{2}\rangle}\tilde{\Gamma}_{0}(\nu)}\right),$ (143)
$\displaystyle\overline{X(\tilde{F},\tilde{h})}=\int_{-\infty}^{\infty}\frac{{\cal
D}\tilde{h}(\cdot)}{\tilde{\cal
N}_{h}}X(\tilde{F},\tilde{h})\exp\left(-\int_{-\infty}^{\infty}\frac{{\rm
d}\nu}{2\pi}
\frac{\tilde{h}(\nu)\tilde{h}(-\nu)}{2\tilde{\sigma}(\nu)}\right),$ (144)

where $\tilde{\cal N}_{F}$ and $\tilde{\cal N}_{h}$ are normalization
constants. Now the averages over $\tilde{F}(\omega,\nu)$ and $\tilde{h}(\nu)$
in Eq. (141) can be done by explicitly performing the path integrals, after
which we obtain

$\displaystyle\rho(A_{2},t_{2};A_{1},t_{1})$
$\displaystyle=\int_{-\infty}^{\infty}\frac{{\rm d}q_{1}}{2\pi}
\int_{-\infty}^{\infty}\frac{{\rm d}q_{2}}{2\pi} e^{\imath q_{1}(A_{1}-\langle
A_{S}\rangle)+\imath q_{2}(A_{2}-\langle A_{S}\rangle)} $
$\displaystyle\times\exp\left(-\int_{-\infty}^{\infty}\frac{{\rm
d}\nu}{4\pi}\tilde{g}(\nu)\tilde{g}(-\nu)\left(\langle({\cal
L}_{0}A_{S}(\omega))^{2}\rangle\tilde{\Gamma}_{0}(\nu)+\tilde{\sigma}(\nu)/M^{2}\right)\right).$
(145)

Inserting the definition of $\tilde{g}(\nu)$ from Eq. (142) we obtain

$\displaystyle\rho(A_{2},t_{2};A_{1},t_{1})$
$\displaystyle=\int_{-\infty}^{\infty}\frac{{\rm d}q_{1}}{2\pi}
\int_{-\infty}^{\infty}\frac{{\rm d}q_{2}}{2\pi} e^{\imath q_{1}(A_{1}-\langle
A_{S}\rangle)+\imath q_{2}(A_{2}-\langle A_{S}\rangle)}
\exp\left(-\int_{-\infty}^{\infty}\frac{{\rm
d}\nu}{4\pi}\tilde{C}(\nu)\sum_{j,k=1}^{2}q_{j}q_{k}e^{\imath\nu(t_{j}-t_{k})}\right),$
(146)

where we used the explicit expression for the Fourier-transformed correlation
function $\tilde{C}(\nu)$ from Eq. (58). Now the Fourier integral in the
exponent can be done and we obtain

$\displaystyle\rho(A_{2},t_{2};A_{1},t_{1})$
$\displaystyle=\int_{-\infty}^{\infty}\frac{{\rm d}q_{1}}{2\pi}
\int_{-\infty}^{\infty}\frac{{\rm d}q_{2}}{2\pi} e^{\imath q_{1}(A_{1}-\langle
A_{S}\rangle)+\imath q_{2}(A_{2}-\langle A_{S}\rangle)}
\exp\left(-\frac{1}{2}\sum_{j,k=1}^{2}q_{j}q_{k}C(t_{k}-t_{j})\right),$ (147)

which can be slightly rewritten as

$\displaystyle\rho(A_{2},t_{2};A_{1},t_{1})$
$\displaystyle=\int_{-\infty}^{\infty}\frac{{\rm d}q_{1}}{2\pi}
\int_{-\infty}^{\infty}\frac{{\rm d}q_{2}}{2\pi}
\exp\left(-\frac{1}{2}\sum_{j,k=1}^{2}q_{j}q_{k}C(t_{k}-t_{j})+\imath\sum_{j=1}^{2}q_{j}(A_{j}-\langle
A_{S}\rangle)\right).$ (148)

Now we use the definition from the main text Eq. (78),
$I_{jk}=C(t_{j}-t_{k})$, and obtain

$\displaystyle\rho(A_{2},t_{2};A_{1},t_{1})$
$\displaystyle=\int_{-\infty}^{\infty}\frac{{\rm d}q_{1}}{2\pi}
\int_{-\infty}^{\infty}\frac{{\rm d}q_{2}}{2\pi}
\exp\left(-\frac{1}{2}\sum_{j,k=1}^{2}q_{j}q_{k}I_{jk}+\imath\sum_{j=1}^{2}q_{j}(A_{j}-\langle
A_{S}\rangle)\right).$ (149)

The Gaussian integrals over $q_{1}$ and $q_{2}$ can be done, after which we
obtain the final result reported in Eq. (IV.5),

$\displaystyle\rho(A_{2},t_{2};A_{1},t_{1})=\frac{\exp\left(-\sum_{j,k=1}^{2}(A_{j}-\langle
A_{S}\rangle)I_{jk}^{-1}(A_{k}-\langle A_{S}\rangle)/2\right)}{\sqrt{\det 2\pi
I}},$ (150)

which is the two-point distribution of a general Gaussian process. Note that
in Eq. (IV.5) the Einstein summation convention is used.

As a final step, we show how to marginalize the two-point distribution. This
is most easily done on the level of Eq. (149) where the positions $A_{1}$ and
$A_{2}$ appear linearly. By integration over one observable value we obtain

$\displaystyle\rho(A_{1},t_{1})$ $\displaystyle=\int_{-\infty}^{\infty}{\rm
d}A_{2}\rho(A_{2},t_{2};A_{1},t_{1})=\int_{-\infty}^{\infty}\frac{{\rm
d}q_{1}}{2\pi} \exp\left(-\frac{q_{1}^{2}C(0)}{2}+q_{1}(A_{1}-\langle
A_{S}\rangle)\right)=\frac{\exp\left(-\frac{(A_{1}-\langle
A_{S}\rangle)^{2}}{2C(0)}\right)}{\sqrt{2\pi C(0)}},$ (151)

which, expectedly, is a Gaussian distribution with a variance corresponding to
the equal-time correlation $C(0)$.

## References

  * Prigogine [1947] I. Prigogine,  _Etude thermodynamic des phenomene irreversibles_ (Desoer, Liege, 1947).
  * Prigogine and Mazur [1953] I. Prigogine and P. Mazur, Sur l’extension de la thermodynamique aux phenomenes irreversibles lies aux degres de liberte internes, Physica XIX, 241 (1953).
  * Lebowitz [1959] J. L. Lebowitz, Stationary nonequilibrium Gibbsian ensembles, Phys. Rev. 114, 1192 (1959).
  * Zwanzig [1960] R. Zwanzig, Ensemble method in the theory of irreversibility, The Journal of Chemical Physics 33, 1338 (1960).
  * de Groot and Mazur [1962] S. R. de Groot and P. Mazur,  _Non-Equilibrium Thermodynamics_ (North-Holland Pub. Co., Amsterdam, 1962).
  * Grabert  _et al._ [1980] H. Grabert, P. Hänggi, and P. Talkner, Microdynamics and nonlinear stochastic processes of gross variables, Journal of Statistical Physics 22, 537 (1980).
  * Risken [1984] H. Risken,  _The Fokker-Planck Equation_ (Springer, Berlin, 1984).
  * Zwanzig [2001] R. Zwanzig,  _Nonequilibrium statistical mechanics_ (Oxford UnivPress, Oxford [u.a.], 2001).
  * Mizuno  _et al._ [2007] D. Mizuno, C. Tardin, C. F. Schmidt, and F. C. MacKintosh, Nonequilibrium mechanics of active cytoskeletal networks, Science 315, 370 (2007).
  * Gomez-Solano  _et al._ [2009] J. R. Gomez-Solano, A. Petrosyan, S. Ciliberto, R. Chetrite, and K. Gawedzki, Experimental verification of a modified fluctuation-dissipation relation for a micron-sized particle in a nonequilibrium steady state, Phys. Rev. Lett. 103, 040601 (2009).
  * Mehl  _et al._ [2010] J. Mehl, V. Blickle, U. Seifert, and C. Bechinger, Experimental accessibility of generalized fluctuation-dissipation relations for nonequilibrium steady states, Phys. Rev. E 82, 032401 (2010).
  * Theurkauff  _et al._ [2012] I. Theurkauff, C. Cottin-Bizonne, J. Palacci, C. Ybert, and L. Bocquet, Dynamic clustering in active colloidal suspensions with chemical signaling, Phys. Rev. Lett. 108, 268303 (2012).
  * Dinis  _et al._ [2012] L. Dinis, P. Martin, J. Barral, J. Prost, and J.-F. Joanny, Fluctuation-response theorem for the active noisy oscillator of the hair-cell bundle, Phys. Rev. Lett. 109, 160602 (2012).
  * Bohec  _et al._ [2013] P. Bohec, F. F. Gallet, C. Maes, S. Safaverdi, P. Visco, and F. van Wijland, Probing active forces via a fluctuation-dissipation relation: Application to living cells, Europhys. Lett. 102, 50005 (2013).
  * Guo  _et al._ [2014] M. Guo, A. J. Ehrlicher, M. H. Jensen, M. Renz, J. R. Moore, R. D. Goldman, J. Lippincott-Schwartz, F. C. MacKintosh, and D. A. Weitz, Cell 158, 822 (2014).
  * Krug [1991] J. Krug, Boundary-induced phase transitions in driven diffusive systems, Phys. Rev. Lett. 67, 1882 (1991).
  * Schmittmann and Zia [1998] B. Schmittmann and R. K. P. Zia, Driven diffusive systems. an introduction and recent developments, Phys. Rep. 301, 45 (1998).
  * Derrida  _et al._ [2001] B. Derrida, J. L. Lebowitz, and E. R. Speer, Free energy functional for nonequilibrium systems: An exactly solvable case, Phys. Rev. Lett. 87, 150601 (2001).
  * Ilg and Barrat [2007] P. Ilg and J. L. Barrat, From single-particle to collective effective temperatures in an active fluid of self-propelled particles, Europhys. Lett. 111, 26001 (2007).
  * Grosberg and Joanny [2015] A. Y. Grosberg and J. F. Joanny, Nonequilibrium statistical mechanics of mixtures of particles in contact with different thermostats, Phys. Rev. E 92, 032118 (2015).
  * Fodor  _et al._ [2016] E. Fodor, C. Nardini, M. E. Cates, J. Tailleur, P. Visco, and F. van Wijland, How far from equilibrium is active matter?, Phys. Rev. Lett. 117, 038103 (2016).
  * Carlon  _et al._ [2018] E. Carlon, H. Orland, T. Sakaue, and C. Vanderzande, Effect of memory and active forces on transition path time distributions, J. Phys. Chem. B 122, 11186 (2018).
  * Lavacchi  _et al._ [2022] L. Lavacchi, J. O. Daldrop, and R. R. Netz, Non-arrhenius barrier crossing dynamics of non-equilibrium non-markovian systems, Europhys. Lett. 139, 51001 (2022).
  * Jarzynski [2000] C. Jarzynski, Hamiltonian derivation of a detailed fluctuation theorem, Journal of Statistical Physics 98, 77 (2000).
  * Hatano and Sasa [2001] T. Hatano and S. I. Sasa, Steady-state thermodynamics of langevin systems, Phys. Rev. Lett. 86, 3463 (2001).
  * Harada and Sasa [2005] T. Harada and S. I. Sasa, Equality connecting energy dissipation with a violation of the fluctuation-response relation, Phys. Rev. Lett. 95, 130602 (2005).
  * Seifert [2005] U. Seifert, Entropy production along a stochastic trajectory and an integral fluctuation theorem, Phys. Rev. Lett. 95, 040602 (2005).
  * Lapolla and Godec [2020] A. Lapolla and A. Godec, Faster uphill relaxation in thermodynamically equidistant temperature quenches, Phys. Rev. Lett. 125, 110602 (2020).
  * Prost  _et al._ [2009] J. Prost, J.-F. Joanny, and J. M. R. Parrondo, Generalized fluctuation-dissipation theorem for steady-state systems, Phys. Rev. Lett. 103, 090601 (2009).
  * Baiesi  _et al._ [2009] M. Baiesi, C. Maes, and B. Wynants, Fluctuations and response of nonequilibrium states, Phys. Rev. Lett. 103, 010602 (2009).
  * Seifert and Speck [2010] U. Seifert and T. Speck, Fluctuation-dissipation theorem in nonequilibrium steady states, Europhys. Lett. 89, 10007 (2010).
  * Willareth  _et al._ [2017] L. Willareth, I. M. Sokolov, Y. Roichman, and B. Lindner, Generalized fluctuation-dissipation theorem as a test of the markovianity of a system, Europhys. Lett. 118, 20001 (2017).
  * Netz [2018] R. R. Netz, Fluctuation-dissipation relation and stationary distribution of an exactly solvable many-particle model for active biomatter far from equilibrium, J. Chem. Phys. 148, 185101 (2018).
  * Netz [2020] R. R. Netz, Approach to equilibrium and nonequilibrium stationary distributions of interacting many-particle systems that are coupled to different heat baths, Phys. Rev. E 101, 022120 (2020).
  * Nakajima [1958] S. Nakajima, On Quantum Theory of Transport PhenomenaSteady Diffusion, Progress of Theoretical Physics 20, 948 (1958).
  * Zwanzig [1961] R. Zwanzig, Memory Effects in Irreversible Thermodynamics, Physical Review 124, 983 (1961).
  * Mori [1965] H. Mori, Transport, Collective Motion, and Brownian Motion, Progress of Theoretical Physics 33, 423 (1965).
  * Ciccotti and Ryckaert [1981] G. Ciccotti and J.-P. Ryckaert, On the derivation of the gen- eralized Langevin equation for interacting Brownian particles, Journal of Statistical Physics 26, 73 (1981).
  * Straub  _et al._ [1987] J. E. Straub, M. Borkovec, and B. J. Berne, Calculation of dynamic friction on intramolecular degrees of freedom, The Journal of Physical Chemistry 91, 4995 (1987).
  * Lange and Grubmüller [2006] O. F. Lange and H. Grubmüller, Collective Langevin dynamics of conformational motions in proteins, The Journal of Chemical Physics 124, 214903 (2006), publisher: American Institute of Physics.
  * Kinjo and Hyodo [2007] T. Kinjo and S. A. Hyodo, Equation of motion for coarse-grained simulation based on microscopic description, Physical Review E 75, 051109 (2007), publisher: American Physical Society.
  * Darve  _et al._ [2009] E. Darve, J. Solomon, and A. Kia, Computing generalized Langevin equations and generalized Fokker–Planck equations, Proceedings of the National Academy of Sciences 106, 10884 (2009).
  * Hijón  _et al._ [2010] C. Hijón, P. Español, E. Vanden-Eijnden, and R. Delgado-Buscalioni, Mori–Zwanzig formalism as a practical computational tool, Faraday Discussions 144, 301 (2010), publisher: Royal Society of Chemistry.
  * Izvekov [2013] S. Izvekov, Microscopic derivation of particle-based coarse-grained dynamics, The Journal of Chemical Physics 138, 134106 (2013), publisher: American Institute of Physics.
  * Lee  _et al._ [2019] H. S. Lee, S.-H. Ahn, and E. F. Darve, The multi-dimensional generalized langevin equation for conformational motion of proteins, The Journal of chemical physics 150, 174113 (2019).
  * Ayaz  _et al._ [2022] C. Ayaz, L. Scalfi, B. A. Dalton, and R. R. Netz, Generalized langevin equation with a nonlinear potential of mean force and nonlinear memory friction from a hybrid projection scheme, Physical Review E 105, 054138 (2022).
  * Vroylandt  _et al._ [2022] H. Vroylandt, L. Goudenège, P. Monmarché, F. Pietrucci, and B. Rotenberg, Likelihood-based non-markovian models from molecular dynamics, Proceedings of the National Academy of Sciences 119, e2117586119 (2022).
  * Plotkin and Wolynes [1998] S. S. Plotkin and P. G. Wolynes, Non-Markovian Configurational Diffusion and Reaction Coordinates for Protein Folding, Physical Review Letters 80, 5015 (1998), publisher: American Physical Society.
  * Satija and Makarov [2019] R. Satija and D. E. Makarov, Generalized Langevin Equation as a Model for Barrier Crossing Dynamics in Biomolecular Folding, The Journal of Physical Chemistry B 123, 802 (2019), publisher: American Chemical Society.
  * Ayaz  _et al._ [2021] C. Ayaz, L. Tepper, F. N. Brünig, J. Kappler, J. O. Daldrop, and R. R. Netz, Non-Markovian modeling of protein folding, Proceedings of the National Academy of Sciences 118, 10.1073/pnas.2023856118 (2021), publisher: National Academy of Sciences Section: Physical Sciences.
  * Dalton  _et al._ [2023] B. A. Dalton, C. Ayaz, H. Kiefer, A. Klimek, L. Tepper, and R. R. Netz, Fast protein folding is governed by memory-dependent friction, Proceedings of the National Academy of Sciences 120, e2220068120 (2023).
  * Bagchi and Oxtoby [1983] B. Bagchi and D. W. Oxtoby, The effect of frequency dependent friction on isomerization dynamics in solution, The Journal of Chemical Physics 78, 2735 (1983).
  * Straub  _et al._ [1986] J. E. Straub, M. Borkovec, and B. J. Berne, Non‐Markovian activated rate processes: Comparison of current theories with numerical simulation data, The Journal of Chemical Physics 84, 1788 (1986).
  * Pollak  _et al._ [1989] E. Pollak, H. Grabert, and P. Hänggi, Theory of activated rate processes for arbitrary frequency dependent friction: Solution of the turnover problem, The Journal of Chemical Physics 91, 4073 (1989).
  * Brünig  _et al._ [2022a] F. N. Brünig, J. O. Daldrop, and R. R. Netz, Pair-Reaction Dynamics in Water: Competition of Memory, Potential Shape, and Inertial Effects, The Journal of Physical Chemistry B 126, 10295 (2022a).
  * Mitterwallner  _et al._ [2020] B. G. Mitterwallner, C. Schreiber, J. O. Daldrop, J. O. Rädler, and R. R. Netz, Non-Markovian data-driven modeling of single-cell motility, Physical Review E 101, 032408 (2020).
  * Tuckerman and Berne [1993] M. Tuckerman and B. Berne, Vibrational relaxation in simple fluids: Comparison of theory and simulation, The Journal of Chemical Physics 98, 7301 (1993).
  * Gottwald  _et al._ [2015] F. Gottwald, S. D. Ivanov, and O. Kühn, Applicability of the Caldeira–Leggett model to vibrational spectroscopy in solution, The Journal of Physical Chemistry Letters 6, 2722 (2015).
  * Brünig  _et al._ [2022b] F. N. Brünig, O. Geburtig, A. von Canal, J. Kappler, and R. R. Netz, Time-Dependent Friction Effects on Vibrational Infrared Frequencies and Line Shapes of Liquid Water, The Journal of Physical Chemistry B 126, 1579 (2022b).
  * Herrera-Delgado  _et al._ [2020] E. Herrera-Delgado, J. Briscoe, and P. Sollich, Tractable nonlinear memory functions as a tool to capture and explain dynamical behaviors, Physical Review Research 2, 043069 (2020).
  * Chorin  _et al._ [2000] A. J. Chorin, O. H. Hald, and R. Kupferman, Optimal prediction and the Mori–Zwanzig representation of irreversible processes, Proceedings of the National Academy of Sciences 97, 2968 (2000), publisher: National Academy of Sciences Section: Physical Sciences.
  * Robertson [1966] B. Robertson, Equations of motion in nonequilibriurn statistical mechanics, Physical Review 144, 151 (1966).
  * Nordholm and Zwanzig [1975] S. Nordholm and R. Zwanzig, A systematic derivation of exact generalized brownian motion theory, Journal of Statistical Physics 13, 347 (1975).
  * Picard and Willis [1977] R. H. Picard and C. R. Willis, Time-dependent projection-operator approach to master equations for coupled systems. ii. systems with correlations, Physical Review A 16, 1625 (1977).
  * Uchiyama and Shibata [1999] C. Uchiyama and F. Shibata, Unified projection operator formalism in nonequilibrium statistical mechanics, Physical Review E 60, 2636 (1999).
  * Koide [2002] T. Koide, Derivation of transport equations using the time-dependent projection operator method, Progress of Theoretical Physics 60, 525 (2002).
  * Latz [2002] A. Latz, Non-equilibrium projection-operator for a quenched thermostatted system, Journal of Statistical Physics 109, 607 (2002).
  * Meyer  _et al._ [2017] H. Meyer, T. Voigtmann, and T. Schilling, On the non-stationary generalized Langevin equation, The Journal of Chemical Physics 147, 214110 (2017), publisher: American Institute of Physics.
  * Cui and Zaccone [2018] B. Cui and A. Zaccone, Generalized langevin equation and fluctuation-dissipation theorem for particle-bath systems in external oscillating fields, Physical Review E 97, 060102(R) (2018).
  * te Vrugt and Wittkowski [2019] M. te Vrugt and R. Wittkowski, Mori-zwanzig projection operator formalism for far-from-equilibrium systems with time-dependent hamiltonians, Physical Review E 99, 062118 (2019).
  * Meyer  _et al._ [2020] H. Meyer, P. Pelagejcev, and T. Schilling, Non-Markovian out-of-equilibrium dynamics: A general numerical procedure to construct time-dependent memory kernels for coarse-grained observables, EPL (Europhysics Letters) 128, 40001 (2020), publisher: IOP Publishing.
  * Dyson [1949] F. J. Dyson, The Radiation Theories of Tomonaga, Schwinger, and Feynman, Physical Review 75, 486 (1949), publisher: American Physical Society.
  * Feynman [1951] R. P. Feynman, An Operator Calculus Having Applications in Quantum Electrodynamics, Physical Review 84, 108 (1951), publisher: American Physical Society.
  * Carof  _et al._ [2014] A. Carof, R. Vuilleumier, and B. Rotenberg, Two algorithms to compute projected correlation functions in molecular dynamics simulations, The Journal of Chemical Physics 140, 124103 (2014), publisher: American Institute of Physics.
  * Jung  _et al._ [2017] G. Jung, M. Hanke, and F. Schmid, Iterative Reconstruction of Memory Kernels, Journal of Chemical Theory and Computation 13, 2481 (2017), publisher: American Chemical Society.
  * Daldrop  _et al._ [2017] J. O. Daldrop, B. G. Kowalik, and R. R. Netz, External Potential Modifies Friction of Molecular Solutes in Water, Physical Review X 7, 041065 (2017), publisher: American Physical Society.
  * Daldrop  _et al._ [2018] J. O. Daldrop, J. Kappler, F. N. Brünig, and R. R. Netz, Butane dihedral angle dynamics in water is dominated by internal friction, Proceedings of the National Academy of Sciences 115, 5169 (2018), publisher: National Academy of Sciences Section: Biological Sciences.
  * Klippenstein and van der Vegt [2021] V. Klippenstein and N. F. A. van der Vegt, Cross-correlation corrected friction in (generalized) Langevin models, The Journal of Chemical Physics 154, 191102 (2021), publisher: American Institute of Physics.