††thanks: Corresponding Author.

# Generalized optimal protocols of Brownian motion in a parabolic potential

Pedro J. Colmenares<EMAIL_ADDRESS>Departamento de Química. Universidad
de Los Andes. Mérida 5101, Venezuela

###### Abstract

The generalized Langevin equation with an exponential kernel is used to
analyze memory effects on the optimal work done by a Brownian particle in a
heat bath and subjected to a harmonic moving potential. The generalized
overdamping scenario is also investigated. Several facts emerge in these more
precise descriptions using the same initial conditions of the Markovian which
lead the particle to do mechanical work against the field. Compared with the
results obtained with the latter, the memory fades the discontinuities
observed in the highly underdamped regime, which suggests that this trait is a
consequence of the Markov approximation as well as the dependence of the
different dynamical susceptibilities with the external field. Unlike the
overdamped Markovian, work is done by the external field in the analog
generalized counterpart. A detailed calculation of the rate of entropy
production gives negative values. It is mathematically correct because the
dynamics deal with a reduced description of the degrees of freedom of the
bath. The theory then requires improving the treatment of them to restore the
second law and thus to get the results with the required thermodynamics
consistency.

Brownian motion, Stochastic processes.

###### pacs:

05.30.?d; 05.40.Jc

##  I Introduction

The thermodynamic work involving the interaction of a Brownian particle with a
heat bath and an external moving potential has been examined from different
perspectives. From experimental setups where the field was provided by an
optical trap displaced at constant speed Imparato  _et al._ (2007); Ciliberto
(2017) up to theoretical descriptions where the external protocol minimizes
the mechanical work produced by the interaction of the particle with the
external field Schmiedl and Seifert (2007); Abreu and Seifert (2011);
Colmenares and Paredes-Altuve (2021) have been the mayor researches in the
field. In the latter, the analysis has been focused on the Markovian Langevin
equation (MLE) in both underdamped and overdamped regimens. However, the
description of non-Markovian systems through a more exact description, such as
the classical generalized Langevin equation (GLE), has received less
attention.

Although the GLE has been successful in describing non-Markovian systems,
there are new results that need to be reconsidered in its formulation. For
instance, it can be highlight the finding by Daldrop et al. Daldrop  _et al._
(2017) in the molecular dynamics simulation of a harmonic particle where the
friction is not a constant but depends on the strength of the applied field.
It must be mentioned that Vroylandt et al. Vroylandt and Monmarché (2022) have
also shown that it depends on the position if a non-linear mean force is
considered. Then, the general consensus is that the field strength must be
included in its interaction with the degrees of freedom of the reservoir Lisý
and Tóthová (2019); Di Cairano (2022). In view of these findings, there are
even sufficient reasons to consider that Kubo’s fluctuation dissipation
theorem (FDT) associated with the GLE is not so accurate Lisý and Tóthová
(2019); Di Cairano (2022). This deficiency has often been overlooked in the
literature, even though it was pointed out decades ago by Balakrishnam
Balakrishnan (1979) and by Cicotti et al Cicotti and Ryckaert (1981). In fact,
the work by Olivares-Rivas et al. Olivares and Colmenares (2016) showed that
for a constant field the noise correlation depends on its squared strength.
Costa et al. Costa  _et al._ (2003) showed that it is even inconsistent in the
linear regime of superdiffusive processes. Additionally, it was also discussed
that a dynamics driven by nonequilibrium steady states, the resulting FDT can
be written as the equilibrium one and an additive correction Seifert and Speck
(2010). It has been also considered adding correlations between the system and
the environment Kutvonen  _et al._ (2015). Despite these basic observations
and potentially the derivation of others that are unknown at the moment, there
has been some progress using the original GLE in the analysis of the subject
matter of this research.

To date, the published works involving the moving harmonic potential
incorporate only the memory effect of the stochastic force due to the
reservoir through an exponential kernel as in Ref. Daldrop  _et al._ (2017).
Thus, Paredes-Altuve et al. Paredes-Altuve  _et al._ (2016) investigated how
the individual measurement of the average position or speed affects the
optimal protocol and therefore the minimum mechanical work. The no measurement
analysis was also included. There, it was shown that the inclusion of memory
due to the noise generated by the thermal bath effectively induces the system
to perform work against the field when the measurement is done. Furthermore,
Di Terlizzi et al. Di Terlizzi  _et al._ (2020) calculated the non-optimal
work and the variance of its distribution for optical tweezers displaced at
constant speed both in the GLE and in the generalized overdamped version
(GOLE).

Under these circumstances, the main objective of this research is to use the
classical GLE to analyze memory effects on the determination of the optimal
protocol and minimum work independently of any measurement at all. In
particular, to compare the new findings with those obtained from the Markovian
case Colmenares and Paredes-Altuve (2021) under the same initial conditions.
To do so, the variational procedure of Ref. Paredes-Altuve  _et al._ (2016)
must be replaced by the more general one implemented in Ref. Colmenares and
Paredes-Altuve (2021). It will be of particular interest to find a procedure
that, for a given friction coefficient, a reasonable value of the memory
kernel parameter be determined that leads the particle to perform mechanical
work against the external field. It will be done for any value of the friction
constant assuming that the particle starts from equilibrium and also from an
arbitrary initial position.

To achieve this goal, it will be shown in Sec. II the general equations
defining the system dynamics and mechanical work. The latter appears to be as
a complex functional of the protocol with local and non-local parts Colmenares
and Paredes-Altuve (2021). The solution for the GLE and GOLE is shown in two
separate subsections where the work functional is explicitly broken down into
its two contributions. For each one, the associated Euler-Lagrange equation is
derived and their sum Ferialdi and Bassi (2012) is solved to get the optimal
driving. It will be shown that for the two regimes, each protocol captures the
memory effect of the thermal bath though a well-defined convoluted integral
equation. They together, Eqs. (18) and (41), are the main results of this
proposal. Next, in Sec. III is presented the procedure to assign meaningful
physical value of the kernel decay constant and the discussion of the results
for the optimal driving and the minimum work. Additionally, the entropy
production associated to the dynamics is determined and discussed as ruled by
the second law. The article concludes with some general remarks in Sec. IV.

##  II General equations

The exact position $q(t)$ of a Brownian particle of mass $M$ submerged in a
bath at temperature $T$ and subjected to a moving harmonic potential with
stiffness $\kappa^{2}$ and external driving $\lambda(t)$, that is,
$V(q,t)=\kappa^{2}(q-\lambda(t))^{2}/2$, is provided by the GLE Zwanzig
(2001); Banik  _et al._ (2002):

$M\,\ddot{q}(t)=-\\!\\!\int_{0}^{t}\\!\\!dy\,\Gamma(t-y)\,q^{\prime}(y)-\kappa^{2}\left(q(t)\\!-\\!\lambda(t)\right)+\xi(t),$
(1)

where the dot and the apostrophe above a function denotes its time and normal
derivative, respectively.

The kernel $\Gamma(t-y)$ compiles the memory or retardation effect on the
movement of the particle due to the collective hydrodynamic response of the
bath and $\xi(t)$ is a zero-mean Gaussian colored noise with correlation
$\left<\xi(t)\xi(y)\right>=k_{{}_{\mathrm{B}}}T\,\Gamma(|t-y|)$ satisfying
Kubo’s FDT Reichl (1998) where $k_{{}_{\mathrm{B}}}$ is the Boltzmann
constant.

The analysis presented here will be framed under the guideline of the original
Kubo theorem with exponential kernel $\Gamma(t)=\tau\,\exp[-\alpha\,t]$ in
order to compare and complement the findings with those of Refs. Paredes-
Altuve  _et al._ (2016) and Di Terlizzi  _et al._ (2020), respectively.
Constant $\tau=\alpha\,\gamma$ being $1/\alpha$ the kernel decay relaxation
time and $\gamma$ the static friction coefficient of the thermal fluid.

Scaling time and length by the factors  $\kappa\,M^{-1/2}$ and
$\kappa(k_{{}_{\mathrm{B}}}T)^{-1/2}$, respectively, the GLE reduces to the
dimensionless equation

$\displaystyle\ddot{q}(t)$ $\displaystyle=$
$\displaystyle-\tau\\!\\!\int_{0}^{t}\\!\\!dy\,f(t-y)\,q^{\prime}(y)-q(t)+\lambda(t)+\xi(t),$
(2)

with $f(t-y)=\exp[-\alpha\,|t-y|]$. Without loss of generality, this is
equivalent to assume $\\{k_{{}_{\mathrm{B}}},\kappa,M\\}=1$.

The thermodynamic work $W$ has been defined according stochastic energetics in
Refs. Sekimoto (2010); Seifert (2012). It is given as a functional of the
protocol and its derivative, i.e.

$W[\lambda,\dot{\lambda}]=\int_{0}^{t_{{}_{\mathrm{f}}}}dt\,\dot{\lambda}(t)\left[\lambda(t)-\left<q[\lambda(t)]\right>\right],$
(3)

where $\left<q[\lambda(t)]\right>$ is the average over the noise distribution
and where it is indicated its explicit functional dependence with the
protocol.

The task is to find by means of variational techniques, the optimal
$\lambda(t)$ which makes $W[\lambda,\dot{\lambda}]$ to reach a minimum taking
into account that the functional has a local (“loc”) and non-local (“nl”)
contributions. Once they are identified, the first variation is

$\delta
W_{r}[\lambda,\dot{\lambda}]\\!=\\!\\!\int_{0}^{t_{{}_{\mathrm{f}}}}\\!\\!\\!\\!\\!dt\\!\left(\\!\frac{\partial
L_{\mathrm{r}}[\lambda,\dot{\lambda},t]}{\partial\lambda}-\\!\frac{d}{dt}\frac{\partial
L_{\mathrm{r}}[\lambda,\dot{\lambda},t]}{\partial\dot{\lambda}}\\!\right)\\!\delta\lambda(t),$
(4)

where $\mathrm{r}=\\{\mathrm{loc},\mathrm{nl}\\}$ and
$L_{\mathrm{r}}[\lambda,\dot{\lambda},t]$ is a definite function to be
determined. The first variation $\delta W$ is then given as the sum of the two
associated Euler-Lagrange equations Ferialdi and Bassi (2012); Colmenares and
Paredes-Altuve (2021) which, after setting it equal to zero will allow to
determine the optimal protocol $\lambda(t)$.

Once the protocol and the average position are derived then $W$ can be
integrated between $\\{0,t_{{}_{\mathrm{f}}}\\}$ where $t_{{}_{\mathrm{f}}}$
is the final time of the protocol application. It renders,

$W=\frac{1}{2}\lambda_{{}_{\mathrm{f}}}^{2}-\lambda_{{}_{\mathrm{f}}}\left<q(t_{{}_{\mathrm{f}}})\right>+\int_{0}^{t_{{}_{\mathrm{f}}}}dy\,\lambda(y)\left<q^{\prime}[\lambda(y)]\right>),$
(5)

where $\lambda_{{}_{\mathrm{f}}}=\lambda(t_{{}_{\mathrm{f}}})$ is a pre-fixed
values fixed by the external agent. We have supposed as in Ref. Colmenares and
Paredes-Altuve (2021) that $\lambda(0)=0$. The first two terms are just the
adiabatic work by changing the protocol to $\lambda_{{}_{\mathrm{f}}}$ and
instantly returning it to the initial value.

The determination of $q(t)$ for the GLE and for the generalized overdamped and
also their associated optimal protocol which makes the work minimal, will be
shown in the next two subsections.

###  II.1 GLE dynamics

The Laplace transform of Eq. (2) gives

$\displaystyle\widehat{q}(s)$ $\displaystyle=$ $\displaystyle
q_{{}_{\mathrm{0}}}\widehat{\chi_{{}_{\mathrm{q}}}\\!}(s)+v_{{}_{\mathrm{0}}}\,\widehat{\chi_{{}_{\mathrm{v}}}\\!}(s)+\widehat{\chi_{{}_{\mathrm{v}}}\\!}(s)\,\widehat{\lambda}(s)+\widehat{\varphi_{{}_{\\!\mathrm{q}}}\\!}(s),$
(6) $\displaystyle\widehat{\chi_{{}_{\mathrm{v}}}\\!}(s)$ $\displaystyle=$
$\displaystyle{\color[rgb]{0,0,0}\frac{1}{s\left(s+\tau\widehat{f}(s)\right)+1},}$
(7) $\displaystyle\widehat{\chi_{{}_{\mathrm{q}}}\\!}(s)$ $\displaystyle=$
$\displaystyle{\color[rgb]{0,0,0}\left(s+\tau\widehat{f}(s)\right)\widehat{\chi_{{}_{\mathrm{v}}}\\!}(s)},$
(8) $\displaystyle\widehat{\varphi_{{}_{\\!\mathrm{q}}}\\!}(s)$
$\displaystyle=$
$\displaystyle\widehat{\chi_{{}_{\mathrm{v}}}\\!}(s)\,\widehat{\xi}(s),$ (9)

where $v_{{}_{\mathrm{0}}}$ and $q_{{}_{\mathrm{0}}}$ are the initial velocity
and position, respectively. Functions $\widehat{\chi_{{}_{\mathrm{v}}}\\!}(s)$
and $\widehat{\chi_{{}_{\mathrm{q}}}\\!}(s)$ are the so-called
susceptibilities.

Inverting Eq. (6) and averaging with the Maxwell velocity distribution gives
the solution of the GLE for any value of $\tau$,

$\displaystyle q(t)$ $\displaystyle=$
$\displaystyle\left<q[\lambda(t)]\right>+\varphi_{{}_{\\!\mathrm{q}}}\\!(t),$
(10) $\displaystyle\left<q[\lambda(t)]\right>)$ $\displaystyle=$
$\displaystyle
q_{{}_{\mathrm{0}}}\,\chi_{{}_{\mathrm{q}}}\\!(t)+{\color[rgb]{0,0,0}\int_{0}^{t}dy\,\chi_{{}_{\mathrm{v}}}\\!(t-y)\,\lambda(y)},$
(11) $\displaystyle\varphi_{{}_{\\!\mathrm{q}}}\\!(t)$ $\displaystyle=$
$\displaystyle{\color[rgb]{0,0,0}\int_{0}^{t}dy\,\chi_{{}_{\mathrm{v}}}\\!(t-y)\,\xi(y)}.$
(12)

The optimal protocol is obtained by minimizing Eq. (3) with respect to
$\lambda(t)$. This is accomplished by substituting it into the work equation,
Eq. (10), to obtain a functional both local, where its value at a given moment
depends only on that point, and non-local through its integration in the
considered time interval. Each contribution has its associated Euler-Lagrange
equation (ELE) that according to Eq. (4) gives the first variation of
$W[\lambda(t)]$ in terms of its functional derivative. Setting it to zero
allows the optimal protocol to be determined. It was already addressed in the
appendix of the previous work Colmenares and Paredes-Altuve (2021) for an
equation similar to Eq. (10) that obeys a dynamic given by the classical
Langevin equation. Thus, fitting the mentioned procedure to the GLE, the first
functional derivative associated to Eq. (4) given by,

$\displaystyle\int_{0}^{t}dy\,\lambda(y)\,A_{{}_{4}}(t,y)$ $\displaystyle=$
$\displaystyle-A_{{}_{1}}(t)+A_{{}_{2}}(t)$ (13) $\displaystyle-$
$\displaystyle\int_{0}^{t_{{}_{\mathrm{f}}}}dy\,\lambda(y)\,A_{{}_{3}}(t,y),$
$\displaystyle A_{{}_{1}}(t)$ $\displaystyle=$ $\displaystyle
q_{{}_{\mathrm{0}}}\,\dot{\chi}_{{}_{\mathrm{q}}}\\!(t),$ (14) $\displaystyle
A_{{}_{2}}(t)$ $\displaystyle=$
$\displaystyle{\color[rgb]{0,0,0}\lambda_{{}_{\mathrm{f}}}\,\chi_{{}_{\mathrm{v}}}\\!(t_{{}_{\mathrm{f}}}-t)},$
(15) $\displaystyle A_{{}_{3}}(t,y)$ $\displaystyle=$
$\displaystyle{\color[rgb]{0,0,0}\chi_{{}_{\mathrm{v}}}\\!^{\prime}(y-t)},$
(16) $\displaystyle A_{{}_{4}}(t,y)$ $\displaystyle=$
$\displaystyle{\color[rgb]{0,0,0}\left[\dot{\chi}_{{}_{\mathrm{v}}}\\!(t-y)-\chi_{{}_{\mathrm{v}}}\\!^{\prime}(y-t)\right]}\\!.$
(17)

Taking the Laplace transform of the convolution and inverting gives the
desired expression for $\lambda(t)$, that is,

$\displaystyle\lambda(t)$ $\displaystyle=$
$\displaystyle\phi_{{}_{1}}(t)-\int_{0}^{t_{{}_{\mathrm{f}}}}dy\,\phi_{{}_{2}}(t,y)\,\lambda(y),$
(18) $\displaystyle\phi_{{}_{1}}(t)$ $\displaystyle=$
$\displaystyle\mathcal{L}^{-1}\left[\frac{\widehat{A}_{2}(s)-\widehat{A}_{1}(s)}{\widehat{A}_{5}(s,t)}\right],$
(19) $\displaystyle\phi_{{}_{2}}(t,y)$ $\displaystyle=$
$\displaystyle\mathcal{L}^{-1}\left[\frac{\widehat{A}_{3}(s,y)}{\widehat{A}_{5}(s)}\right],$
(20) $\displaystyle\widehat{A}_{5}(s)$ $\displaystyle=$ $\displaystyle
s\left(\widehat{\chi_{{}_{\mathrm{v}}}\\!}(s)-\mathcal{L}\left\\{\chi_{{}_{\mathrm{v}}}\\!(-t)\right\\}\right).$
(21)

Equation 18 incorporates the whole history of the dynamics when the field is
turned on. Functions $\phi_{{}_{i}}(t)$ are well-defined but unfortunately are
long expressions to show them explicitly. They compile together the effect of
the memory kernel on the protocol. Fortunately, it appears to be defined as a
Fredholm integral equation of the second kind for which the Newman method
Arfken (1970) is suitable to solve it using $\phi_{{}_{1}}(t)$ as the seed of
the iteration process. The convergence is reached after a convenient number of
steps as it will be shown in Sec. III. This is the first mayor result.

Once Eqs. (18 is numerically solved and substituted back into Eq. (11), then
the mechanical work can be obtained from Eq. (5).

###  II.2 GOLE dynamics

The generalized overdamped equation is obtained by setting $\ddot{q}(t)=0$ in
Eq. (2). Its solution can be found in a similar way as for the GLE by
transforming it to the Laplace space. It gives
$\widehat{q}(s)={\color[rgb]{0,0,0}\widehat{\chi}_{{}_{\mathrm{ov}}}(s)(q_{{}_{\mathrm{0}}}\,\tau\,\widehat{f}(s)+\widehat{\lambda}(s)+\widehat{\varphi}_{{}_{\mathrm{ov}}}(s))}$
with
$\widehat{\chi}_{{}_{\mathrm{ov}}}(s)={\color[rgb]{0,0,0}1/(\tau\,s\,\widehat{f}(s)+1})$
being the corresponding susceptibility. From the latter
${\color[rgb]{0,0,0}\tau\,\widehat{f}(s)\widehat{\chi}_{{}_{\mathrm{ov}}}(s)}=\left(1-\widehat{\chi}_{{}_{\mathrm{ov}}}(s)\right)/s$
which after substituting it into $\widehat{q}(s)$ and inverting finally
renders

$\displaystyle q_{{}_{\mathrm{ov}}}\\!(t)$ $\displaystyle=$
$\displaystyle\left<q_{{}_{\mathrm{ov}}}\\![\lambda(t)]\right>+\varphi(t),$
(22) $\displaystyle\left<q_{{}_{\mathrm{ov}}}\\![\lambda(t)]\right>$
$\displaystyle=$ $\displaystyle q_{{}_{\mathrm{0}}}\left(1-\chi(t)\right)$
(23) $\displaystyle+$
$\displaystyle\int_{0}^{t}dy\,\chi_{{}_{\mathrm{ov}}}(t-y)\,\lambda(y),$
$\displaystyle\varphi(t)$ $\displaystyle=$
$\displaystyle\int_{0}^{t}dy\,\chi_{{}_{\mathrm{ov}}}(t-y)\,\xi(y),$ (24)
$\displaystyle\chi_{{}_{\mathrm{ov}}}(t)$ $\displaystyle=$
$\displaystyle{\color[rgb]{0,0,0}\mathcal{L}^{-1}\left[\frac{1}{\tau\,s\,\widehat{f}(s)+1}\right],}$
(25) $\displaystyle\chi(t)$ $\displaystyle=$
$\displaystyle\int_{0}^{t}dy\,\chi_{{}_{\mathrm{ov}}}(y).$ (26)

This solution, however, has an associated stationary PDF inconsistent with the
Boltzmann-Gibbs distribution as was proved by Nascimiento and Morgado
Nascimento and Morgado (2019). The reason lies on the physical consequences of
the the memory kernel on the dynamics. Physically, it gives an account of the
particle movement through the flow lines of the bath fluid. However, for
overdamped systems the inertia is irrelevant since the particle velocity
relaxation is very fast reaching the same of the reservoir. The interaction of
particle with the bath is almost instantaneous which is not considered in the
GLE by merely setting the inertial term equals to zero. In order to correct
this physical deficiency, these authors proposed to add an extra white noise
$\eta(t)$ with intensity $\gamma_{{}_{0}}$ and correlation
$\left<\eta(t)\eta(y)\right>=2\gamma_{{}_{0}}T\delta(t-y)$ to the non-inertial
GLE. In this circumstance, GOLE’s correct noise must be

$\varphi_{{}_{\\!\mathrm{ov}}}(y)=\int_{0}^{t}dy\,\chi_{{}_{\mathrm{ov}}}(t-y)\left(\xi(t)+\eta(t)\right),$
(27)

so that according to the FDT and that Eq. (22) reaches the Boltzmann-Gibbs
distribution at equilibrium, the extra term $\gamma_{{}_{0}}\,\delta(t)$ must
be added to $f(t)$ Nascimento and Morgado (2019). These arguments supported by
a detailed mathematical analysis were subscribed by Di Terlizzi et al.Di
Terlizzi  _et al._ (2020) in their GOLE assuming the optical trap is moved at
constant rate. They used $2\,\gamma_{{}_{0}}\,\delta(t)$ instead to give
account of the experiment by Berner et al. in viscoelastic fluids Berner  _et
al._ (2018).

Then, assuming the particle is rigid, replacing Eq. (23) into Eq. (3) and
recalling that $\chi_{{}_{\mathrm{ov}}}(t)$ is determined from a kernel given
by the sum of $\gamma_{{}_{0}}\,\delta(t)$ plus the exponential, the
corresponding generalized overdamped functions appearing in Eq. (4) are given
by:

$\displaystyle L_{\mathrm{loc}}(\lambda,\dot{\lambda},t)$ $\displaystyle=$
$\displaystyle\dot{\lambda}(t)\lambda(t)-q_{{}_{\mathrm{0}}}\\!\left(1-\chi(t)\right)\dot{\lambda}(t),$
(28) $\displaystyle L_{\mathrm{nl}}(\lambda,\dot{\lambda},t)$ $\displaystyle=$
$\displaystyle-\dot{\lambda}(t)\\!\int_{0}^{t}dy\,\chi_{{}_{\mathrm{ov}}}\\!(t-y)\,\lambda(y).$
(29)

Thus, the local part of $\delta W$ is simply

$\displaystyle\delta W_{\mathrm{loc}}$ $\displaystyle=$
$\displaystyle\int_{0}^{t_{{}_{\mathrm{f}}}}dt\,\dot{\lambda}(t)\,\delta\lambda(t)$
(30) $\displaystyle+$
$\displaystyle\int_{0}^{t_{{}_{\mathrm{f}}}}\,dt\left(\lambda(t)-q_{{}_{\mathrm{0}}}\left(1-\chi(t)\right)\right)\delta\dot{\lambda}(t),$

where the functional dependence was omitted for shortness.

Similarly, the non-local contribution is obtained using the general guidelines
of the procedure applied to the underdamped. It reads

$\displaystyle\delta W_{\mathrm{\\!nl}}$ $\displaystyle=$
$\displaystyle-\int_{0}^{t_{{}_{\mathrm{f}}}}dt\,\delta\lambda(t)\int_{t}^{t_{{}_{\mathrm{f}}}}dy\,\chi_{{}_{\mathrm{ov}}}(y-t)\,\lambda^{\prime}(y)$
(31) $\displaystyle-$
$\displaystyle\int_{0}^{t_{{}_{\mathrm{f}}}}dt\,\delta\dot{\lambda}(t)\int_{0}^{t}dy\,\chi_{{}_{\mathrm{ov}}}(t-y)\,\lambda(y).$

Setting to zero the sum of the associated Euler-Lagrange equations of these
two contributions and splitting the integral involving the limits
$\\{t,t_{{}_{\mathrm{f}}}\\}$ gives the final result

$\displaystyle\int_{0}^{t}dy\,A_{{}_{1}}(t,y)\\!$ $\displaystyle-$
$\displaystyle\\!\int_{0}^{t_{{}_{\mathrm{f}}}}\\!\\!dy\,A_{{}_{2}}(t,y)\lambda^{\prime}(y)\\!+\\!A_{{}_{3}}(t)=0,$
(32) $\displaystyle A_{{}_{1}}(t,y)$ $\displaystyle=$
$\displaystyle\lambda^{\prime}(y)\,\chi_{{}_{\mathrm{ov}}}(y\\!-\\!t)\\!+\\!\lambda(y)\,\dot{\chi}_{{}_{\mathrm{ov}}}(t\\!-\\!y),$
(33) $\displaystyle A_{{}_{2}}(t,y)$ $\displaystyle=$
$\displaystyle\chi_{{}_{\mathrm{ov}}}(y-t),$ (34) $\displaystyle
A_{{}_{3}}(t)$ $\displaystyle=$ $\displaystyle
q_{{}_{\mathrm{0}}}\,\dot{\chi}(t)+\chi_{{}_{\mathrm{ov}}}(0)\,\lambda(t).$
(35)

To proceed further, the Laplace transform of Eq. (32) renders

$\displaystyle\widehat{\lambda}(s)$ $\displaystyle=$
$\displaystyle\widehat{\Psi}_{{}_{1}}(s)+\int_{0}^{t_{{}_{\mathrm{f}}}}dy\,\widehat{\Psi}_{{}_{2}}(s,y)\,\lambda^{\prime}(y),$
(36) $\displaystyle\widehat{\Psi}_{{}_{1}}(s)$ $\displaystyle=$
$\displaystyle\frac{\lambda_{{}_{\mathrm{0}}}\,\widehat{\chi}_{{}_{\mathrm{ov}}}^{\,\star}(s)-q_{{}_{\mathrm{0}}}\left(s\,\chi(s)-\chi(0)\right)}{\widehat{A}_{{}_{4}}(s)},$
(37) $\displaystyle\widehat{\Psi}_{{}_{2}}(s,y)$ $\displaystyle=$
$\displaystyle\frac{\widehat{A}_{2}(s,y)}{\widehat{A}_{4}(s)},$ (38)
$\displaystyle\widehat{A}_{{}_{4}}(s)$ $\displaystyle=$ $\displaystyle
s\left(\widehat{\chi}_{{}_{\mathrm{ov}}}^{\,\star}(s)+\widehat{\chi}_{{}_{\mathrm{ov}}}(s)\right),$
(39) $\displaystyle\widehat{\chi}_{{}_{\mathrm{ov}}}^{\,\star}(s)$
$\displaystyle=$ $\displaystyle\mathcal{L}\\{\chi_{{}_{\mathrm{ov}}}(-t)\\}.$
(40)

Inverting Eq (36) finally gives

$\lambda(t)=\Psi_{{}_{\\!\\!1}}(t)+\int_{0}^{t_{{}_{\mathrm{f}}}}dy\,\Psi_{{}_{\\!\\!2}}(t,y)\,\lambda^{\prime}(y),$
(41)

with the two functions $\Psi_{{}_{\\!\\!1}}(t)$ and $\Psi_{{}_{\\!\\!2}}(t,y)$
given by Mat :

$\displaystyle\Psi_{{}_{\\!\\!1}}(t)$ $\displaystyle=$
$\displaystyle\frac{1}{2}\left[\lambda_{{}_{\mathrm{0}}}-q_{{}_{\mathrm{0}}}+\left(q_{{}_{\mathrm{0}}}+\lambda_{{}_{\mathrm{0}}}\right)\left(\frac{\alpha}{\omega_{{}_{1}}}t+\frac{\left(1+\omega_{{}_{1}}^{2}\right)}{\omega_{{}_{1}}^{2}\,\sqrt{\alpha\,\gamma_{{}_{0}}}}\,\sinh\left(\sqrt{\frac{\alpha}{\gamma_{{}_{0}}}}\,\omega_{{}_{1}}\,t\right)\right)\right],$
(42) $\displaystyle\Psi_{{}_{\\!\\!2}}(t,y)$ $\displaystyle=$
$\displaystyle\frac{\exp\left[-\frac{(1+\omega_{{}_{1}}^{2})}{2\,\gamma_{{}_{0}}}y\right]}{2\,\sqrt{\alpha\,\gamma_{{}_{0}}}\,\omega_{{}_{1}}^{2}\,\omega_{{}_{2}}}\Bigg{\\{}\cosh\left(\frac{\omega_{{}_{2}}}{2\,\gamma_{{}_{0}}}y\right)\left[\omega_{{}_{1}}\,\omega_{{}_{2}}\,\sqrt{\alpha\,\gamma_{{}_{0}}}\left(\omega_{{}_{1}}^{2}+\alpha\,t\right)+\left(1+\omega_{{}_{1}}^{2}\right)^{2}\sinh\left(\sqrt{\frac{\alpha}{\gamma_{{}_{0}}}}\,\omega_{{}_{1}}\,t\right)\right]$
(43) $\displaystyle+$
$\displaystyle\left.\sinh\left(\frac{\omega_{{}_{2}}}{2\,\gamma_{{}_{0}}}y\right)\left[\sqrt{\alpha\,\gamma_{{}_{0}}}\,\,\omega_{{}_{1}}\,\left(1+\omega_{{}_{1}}^{2}+\alpha^{2}\,\gamma_{{}_{0}}\left(\gamma_{{}_{0}}+t\right)\right)-2\,\sqrt{\alpha\,\gamma_{{}_{0}}}\,\,\omega_{{}_{1}}\left(1+\omega_{{}_{1}}^{2}\right)\cosh\left(\sqrt{\frac{\alpha}{\gamma_{{}_{0}}}}\,\omega_{{}_{1}}\,t\right)\right.\right.$
$\displaystyle-$
$\displaystyle\left(1+\omega_{{}_{1}}^{2}\right)^{2}\sinh\left(\sqrt{\frac{\alpha}{\gamma_{{}_{0}}}}\,\omega_{{}_{1}}\,t^{2}\right)\bigg{]}\Bigg{\\}},$

where $\omega_{{}_{1}}=\sqrt{1+\alpha\,\gamma_{{}_{0}}}$ and
$\omega_{{}_{2}}=\sqrt{4+\alpha^{2}\,\gamma_{{}_{0}}^{2}}$. The integral
equation (36) converges after a moderate number of iterations if
$\Psi_{{}_{\\!\\!1}}^{\prime}(y)$ is used as the initial guess for
$\lambda^{\prime}(y)$. Functions
$\left<q_{{}_{\mathrm{ov}}}\\![\lambda(t)]\right>$,
$\left<\dot{q}_{{}_{\mathrm{ov}}}\\![\lambda(t)]\right>$ and $W$ are long
expressions to display so they are numerically determined.

In the next section, it will be presented an auto-consistent procedure to
calculate the optimal protocol and minimum work for GLE and GOLE. They will be
compared with those obtained for the same system using the MLE Colmenares and
Paredes-Altuve (2021) as a standard or point of reference. Likewise, it will
discuss the results in a thermodynamic basis.

##  III Discussion of results

The chosen initial conditions are the same as in Ref. Colmenares and Paredes-
Altuve (2021) where the Markovian equations lead to a negative mechanical
work. Namely, $\\{q_{{}_{\mathrm{0}}},\lambda_{{}_{\mathrm{f}}}\\}=\\{1,0\\}$.
The final time of the protocol application $t_{{}_{\mathrm{f}}}$ is 12 and the
friction coefficient for the GLE ranges from the very underdamped $\gamma=1$
up to the overdamped described by GOLE. The decay constant $\alpha$ in the GLE
is assumed to decrease with the friction to incorporate intuitively the
retardation in the particle translation as $\gamma$ gets higher. It was chosen
to get $W$ of the same order of magnitude for the considered initial set and
all $\gamma$ values. It should be mentioned that the GLE algorithm is unstable
for high $\tau$ due to high oscillations in the functions $\phi_{1}(t)$ and
$\phi_{2}(t)$ and because of it, it is replaced for that of the GOLE.
Additionally, the selected $\alpha$ for a given $\gamma$ is critical otherwise
the oscillations mentioned before show unbound values.

The plots were calculated for increasing values of $\tau$ and segmented in
terms of $\gamma$ and $\alpha$. The $\alpha$ for a given $\gamma$ was
arbitrarily chosen in such a way that it generate mechanical work done by the
particle in the first three cases, $\gamma=\\{1,3,6\\}$ and done by the field
in GOLE. Besides the increasing $\tau$ condition, an extra restriction was
used to tune up the selected $\alpha$. It is related to the total work
involved over the entire time range such that for the assigned decay constant,
this quantity stays around a fixed meaningful value. It was fixed to be around
$-0.94$ for $\gamma=\\{1,3,6\\}$ while for GOLE was around 15. In this way,
the color tags used in the plots are for $\gamma$ of 1 (black), 3 (red), 6
(blue) and GOLE (green) with their respective associated $\alpha$ values
$\\{0.16,0.13,0.09\\}$ for the first three cases. For GOLE a high $\gamma$ and
a low $\alpha$ is required, thus we set $\alpha=0.001$, $\gamma=1/\alpha$ and
$\gamma_{{}_{0}}=5$. Using the above procedure to select $\alpha$, the
numerical determination of $\lambda(t)$ and the desired $W$ is very efficient
in the GLE.

It is important to mention that GOLE’s calculations are very sensitive to  the
value of $\gamma_{{}_{0}}$. Convergence is reached for a given value in such a
way that the memory effect avoids taking over the non-inertial dynamics. For
other values besides this critical value, the solution is out of the
prescribed bounds with no physical sense at all.

Since no physical arguments have been used to assign the correct value of the
decay constant, computer simulation seems to be the only way to get the proper
information about its magnitude. Various methods could be used to accomplish
this purpose. The first could be the technique developed by Ceriotti et al.
Ceriotti  _et al._ (2009) to use the GLE to perform an isothermal molecular
dynamics simulation and the second, that from Stella et al. Stella  _et al._
(2014) where the colored exponential noise is consistently generated.

Based on this reasoning, the optimal driving are depicted in Fig. 1.

Figure 1: Time- dependent optimal protocol in terms of their defining
parameters $\\{\gamma,\alpha\\}$ of $\\{1,0.16\\}$ (black), $\\{3,0.13\\}$
(red) and $\\{6,0.09\\}$ (blue). The green curve is for GOLE with
$\alpha=0.001$ and $\gamma=1/\alpha$ to fulfill the non-inertial character of
the description. The set
$\\{q_{{}_{\mathrm{0}}},\lambda_{{}_{\mathrm{f}}}\\}=\\{1,0\\}$. See text for
details.

For the underdamped, they are non-linear functions of time unlike the straight
lines predicted in the Markovian case Colmenares and Paredes-Altuve (2021).
GOLE shows an opposite behavior in comparison with the overdamped MLE where
$\lambda(t)$ decreases with time. Together, they have a non-vanishing value at
$t=0$ not shown because of the scale.

The minimum work is shown in Fig. 2 for the optimal protocols of Fig. 1.

Figure 2: Minimum work for the protocols of Fig. 1.

The differences between the results shown in this graph compared to those
obtained with the MLE counterpart become now more than evident. Particularly,
there is no observation of discontinuities in the very high underdamped as
seen in Fig. 1 of Ref. Colmenares and Paredes-Altuve (2021) with the
$\gamma=1$ curve (black). The similarity in behavior of the GLE underdamped
curves is due to the requirement that $W$ has a fixed value around $-0.94$
throughout the time interval. Furthermore, excluding the overdamped, GLE does
not predict a smooth and negative $W$ as in MLE, but has an almost periodic-
like behavior with positive values in some ranges of the protocol application
time. Extractable work of these characteristics were also found by Bylicka et
al. Bylicka  _et al._ (2016) in open quantum systems with memory effects.
Furthermore, GOLE predicts a behavior contrary to the Markovian.

One way to analyze it is through the Fokker-Planck formalism. It requires
knowing the reduced conditional probability distribution function (PDF),
$p(q,t|q_{{}_{\mathrm{0}}})$. The method to obtain it is fully described in
Secs. 2 y 3.1 of Ref. Olivares and Colmenares (2016). Making it suitable to
the problems at hand, it is a Gaussian given by

$\displaystyle p(q,t|q_{{}_{\mathrm{0}}})$ $\displaystyle=$
$\displaystyle\frac{1}{\sqrt{2\,\pi\,\sigma^{2}(t)}}\exp\left[-\frac{\left(q-H(t)\right)^{2}}{2\,\sigma^{2}(t)}\right],$
(44)

where the standard deviation $\sigma_{\mathrm{q}}^{2}(t)$ and function $H(t)$
are defined respectively as:

$\displaystyle\hskip 28.45274pt\mathrm{GLE}$ $\displaystyle\sigma^{2}(t)$
$\displaystyle=$ $\displaystyle
2\int_{0}^{t}\\!\\!\\!\\!dy\int_{0}^{y}\\!\\!\\!\\!dz\,\left<\varphi_{{}_{\\!\mathrm{v}}}\\!(y)\varphi_{{}_{\\!\mathrm{v}}}\\!(z)\right>+T\chi_{{}_{\mathrm{v}}}^{2}(t),$
(45) $\displaystyle H(t)$ $\displaystyle=$
$\displaystyle\left<q[\lambda(t)]\right>),$ (46)
$\displaystyle\left<\varphi_{{}_{\\!\mathrm{v}}}\\!(t)\varphi_{{}_{\\!\mathrm{v}}}\\!(s)\right>$
$\displaystyle=$ $\displaystyle
2\,\gamma\,\alpha\,T\,\int_{0}^{t}dy\,\dot{\chi}_{{}_{\mathrm{v}}}\\!(t-y)$
(47) $\displaystyle\times$
$\displaystyle\int_{0}^{s}dz\,\dot{\chi}_{{}_{\mathrm{v}}}\\!(s-z)\,f(|y-z|),$
$\displaystyle\hskip 28.45274pt\mathrm{GOLE}$ $\displaystyle\sigma^{2}(t)$
$\displaystyle=$ $\displaystyle
2\int_{0}^{t}dy\int_{0}^{y}dz\,\left<\varphi_{{}_{\\!\mathrm{ov}}}(y)\varphi_{{}_{\\!\mathrm{ov}}}(z)\right>,$
(48) $\displaystyle H(t)$ $\displaystyle=$
$\displaystyle\left<q_{{}_{\mathrm{ov}}}\\![\lambda(t)]\right>$ (49)

The correlation
$\left<\varphi_{{}_{\\!\mathrm{ov}}}(t)\varphi_{{}_{\\!\mathrm{ov}}}(s)\right>$
can be determined assuming the noises are uncorrelated with
$\left<\eta(t)\eta(s)\right>=\gamma_{{}_{0}}\,\delta(t-s)$ and
$\left<\xi(t)\xi(s)\right>=f(|t-s))$.

The Fokker-Planck equation whose solution is the above Gaussian is determined
by applying the procedure originally designed by Adelman and Garrison Adelman
(1976) and detailed shown in Ref. Colmenares (2018). It renders,

$\displaystyle\frac{\partial p(q,t)}{\partial t}$ $\displaystyle=$
$\displaystyle-\frac{\partial J(t,q)}{\partial q},$ (50) $\displaystyle
J(q,t)$ $\displaystyle=$
$\displaystyle\Omega(t)\,q\,p(q,t)-\frac{1}{2}D(t)\frac{\partial
p(q,t)}{\partial q},$ (51) $\displaystyle\Omega(t)$ $\displaystyle=$
$\displaystyle\frac{\dot{H}(t)}{H(t)},$ (52) $\displaystyle D(t)$
$\displaystyle=$
$\displaystyle\dot{\sigma}^{2}(t)-2\,\sigma^{2}(t)\,\Omega(t),$ (53)

where $D(t)$ is the effective diffusion constant. It is shown in Fig. 3 for
$T=10$ along with the average position and standard deviation.

Figure 3: Function $H(t)=\left<q(t)\right>$, effective diffusion constant
$D(t)$ and standard deviation $\sigma^{2}(t)$ for GLE and GOLE. Color scheme
as of Fig. 1.

Regarding the average position $H(t)$, we see that  all curves tend to diverge
from the potential center being more evident in GOLE due to the instantaneous
hits received from the particles of the thermal bath. This is also manifested
in the protocols shown in Fig. 1 through the final jumps.

The diffusion constant appears to be a smooth function for all the considered
friction constants unlike the strong Markovian underdamped which shows
discontinuities in $D(t)$ because of the ratio $\Omega(t)$. This was
extensively discussed in a previous paper Colmenares (2019) for the Markovian
harmonic oscillator. The effect of average velocity and position on the above
relationship makes the GLE’s $D(t)$ a continuous function of time. Then we can
conclude that the MLE applied to the highly underdamped case cannot cope with
the cusps in $D(t)$ because the averages describe a physical situation without
the translational memory. This means that the mechanical work has no clear
physical meaning because it also has unexpected cusps Colmenares and Paredes-
Altuve (2021).

The standard deviation $\sigma_{\mathrm{q}}^{2}(t)$ results of the combined
effects of the noise correlation function, the effective diffusion constant
and the ratio $\Omega(t)$. Its non-monotonicity is primarily due to the
effects of the oscillations shown by the susceptibilities and the functions
$\varphi_{1}(t)$ and $\varphi_{2}(t,y)$. They are canceled out in the
overdamped giving a linear function because of the preeminence of the
intensity of white noise over the memory kernel $f(t)$. What is important to
notice is that GOLE being a more accurate theory, it gives a totally different
$\sigma^{2}(t)$ compared with the saturated curve of the MLE description
Colmenares (2019).

###  III.1 Entropy production rate

To give an account of the thermodynamic cost invested in doing work by or
against the field, entropy $S=-\int
dq\,p(q,t|q_{{}_{\mathrm{0}}})\,\ln[p(q,t|q_{{}_{\mathrm{0}}})]$ is the proper
property to be called upon. This property was detailed analyzed by Seifert
along the trajectory of the Brownian particle Seifert (2005) and subsequently
by van den Broeck and Esposito in terms of the processes responsible for the
breaking of the detailed balance once the protocol is initiated Van den Broeck
and Esposito (2010).

The standard deviation is crucial on the behavior of the total entropy $S(t)$
and its derivative, the entropy rate $\Sigma(t)$. They can be easily
determined from Eq. (44) and plotted in Fig. 4.

Figure 4: Total Entropy $S(t)$ and entropy production rate $\Sigma(t)$ for GLE
and GOLE probability distributions. Color scheme as of Fig. 1.

It is observed a negativity temporarily character of $\Sigma(t)$. We must
emphasize, however, they are legitimate outcomes own of the non-Markovian
dynamics despite of the provisions of the Second Law. In fact, because the
bath degrees of freedom (BDF) in the theory are not explicitly considered but
have been described by the static friction coefficient, their real
contribution to the entropy balance are excluded so their explicit
incorporation to the total will provide the required positivity. This kind of
“thermodynamical inconsistent” behavior is inherent to non-Markovian systems.
It was proved by Marcantoni et al. Marcantoni  _et al._ (2017) for open
quantum systems described by non-Markovian maps and by Strasberg and Esposito
Strasberg and Esposito (2019) for a dynamics based on an effective master
equation where the non-Markovian character is incorporated through an
transition matrix depending on the initial state having negative entries and
fulfilling the Chapmann-Kolmogorov equation. Moreover, Bonança et al. Bonança
_et al._ (2021) have shown for the electrical conduction in metals that
temporal negative entropy production also spring out as a direct consequence
of Ohm’s law. In the language of the GLE, the Ohmic driving of the dynamics is
the non-Markovian dissipation term Schramm  _et al._ (1985).

Therefore, it is plausible to assume that the control of the thermodynamic
outcome by the standard deviation shown in Fig. 3 along with the time-
dependent diffusion term, are responsible for the system showing a negative
rate of entropy production. This property will be further analyzed next.

In a system in equilibrium, the correlation of the fluctuations for small
external disturbances is given in terms of the response function. This is not
the case in stationary systems out of equilibrium because the detailed balance
relation is broken and therefore there is a continuous degradation of energy
to the thermal reservoir. This behavior occurs along the trajectory which in
turn reduces to the conventional FDT if the velocity fluctuations are given in
terms of the local mean velocity Speck and Seifert (2006).

In the detailed balance restoration, van den Broeck and Sposito Van den Broeck
and Esposito (2010) claim that entropy production involves the interplay of
two different mechanisms. One is the adiabatic (ad) recovering of the temporal
symmetry due to any given constraint and the second one, the persistent non-
adiabatic (na) influence of the field thought the driving.

The two mentioned contributions to the entropy rate are given by Eqs. (27) and
(28) of Ref. Van den Broeck and Esposito (2010), which applied to Eq. (51)
gives,

$\displaystyle\Sigma_{\mathrm{na}}(t)$ $\displaystyle=$
$\displaystyle\frac{2}{D(t)}\int_{-\infty}^{\infty}dq\,p(q,t|q_{{}_{\mathrm{0}}})\left(\frac{J(q,t)}{p(q,t|q_{{}_{\mathrm{0}}})}\right)^{\\!2},$
(54) $\displaystyle=$
$\displaystyle\frac{D(t)}{2\,\sigma^{2}(t)}\\!+\\!2\,\Omega(t)\\!\left(\\!1\\!+\\!\Omega(t)\frac{H^{2}(t)\\!+\\!\sigma^{2}(t)}{D(t)}\right)\\!\\!,$
(55) $\displaystyle\Sigma_{\mathrm{ad}}(t)$ $\displaystyle=$ $\displaystyle
0.$ (56)

There is not an adiabatic contribution since the term
$J^{\mathrm{st}}(q,t)/p^{\mathrm{st}}(q,t)$ vanishes because of the stationary
(st) PDF

$p^{\mathrm{st}}(q,t)=\sqrt{-\frac{\Omega(t)}{\pi\,D(t)}}\exp\left[\frac{\Omega(t)}{D(t)}q^{2}\right].$
(57)

The negativity of $\Omega(t)/D(t)$ is guaranteed so the system reaches the
steady state. It is important to recall that $\Sigma_{\mathrm{na}}(t)$ agrees
with that by Seifert Seifert (2005) in which the rate is calculated along the
trajectory.

Hence, the protocol is the only source of entropy production to restore the
symmetry breaking since there are no constrictions of any kind on the
dynamics.

The entropy rate $\Sigma_{\mathrm{na}}(t)$ is a complicated function of time
showing sudden changes at given times for all values of the friction constant.
Furthermore, $\int_{0}^{t_{{}_{\mathrm{f}}}}dt\,\Sigma_{\mathrm{na}}(t)$ is
positive. These two facts are due to the non-rigorous entropy balance from the
effective management of the BDFs. For the matter of presentation, it is only
shown in Fig. 5 the results for $\gamma=3$ and GOLE.

Figure 5: Non-adiabatic entropy production rate $\Sigma_{\mathrm{na}}(t)$ for
$\gamma=3$ and GOLE.

It can be seen, that indeed, the temporally negative entropy production rates
in discrete non-Markovian systems also appear in the non-Markovian Langevin
equation in both the inertial and underdamped descriptions.

In order to test this result then it is important to determine the entropy
production, $\Sigma_{\mathrm{ex}}(t)$, due to the heat dissipated to the bath.
Using any of the approaches Seifert (2005); Van den Broeck and Esposito (2010)
it is found that

$\displaystyle\Sigma_{\mathrm{ex}}(t)$ $\displaystyle=$
$\displaystyle-2\frac{\Omega(t)}{D(t)}\int_{-\infty}^{\infty}dq\,q\,J(q,t),$
(58) $\displaystyle=$
$\displaystyle-\frac{\Omega(t)}{D(t)}\left(D(t)\\!+\\!2\,\Omega(t)\left(H^{2}(t)\\!+\\!\sigma(t)\right)\right)\\!,$
(59)

which added to $\Sigma_{\mathrm{na}}(t)$ gives $\Sigma(t)$ as it should be.

These results do not give any clues about the effects that the actual
production of entropy has on the ability of the system to do work against the
external field. The negative entropy rate is a consequence of the ineffective
management of the bath degrees of freedom requiring a reformulation of the FDT
as was mentioned in the introduction. In particular, the observation by
Kutvonen et al. Kutvonen  _et al._ (2015) that Markov limit of non-Markovian
systems leads to additional entropy components. The absence of it in the MLE
can provide a secondary route to thermodynamically understand the cusp
appearing in the mechanical work mentioned above. Recently, Hsiang and Hu
Hsiang and Hu (2021) derived the FDT for systems subjected to non-equilibrium
dynamics driven by an external protocol.

Despite these arguments, the results provided by GLE and GOLE are
mathematically correct. It corroborates the work of Van den Broeck and
Esposito, where a dissipative-driven dynamics of coarse-grained non-Markov
systems, the rate of entropy production is naturally negative Strasberg and
Esposito (2019).

However, since the classical GLE is only an approximation of a real system
based on an FDT framed in linear response theory, then it is plausible to
assume that a ”more consistent” formulation of it will provide results to be
in accordance with the prescriptions of the second law.

##  IV Concluding remarks

Summarizing, the auto-consistent inclusion of memory effects definitely
changes the expected value of the work done by the particle in comparison with
the Markovian counterpart Colmenares and Paredes-Altuve (2021). The initial
position other than the equilibrium values is responsible for the particle to
do work against the external field except in the generalized overdamped. The
calculations for other initial conditions will show as in Colmenares and
Paredes-Altuve (2021) that work is done by the field. The generalized
overdamped protocol is, as in the Markovian, a general result no matter the
friction coefficient values.

Unlike the results for the Markovian where the optimal protocol is given by a
differential equation, the inclusion of memory leads to an integral equation
for the driving in either description.

It has been provided explicit calculations for the entropy rate of non-
Markovian dynamics extending those of Ref. Van den Broeck and Esposito (2010)
and reinforcing the findings of Ref. Strasberg and Esposito (2019).

It was found the protocol is the only source the system to show a negative
entropy rate which clearly violates the second law. It is due to consider in
the formulation of the theory that the effects of the degrees of freedom of
the bath are described as a unique function of the temperature and as a
consequence putting aside important contributions to the total entropy
balance.

A reformulation of the FDT is called upon to make the GLE approach
thermodynamically consistent.

## References

  * Imparato  _et al._ (2007) A. Imparato, L. Peliti, G. Pesce, G. Rusciano, and A. Sasso, Phys. Rev. E 76, 050101(R) (2007).
  * Ciliberto (2017) S. Ciliberto, Phys. Rev. X 7, 021051 (2017).
  * Schmiedl and Seifert (2007) T. Schmiedl and U. Seifert, Phys. Rev. Letts. 98, 108301 (2007).
  * Abreu and Seifert (2011) D. Abreu and U. Seifert, Europhys. Lett. 94, 10001 (2011).
  * Colmenares and Paredes-Altuve (2021) P. J. Colmenares and O. Paredes-Altuve, Phys. Rev. E 104, 034115 (2021).
  * Daldrop  _et al._ (2017) J. O. Daldrop, B. G. Kowalik, and R. R. Netz, Phys. Rev. X 7, 041065 (2017).
  * Vroylandt and Monmarché (2022) H. Vroylandt and P. Monmarché, J. Chem. Phys. 156, 244105 (2022).
  * Lisý and Tóthová (2019) V. Lisý and J. Tóthová, Results Phys. 12, 1212 (2019).
  * Di Cairano (2022) L. Di Cairano, J. Phys. Commun. 6, 015002 (2022).
  * Balakrishnan (1979) B. Balakrishnan, Pramāna 12, 301 (1979).
  * Cicotti and Ryckaert (1981) G. Cicotti and J. P. Ryckaert, J. Stat. Phys. 26, 73 (1981).
  * Olivares and Colmenares (2016) W. Olivares and P. J. Colmenares, Physica A 458, 76 (2016).
  * Costa  _et al._ (2003) I. V. L. Costa, R. Morgado, M. V. B. T. Lima, and F. A. Oliveira, Europhys. Lett. , 173 (2003).
  * Seifert and Speck (2010) U. Seifert and T. Speck, Europhys. Lett. 89, 10007 (2010).
  * Kutvonen  _et al._ (2015) A. Kutvonen, T. Ala-Nissila, and J. Pekola, Phys. Rev. E 92, 012107 (2015).
  * Paredes-Altuve  _et al._ (2016) O. Paredes-Altuve, E. Medina, and P. J. Colmenares, Phys. Rev. E 94, 062111 (2016).
  * Di Terlizzi  _et al._ (2020) I. Di Terlizzi, F. Ritort, and M. Baiesi, J. Stat. Phys. 181, 1609 (2020).
  * Ferialdi and Bassi (2012) L. Ferialdi and A. Bassi, Europhys. Lett. 98, 30009 (2012).
  * Zwanzig (2001) R. Zwanzig,  _Nonequilibrium Statistical Mechanics_ (Oxford University Press, 2001).
  * Banik  _et al._ (2002) S. K. Banik, B. C. Bag, and D. S. Ray, Phys. Rev. E 65, 051106 (2002).
  * Reichl (1998) L. E. Reichl,  _A Modern Course in Statistical Physics_ , 2nd ed. (John Wiley & Sons, New York, 1998).
  * Sekimoto (2010) K. Sekimoto,  _Stochastic Energetics._, Lecture Notes in Physics 799. (Springer, Heidelberg, 2010).
  * Seifert (2012) U. Seifert, Rep. Prog. Phys. 75, 126001 (2012).
  * Arfken (1970) G. Arfken,  _Mathematical methods for physicists_ , 2nd ed. (Academic Press, New York-London, 1970).
  * Nascimento and Morgado (2019) E. S. Nascimento and W. A. M. Morgado, Europhys. Lett. 126, 10002 (2019).
  * Berner  _et al._ (2018) J. Berner, B. Müller, J. R. Gomez-Solano, M. Krüger, and C. Bechinger, Nat. Commun. 9, 1 (2018).
  * (27) “Mathematica™ 11.2.0.0 computer package was used for the algebraic, numerical and graphics manipulations,” .
  * Ceriotti  _et al._ (2009) M. Ceriotti, G. Bussi, and M. Parrinello, Phys. Rev. Letts. 102, 020601 (2009).
  * Stella  _et al._ (2014) L. Stella, C. D. Lorenz, and L. Kantorovich, Phys. Rev. B 89, 134303 (2014).
  * Bylicka  _et al._ (2016) B. Bylicka, M. Tukiainen, D. Chruściński, J. Piilo, and S. Maniscalco, Scientific Reports 6, 1 (2016).
  * Adelman (1976) S. A. Adelman, J. Chem. Phys. 64, 124 (1976).
  * Colmenares (2018) P. J. Colmenares, Phys. Rev. E 97, 052126 (2018), There is a misprint in Eq. (12), $\xi(t)$ must be $\varphi_{\mathrm{v}}(t)$.
  * Colmenares (2019) P. J. Colmenares, Phys. Rev. E 100, 052117 (2019), The term inside the summation in Eq. (10d) must be $(\nu_{n}\exp[-\nu_{n}t])[(\nu_{n}^{2}+\lambda_{1}^{2}](\nu_{n}^{2}+\lambda_{2}^{2}))]^{-1}$.
  * Seifert (2005) U. Seifert, Phys. Rev. Letts. 95, 040602 (2005).
  * Van den Broeck and Esposito (2010) C. Van den Broeck and M. Esposito, Phys. Rev. E 82, 011144 (2010).
  * Marcantoni  _et al._ (2017) S. Marcantoni, S. Alipour, F. Benatti, R. F. Floreanini, and A. T. Rezakhani, Scientific Reports 7, 1 (2017).
  * Strasberg and Esposito (2019) P. Strasberg and M. Esposito, Phys. Rev. E 99, 012120 (2019).
  * Bonança  _et al._ (2021) M. V. S. Bonança, P. Nazé, and S. Deffner, Phys. Rev. E 103, 012109 (2021).
  * Schramm  _et al._ (1985) P. Schramm, R. Jung, and H. Grabert, Phys. Lett. A 107, 385 (1985).
  * Speck and Seifert (2006) T. Speck and U. Seifert, Europhys. Lett. 74, 391 (2006).
  * Hsiang and Hu (2021) J.-T. Hsiang and B.-L. Hu, Annals of Physics 433, 168594 (2021).