# Analytical Insights and Universal Behavior in Fast Thermal Equilibration
Protocols

Diego Rengifo<EMAIL_ADDRESS>Gabriel Téllez
<EMAIL_ADDRESS>Physics Department, Universidad de los Andes, Bogotá,
Colombia.

###### Abstract

When a system deviates from equilibrium, it is possible to manipulate and
control it to drive it towards equilibrium within a finite time $t_{f}$, even
reducing its natural relaxation time scale ${\tau_{\text{relax}}}$. Although
numerous theoretical and experimental studies have explored these shortcut
protocols, few have yielded analytical results for the probability
distribution of work, heat and entropy production. In this paper, we propose a
two-step protocol that captures the essential characteristics of more general
protocols and has analytical solution for the relevant thermodynamic
probability distributions. Additionally, we demonstrate that for very short
protocol duration $t_{f}\ll{\tau_{\text{relax}}}$, all protocols exhibit
universal behavior in their shape and for the ratio of probability
distribution functions of positive and negative work, heat and the entropy
production.

Stochastic thermodynamics, control, fast thermalization, brownian motion

##  I Introduction

Almost every thermodynamic system in nature is out of equilibrium. Equilibrium
states are therefore not common but desirable. When a system is out of
equilibrium and left without any external intervention, it takes a finite time
to reach an equilibrium state. Let us call ${\tau_{\text{relax}}}$ the time
scale of relaxation. This time is an intrinsic characteristic of any physical
system and depends on various factors, such as the underlying interactions,
which are encoded in the friction, transport coefficients, the external
parameters (if any), and the temperature [1].

Technology is advancing very rapidly, and one of its trends is to create and
control smaller devices. This miniaturization of devices has led to an
increasing interest in the development of engineered techniques that can
shorten the natural time scale for relaxation between equilibrium states.
These procedures are designed to connect equilibrium states through a
particular protocol that is considerably shorter than the natural
equilibration time. Such techniques have been inspired by the so-called
shortcut to adiabaticity [2, 3, 4]. Since then, the term Engineered Swift
Equilibration (ESE) has been coined to describe such protocols [5]. These
protocols are also known as “shortcut to isothermality” [6] or “swift state-
to-state transformations” [7].

Several protocols have been established for the ESE process, including those
for frictionless atom cooling in harmonic traps [8, 9] and nanosystems as
micromechanical oscillators in contact with a thermostat, both in the
overdamped or underdamped regime [5, 10]. These procedures enable the creation
of a non-equilibrium state that can be controlled and manipulated, thereby
allowing the exploration of novel physical phenomena [7].

Overall, the ESE process has shown great promise in the field of non-
equilibrium statistical mechanics and has opened up new avenues for research
in the development of novel techniques for controlling and manipulating the
dynamics of physical systems, especially for nanodevices. Due to the
miniaturization of chips, robots, and devices, understanding controlled
dynamics has become mandatory in modern technology. The suitable framework to
tackle such small systems is stochastic thermodynamics.

To advance our understanding of how protocols accelerate system equilibration,
it is crucial to comprehend the stochastic thermodynamics of small systems.
Therefore, it is vital to calculate the probability distribution functions
(PDF) of work, heat, and entropy production for these distinct processes.

In this paper, we consider a Brownian particle, the system, immerse in a
viscose medium, the environment at temperature $T$, whose particles have a
smaller size than the Brownian particle. This justifies the validity of the
commonly called overdamped regime, which is particularly applicable to
colloidal particles. The particle is subjected to a time-dependent and
externally controlled harmonic-type potential, which allows for precise
manipulation and control. An experimental realization of this system has been
made with colloidal particles trapped with laser beams [11, 12, 5] and has
been used to build microscopic heat engines [13, 14].

When the stiffness of the potential is controlled through ESE protocols, the
probability distributions of the relevant thermodynamic quantities like work
and heat cannot generally be calculated analytically. Therefore, their
theoretical study is limited to numerical simulations. In this paper, we
introduce the two-step protocol as a novel approach. Our protocol enables us
to derive exact analytical results. Remarkably, the analytical solution of the
two-step protocol captures the essential characteristics of broader protocols,
exhibiting a certain level of universality.

The structure of this article is as follows. In the first part, inspired by
the ideas presented in [5], we have developed a protocol that establishes a
connection between two equilibrium states, theoretically reducing the
equilibration time. The work from Ref. [5] also proposed a similar protocol
and carried out experimental realizations, achieving a two-order-of-magnitude
reduction in equilibration time. However, due to the complex energetics
involved, a detailed study of these protocols has proven to be extremely
challenging [15].

In the second part, we propose a simplified toy model called the Two-Step
Protocol (TSP), which allows for an analytical solution of the probability
distribution functions of relevant thermodynamic quantities. We calculate the
work, heat, and entropy production PDF for the TSP and we check the validity
of the Jarzynski equality [16, 17], Crooks relation [18, 19], and the entropy
production fluctuation theorem [20] for this specific protocol. A
comprehensive comparison with more general protocols reveals some universal
features of the shape of the protocols as well as in the ratio of PDF
$P_{A}(A)/P_{A}(-A)$ (with $A$ being the work, heat and entropy production).
These relations are different from the usual fluctuations theorems [20, 21,
22, 16, 17, 18, 19], as both $P_{A}(A)$ and $P_{A}(-A)$ refer to the PDF of a
forward process (we do not consider a backward process with the direction of
time reversed).

##  II ESE Protocol

In this section, we present a review of the concepts underlying the ESE
protocols [5, 7]. The central idea behind an ESE process is to construct a
customized time-dependent protocol, denoted as $\lambda(t)$, for the
externally controlled parameter. This protocol is specifically designed to
guide the system from the initial equilibrium state characterized by
$\lambda_{i}$ to the desired final equilibrium state characterized by
$\lambda_{f}$ within a finite time interval $t_{f}$ that is shorter than the
system’s equilibration time ${\tau_{\text{relax}}}$.

Let us consider a Brownian particle trapped in a harmonic potential with time-
dependent stiffness, given by

$U(x,t)=\frac{1}{2}k(t)x^{2}.$ (1)

Here $x$ is the position of the particle and the protocol $\lambda(t)$ is
essentially characterized by the stiffness $k(t)$. Initially, the particle is
in equilibrium with stiffness $k_{i}$. By controlling the stiffness of the
potential $k(t)$, the system reaches a final equilibrium state with stiffness
$k_{f}$. However, there are infinitely many possible protocols that can
achieve this transformation, so we can impose constraints to obtain the
desired solution. Since the goal is to reduce the equilibration time through
external control, the problem falls within the realm of optimization theory
[23, 24].

A solution to this problem has been presented in [5], where both experimental
and theoretical results have been obtained in the overdamped limit, i.e., the
acceleration term in the Langevin equation has been neglected. In this limit,
the corresponding Langevin equation is given by

$\dot{x}=-\frac{k(t)}{\gamma}x+\sqrt{2\mathcal{D}}\xi(t),$ (2)

where $\xi(t)$ is a Gaussian white noise with zero average and auto-
correlation function
$\expectationvalue{\xi(t)\xi(t^{\prime})}=\delta(t-t^{\prime})$. Here,
$\beta=1/(k_{B}T)$ and $\gamma$ are the inverse temperature and friction
coefficient, respectively, and $\mathcal{D}$ is the diffusion constant
$\mathcal{D}=k_{B}T/\gamma$. From now on, we will assume that the overdamped
limit is valid.

From the Langevin equation, we can infer that the natural time scale is given
by ${\tau_{\text{relax}}}=\gamma/k_{f}$. It is useful to use a set of
adimensional variables associated with this time scale and the final stiffness
$k_{f}$ of the protocol: $\tilde{t}=t/{\tau_{\text{relax}}}$,
$\tilde{x}(\tilde{t})=x(t)/\sqrt{\mathcal{D}{\tau_{\text{relax}}}}$, and
$\tilde{k}(\tilde{t})=k(t)/k_{f}$. The natural energy scale is $k_{B}T$, given
by $\tilde{U}=U/(k_{B}T)=(1/2)\tilde{k}(\tilde{t})(\tilde{x}(\tilde{t}))^{2}$.
With this set of units, the reduced Langevin equation reads

$\frac{d\tilde{x}}{d\tilde{t}}=-\tilde{k}(\tilde{t})\tilde{x}(\tilde{t})+\sqrt{2}\tilde{\xi}(\tilde{t}),$
(3)

where $\tilde{\xi}(\tilde{t})=\sqrt{{\tau_{\text{relax}}}}\xi(t)$, which
satisfies $\langle\tilde{\xi}(\tilde{t})\rangle=0$ and
$\langle\tilde{\xi}(\tilde{t})\tilde{\xi}(\tilde{t}^{\prime})\rangle=\delta(\tilde{t}-\tilde{t}^{\prime})$.
Using these units, it is clear that there are only two parameters for our
problem: the initial stiffness with respect to the final one,
$\tilde{k}_{i}=k_{i}/k_{f}$, and the target duration of the protocol compared
to the relaxation time scale, $\tilde{t}_{f}=t_{f}/{\tau_{\text{relax}}}$.
From now on, we will use these non-dimensional units and remove the tilde to
lighten the notation.

The system is in an equilibrium state at the beginning, $t_{i}=0$, and the
goal of the ESE protocol is that it reaches a new equilibrium state at time
$t_{f}$. As a result, the position probability distributions can be described
by Gaussian distributions, which are characterized by

$\displaystyle\expectationvalue{x(0)}$ $\displaystyle=0$ (4)
$\displaystyle\sigma_{i}^{2}=\expectationvalue{x^{2}(0)}$
$\displaystyle=\frac{1}{k_{i}}$ (5)
$\displaystyle\sigma_{f}^{2}=\expectationvalue{x^{2}(t_{f})}$
$\displaystyle=\frac{1}{k_{f}}=1.$ (6)

As is well-known [25], the Langevin equation can be mapped to a Fokker-Planck
equation

$\frac{\partial P}{\partial t}=\frac{\partial}{\partial
x}\left[k(t)xP\right]+\frac{\partial^{2}P}{\partial x^{2}}$ (7)

for the probability density function $P(x,t)$ of position. To solve this
equation, we performed a Fourier transform

$\hat{G}(p,t)=\int_{-\infty}^{\infty}P(x,t)e^{-ipx}\,dx$ (8)

leading to an equation of the form

$\frac{\partial}{\partial t}(\ln\hat{G})=k(t)p\frac{\partial}{\partial
p}(\ln\hat{G})-p^{2}.$ (9)

The combination $\ln(\hat{G})$ is the cumulant generating function whose
Taylor series is

$\ln{\hat{G}(p,t)}=\sum_{n=0}^{\infty}\chi_{n}(t)\frac{(-ip)^{n}}{n!},$ (10)

where $\chi_{n}(t)$ are the cumulants of $P(x,t)$. The average value and
variance correspond to $n=1$ and $n=2$, respectively. The insertion of this
expansion leads to the derivation of a set of ordinary differential equations
governing the time evolution of all cumulants for the position. Notably, the
only non-trivial equation corresponds to $n=2$ ($\chi_{2}(t)=\langle
x(t)^{2}\rangle$)

$\dot{\chi}_{2}(t)+2k(t)\chi_{2}(t)=2.$ (11)

The essence of the ESE protocols is to propose a particular functional form
for the variance $\chi_{2}$ such that it satisfies the initial (5) and final
conditions (6): $\chi_{2}(0)=1/k_{i}$ and $\chi_{2}(t_{f})=1/k_{f}=1$. After
that, the appropriate stiffness $k(t)$ can be extracted from equation (11).
Thus, the process of finding the appropriate stiffness is done in a reverse-
engineered way.

To ensure a smooth transition to equilibrium, two additional conditions can be
optionally imposed on $\chi_{2}(t)$, given by

$\dot{\chi}_{2}(0)=\dot{\chi}_{2}(t_{f})=0.$ (12)

We have four conditions that must be satisfied, and hence, we need four
parameters to adjust. In addition, we aim to optimize the work done on the
system during the process by adding an extra fifth parameter that allows us to
tune it in such a way that the work is minimized. Based on this consideration,
we propose a solution in the form of a fourth-degree polynomial for
$\chi_{2}$:

$\chi_{2}(t)=A_{0}+A_{1}t+A_{2}t^{2}+A_{3}t^{3}+A_{4}t^{4}.$ (13)

Finding the parameters $A_{i}$ using the boundary conditions (5)-(6), the
variance is given by

$\chi_{2}(t)=\frac{1}{k_{i}}-\frac{\Delta
k}{k_{i}k_{f}}\left(3s^{2}-2s^{3}\right)+A_{4}t_{f}^{4}\left(s^{2}-2s^{3}+s^{4}\right),$
(14)

where $s=t/t_{f}$ is the reduced time and $\Delta k=k_{f}-k_{i}$. Replacing
this in (11), the function $k(t)$ must be

$k(t)=k_{i}\frac{1+\frac{3}{t_{f}}\frac{\Delta
k}{k_{i}}(s-s^{2})-\epsilon\frac{k_{f}}{k_{i}}\frac{1}{t_{f}}(s-3s^{2}+2s^{3})}{1-\frac{\Delta
k}{k_{f}}(3s^{2}-2s^{3})+\epsilon(s^{2}-2s^{3}+s^{4})},$ (15)

where $\epsilon=A_{4}k_{i}t_{f}^{4}$. The stochastic work [26, 27] done on the
system is given by

$W=\frac{1}{2}\int_{0}^{t_{f}}x^{2}\frac{dk}{dt}dt.$ (16)

Computing the average and replacing (15), we obtain

$\displaystyle\expectationvalue{W}$
$\displaystyle=\frac{1}{2}\ln\left(\frac{k_{f}}{k_{i}}\right)+\frac{1}{4t_{f}}\frac{k_{f}}{k_{i}}\eta,$
$\displaystyle=\Delta F+\expectationvalue{W_{\text{irr}}}$ (17)

where $\eta$ has the expression

$\eta=\int_{0}^{1}\frac{\left[\frac{-6\Delta
k}{k_{f}}(s-s^{2})+\epsilon(2s-6s^{2}+4s^{3})\right]^{2}}{1-\frac{\Delta
k}{k_{f}}(3s^{2}-2s^{3})+\epsilon(s^{2}-2s^{3}+s^{4})}ds,$ (18)

and $\Delta F=\frac{1}{2}\ln(k_{f}/k_{i})$ is the free energy difference
between the final and initial state and $\expectationvalue{W_{\text{irr}}}$ is
the irreversible work. The average work can be interpreted as a function of
$\epsilon$, and numerical methods can be used to find the value of $\epsilon$
that minimizes the average work. It is worth noting that the irreversible work
$\expectationvalue{W_{\text{irr}}}$ exhibits an inverse relationship with the
duration of the protocol, $t_{f}$. This inverse correlation is consistent with
our expectation for a process aimed at accelerating equilibration. The
external control exerted on the system necessitates additional work to achieve
equilibrium within a shorter timeframe.

Since the choice for the functional form of $\chi_{2}$ is arbitrary, there
exist an infinite number of solutions to the fast equilibration problem. For
example, in [5], a polynomial of degree 3 form was proposed for the inverse of
$\chi_{2}(t)$, and the corresponding experimental realization using optical
tweezers has been done. The experimental realization has proved that their
protocol has shortened the relaxation time by two orders of magnitude. In
order to make a comparison, we use the same values reported in [5], which in
adimensional units are $k_{i}=0.5$ and $t_{f}=1/30$, as shown in Figure 1. The
average work for the protocol $k(t)$ defined in Eq. (15), optimizing
$\epsilon$ to minimize the average work, is $\expectationvalue{W}=6.52$, while
the result published in [5] is $\expectationvalue{W}=6.71$. It is worth noting
that the optimal work protocol can be obtained using Pontryagin’s principle
(see Eq. (20)), with a result of $\expectationvalue{W}_{\text{opt}}=5.49$ [5,
23, 28, 29, 24].

Figure 1: Plots of protocols (15), Martinez et al. ([5]) and TSP (23). Note
that the maximum value of $k(t)$ occurs around $k_{\text{max}}=16k_{f}$. This
high value is necessary to reduce the equilibration time. The central step in
the protocol is $k_{m}$ that is solution of eq (26). The difference between
the stiffness average of the Martinez et al. [5], the stiffness average of the
protocol $k(t)$ Eq. (15) and $k_{m}$ are small, and it is not visible in the
plot. The parameters used are $k_{i}=1/2$ and $t_{f}=1/30$.

Once we have proposed a form for $\chi_{2}(t)$, it is possible to find the
corresponding protocol $k(t)$ using Eq. (11). It is clear that there are
infinitely many protocols with different characteristics. For example, for a
variance of the form

$\chi_{2}(t)=\frac{1}{k_{i}}(1+ct)^{2},$ (19)

the protocol is

$k_{\text{op}}(t)=\frac{k_{i}}{(1+ct)^{2}}-\frac{c}{1+ct},$ (20)

where $c=\frac{1}{t_{f}}(\sqrt{k_{i}/k_{f}}-1)$ and time interval
$t\in(0,t_{f})$. This protocol has the peculiarity that it minimizes the work,
as it has been shown in Refs. [28, 23]. A second example is a linear relation
for the variance $\chi_{2}(t)$ that satisfy the boundary conditions
Eqs.(5)-(6):

$\chi_{2}(t)=\frac{1}{k_{i}}+(\frac{1}{k_{f}}-\frac{1}{k_{i}})\frac{t}{t_{f}},$
(21)

whose protocol $k_{L}(t)$ is of the form for $t\in(0,t_{f})$

$k_{L}(t)=\frac{2-(\frac{1}{k_{f}}-\frac{1}{k_{i}})\frac{1}{t_{f}}}{2\left(\frac{1}{k_{i}}+(\frac{1}{k_{f}}-\frac{1}{k_{i}})\frac{t}{t_{f}}\right)}.$
(22)

Both protocols are discontinuous at $t=0$ and $t=t_{f}$ and they will be
called optimal and linear, respectively. Studying the energetics [26] of these
protocols, as well as almost all other protocols, can only be achieved through
numerical simulations [15]. This raises the question of whether it is possible
to design a protocol for which analytical solutions can be obtained for the
PDF of work, heat, and entropy production. The next section addresses this
question.

##  III  Two-step protocol

In this section, we propose a new fast thermalization protocol that provides
analytical solutions for the PDF of relevant thermodynamic quantities. Upon
careful examination of Figure 1, it is evident that the stiffness must be
significantly greater than both the initial and final values of $k$. This is
because, in order to accelerate the equilibration process, the stiffness needs
to be increased to reduce the characteristic time scale evolution of the
system. Therefore, it is reasonable to conclude that the large values of
stiffness play a critical role in achieving a shortcut to adiabaticity.

Keeping this in mind, we propose a two-step protocol (TSP) defined as follows:

$\displaystyle k(t)=\begin{cases}k_{i}&\text{if }\,t\leq 0,\\\ k_{m}&\text{if
}\,0<t<t_{f},\\\ k_{f}&\text{if }\,t\geq t_{f}.\end{cases}$ (23)

As we will see, this toy model captures the essential features of general
protocols. Using this stiffness, the variance in the time interval $0\leq
t\leq t_{f}$ is obtained by solving Eq. (11) with the initial condition (5):

$\chi_{2}(t)=\sigma_{X}(t)^{2}=\frac{1}{k_{m}}+\left(\frac{1}{k_{i}}-\frac{1}{k_{m}}\right)e^{-2k_{m}t}.$
(24)

Using this variance, the expression for the probability distribution of the
position at any time ($0\leq t\leq t_{f}$) can be written as:

$P_{X}(x_{t},t)=\frac{1}{\sqrt{2\pi}\sigma_{X}(t)}\exp\left(-\frac{x_{t}^{2}}{2\sigma_{X}^{2}(t)}\right).$
(25)

Now, by forcing the system to arrive at a final equilibrium state at
$t=t_{f}$, the variance $\chi_{2}(t_{f})$ has to be $1/k_{f}$. This leads us
to a self-consistency equation for $k_{m}$ given by:

$\exp(-2k_{m}t_{f})=\frac{\displaystyle\frac{1}{k_{f}}-\frac{1}{k_{m}}}{\displaystyle\frac{1}{k_{i}}-\frac{1}{k_{m}}}.$
(26)

This equation is transcendental and can be solved numerically to determine the
value of $k_{m}$. Utilizing this value, the TSP is a control protocol that
joins the equilibrium state at $k_{i}$ with the equilibrium state at $k_{f}$
in a finite time $t_{f}$. This equation is not symmetric under time inversion
of the protocol, i.e., if the system evolves from $k_{f}$ to $k_{i}$, and the
equilibrium distributions are desired at initial and final times, the value of
$k_{m}$ changes because the roles of $k_{i}$ and $k_{f}$ need to be
interchanged in Eq. (26), and therefore the system does not follow the same
evolution under the reverse direction of time. This fact is related to the
entropy production, as discussed below.

This approach differs from that of the previous section and previous works [5,
30], as here we propose a specific form of the protocol and determine the
corresponding variance. Thus, this process can be referred to as direct
engineering.

Although Eq. (26) cannot be solved analytically, we can derive the asymptotic
behavior of $k_{m}$ for large and small values of the target final time
$t_{f}$. If $t_{f}\gg 1$, $\exp(-2k_{m}t_{f})\to 0$. From the right-hand side
of (26), we find $k_{m}=k_{f}=1$. This is to be expected as if we have an
infinite amount of time to let the system equilibrate ($t_{f}\to\infty$), we
can fix $k_{m}=k_{f}$ and wait for the system to reach equilibrium at this new
value of the stiffness.

More interesting and pertinent is the behavior when $t_{f}\to 0$. In this
limit, $k_{m}\to\infty$, which means that the shorter the target time is, the
larger the stiffness needs to be. In this limit, the right-hand side of (26)
converges to a finite value $k_{i}/k_{f}$. Therefore, we deduce that
$k_{m}t_{f}$ must remain finite, and we have:

$k_{m}\sim\frac{1}{2t_{f}}\ln\frac{k_{f}}{k_{i}}.$ (27)

We observe that $k_{m}$ is inversely proportional to $t_{f}$. This quantifies
the previous observation that to achieve fast relaxation to equilibrium, it is
necessary to significantly increase the stiffness. Since the typical
relaxation time at a fixed stiffness $k_{m}$ is of the order $1/k_{m}$, the
TSP matches a transient relaxation time of the order $1/k_{m}$ to the target
duration $t_{f}$ of the accelerated protocol. The specific relation is given
by (27).

By expanding the right-hand side of (26) in powers of $1/k_{m}$, one can
systematically compute the corrections to (27). For example, the next-to-
leading order term is given by:

$k_{m}t_{f}=\frac{1}{2}\ln\frac{k_{f}}{k_{i}}+t_{f}\frac{k_{f}-k_{i}}{\ln(k_{f}/k_{i})}+O(t_{f}^{2}).$
(28)

It turns out that the leading order (27) is universal for more general
protocols, if we replace $k_{m}$ by the average value of $k(t)$ over the
protocol duration

$k_{m}=\frac{1}{t_{f}}\int_{0}^{t_{f}}k(t)\,dt.$ (29)

Indeed, using Eq. (11) for a general protocol, we find that the average of
$k(t)$ is

$k_{m}=\frac{1}{2t_{f}}\ln\frac{\chi_{2}(0)}{\chi_{2}(t_{f})}+\int_{0}^{t_{f}}\frac{1}{\chi_{2}(t)}\frac{dt}{t_{f}}.$
(30)

Since $\chi_{2}$ is a continuous function of $t$, applying the mean value
theorem and using the initial and final conditions (5)-(6) leads to

$k_{m}=\frac{1}{2t_{f}}\ln\frac{k_{f}}{k_{i}}+\frac{1}{\chi_{2}(t^{*})}$ (31)

for some $t^{*}\in[0,t_{f}]$. Figure 2 illustrates this by showing a plot of
$k_{m}t_{f}$ as a function of $t_{f}$ for different values of $k_{i}$ and the
ESE protocols discussed in the previous section and Ref. [5]. As $t_{f}\to 0$,
all the protocols converge to the same value given by (27). However, the next-
to-leading corrections differ from one protocol to another.

Figure 2:  The average value $k_{m}$ of the stiffness multiplied by $t_{f}$ is
plotted for several protocols: the Two-Step Protocol (TSP) given by Eq.(23),
the protocol proposed by Martinez et al.[5], and the protocol $k(t)$ given by
Eq.(15). It can be observed that at short target times $t_{f}$, all protocols
exhibit the same asymptotic behavior, as predicted by Eq.(27).

###  III.1 Work distribution for two-step protocol

Since the TSP seems to be analytically tractable and exhibits some general
features of more complex ESE protocols, we now proceed to use it as a basis
for understanding the stochastic energetics of the fast relaxation protocols.
More specifically, in this section, we want to calculate the work probability
distribution. Using (16), we obtain the work done on the system $W_{D}$ during
the process

$\displaystyle
W_{D}=\left\\{\begin{array}[]{lcc}\frac{1}{2}(k_{m}-k_{i})x_{i}^{2}&\quad
0<t<t_{f},\\\ \\\
\frac{1}{2}(k_{m}-k_{i})x_{i}^{2}+\frac{1}{2}(k_{f}-k_{m})x_{f}^{2}&\quad
t>t_{f}.\end{array}\right.$ (35)

Since the initial and final positions are stochastic variables with normal
distributions, the work done is also a stochastic variable. To compute the
work distribution function, we consider the following cases:

####  III.1.1 For $0<t<t_{f}$

To study this time interval, we need to determine the position probabilities
at two points described by $(x_{i},t_{i}=0)$ and $(x_{t},t)$. Therefore, the
joint probability distribution (JPD) can be expressed as follows:

$P_{2}(x_{i},t_{i};x_{t},t)=P(x_{t},t|x_{i},t_{i})P(x_{i},t_{i}).$ (36)

Computing $P_{2}$ is a straightforward task; $P(x_{t},t|x_{i},t_{i})$ follows
an Ornstein-Uhlenbeck process distribution [31]:

$P(x_{t},t_{t}|x_{i},t_{i})=\frac{1}{\sqrt{2\pi\sigma_{m}(t)^{2}}}\exp{-\frac{(x_{t}-x_{i}e^{-k_{m}t})^{2}}{2\sigma_{m}(t)^{2}}}$
(37)

where

$\sigma_{m}(t)^{2}=\frac{1-e^{-2k_{m}t}}{k_{m}}.$ (38)

The initial Gaussian probability distribution is

$P(x_{i},t_{i})=\frac{1}{\sqrt{2\pi}\sigma_{i}}\exp{-\frac{x_{i}^{2}}{2\sigma_{i}^{2}}}$
(39)

with $\sigma_{i}^{2}=1/k_{i}$ (eq. (25) at $t=0$).

To compute the probability density function of the work, we adopt the
following approach: while there exist an infinite number of paths connecting
the points $(x_{i},t_{i})$ and $(x_{t},t)$, we focus solely on the paths that
yield a specific amount of work, denoted as $W$, regardless of the initial and
final points. Hence, the PDF of the work can be expressed as follows:

$P_{W}(W,t)=\int_{-\infty}^{\infty}dx_{i}dx_{t}\delta(W-W_{D})P_{2}(x_{i},t_{i};x_{t},t).$
(40)

Since the work $W_{D}$ is independent of $x_{t}$ for $t<t_{f}$, the work
probability is stationary, and the previous equation simplifies to:

$P_{W}(W,t)=\int_{-\infty}^{\infty}dx_{i}\delta\left(W-\frac{(k_{m}-k_{i})x_{i}^{2}}{2}\right)P(x_{i},t_{i}).$
(41)

If $W/(k_{m}-k_{i})<0$, the argument of the Dirac delta distribution never
vanishes, resulting in $P(W,t)=0$. However, when $W/(k_{m}-k_{i})\geq 0$,
there are two roots to the equation $W-(k_{m}-k_{i})x_{i}^{2}/2=0$, yielding:

$P_{W}(W,t)=\sqrt{\frac{k_{i}}{\pi(k_{m}-k_{i})W}}\exp\left(-\frac{k_{i}}{k_{m}-k_{i}}W\right).$
(42)

This result can also be obtained using the methods described in [15, 27].
Returning to (41), one can use the Fourier representation of the Dirac
distribution

$\delta(W-(k_{m}-k_{i})x_{i}^{2}/2)=\int_{-\infty}^{+\infty}e^{iz(W-(k_{m}-k_{i})x_{i}^{2}/2)}\,\frac{dz}{2\pi}.$
(43)

to compute the characteristic function of $P_{W}(W,t)$. Using (43) into (41),
the resulting gaussian integral over $x_{i}$ can be performed to obtain

$P_{W}(W,t)=\frac{1}{2\pi}\int_{-\infty}^{\infty}dz\hat{P}_{W}(z)e^{iWz}$ (44)

where the characteristic function is

$\hat{P}_{W}(z)=\frac{1}{\sqrt{1+ia_{1}z}}$ (45)

with $a_{1}=(k_{m}-k_{i})/k_{i}$. From this we can obtain all the cumulants
$c_{n}$ of $W$,

$\displaystyle\ln\hat{P}_{W}(z)$
$\displaystyle=\sum_{n=1}^{\infty}\frac{(-iz)^{n}}{n!}c_{n}$
$\displaystyle=\sum_{n=1}^{\infty}\frac{(-ia_{1}z)^{n}}{n}.$ (46)

Thus

$c_{n}=\frac{1}{2}(n-1)!a_{1}^{n}.$ (47)

In particular, the average work and variance for this time interval are

$\displaystyle\expectationvalue{W}$ $\displaystyle=$
$\displaystyle\frac{k_{m}-k_{i}}{2k_{i}},$ (48) $\displaystyle\sigma_{W}^{2}$
$\displaystyle=$ $\displaystyle\frac{(k_{m}-k_{i})^{2}}{2k_{i}^{2}}.$ (49)

Note that $\sigma_{W}^{2}=2\expectationvalue{W}^{2}$. From the characteristic
function one can check that the work distribution satisfies Jarzynski equality
[16] as it should:

$\left<e^{-W}\right>=\hat{P}_{W}(i)=\sqrt{\frac{k_{i}}{k_{m}}}=e^{-\Delta F},$
(50)

where $\Delta F=\frac{1}{2}\ln(k_{m}/k_{i})$ is the free energy difference for
the oscillator with final stiffness $k_{m}$ (in the time interval $0<t<t_{f}$)
and initial stiffness $k_{i}$.

####  III.1.2 For $t\geq t_{f}$

In the temporal range $t\geq t_{f}$, it is necessary consider the three-point
joint probability distribution

$P_{3}(x_{i},t_{i};x_{f},t_{f};x,t)=P(x,t|x_{f},t_{f})P_{2}(x_{f},t_{f};x_{i},t_{i}).$
(51)

Again, the work distribution is stationary, hence we only need the two-points
JPD

$\displaystyle P_{2}(x_{i},t_{i};x_{f},t_{f})$
$\displaystyle=\frac{1}{2\pi\sigma_{i}\sigma_{m}(t)}\exp{-\frac{(x_{f}-x_{i}e^{-k_{m}t_{f}})^{2}}{2\sigma_{m}(t_{f})^{2}}}$
$\displaystyle\times\exp{-\frac{x_{i}^{2}}{2\sigma_{i}^{2}}}$ (52)

here $\sigma_{m}(t)$ is given by (38). Inserting the two-point JPD in an
expression similar to (40) and using the Fourier representation of the Dirac
distribution, we obtain a Gaussian integral in two variables that can easily
be computed leading us to a characteristic function of the form

$\hat{P}_{W}(z)=\frac{1}{\sqrt{1+a_{1}z+a_{2}z^{2}}}$ (53)

where the coefficients are given by

$\displaystyle a_{1}$ $\displaystyle=$
$\displaystyle\frac{i(k_{f}-k_{i})}{k_{i}}+\frac{i(k_{f}-k_{m})(k_{i}-k_{m})\sigma_{m}^{2}(t_{f})}{k_{i}}$
(54) $\displaystyle a_{2}$ $\displaystyle=$
$\displaystyle\frac{i(k_{f}-k_{m})(k_{i}-k_{m})\sigma_{m}^{2}(t_{f})}{k_{i}}$
(55)

or utilizing the consistency equation (26), these expressions can be
simplified to

$\displaystyle a_{1}$ $\displaystyle=$
$\displaystyle\frac{i(k_{f}-k_{i})k_{m}}{k_{i}k_{f}}$ (56) $\displaystyle
a_{2}$ $\displaystyle=$
$\displaystyle\frac{(k_{f}-k_{i})(k_{m}-k_{f})}{k_{f}k_{i}}$ (57)

Clearly, $P_{W}(W)$ is normalized ($\hat{P}_{W}(0)=1$) and satisfies the
Jarzynski equality [16]

$\expectationvalue{e^{-W}}=\hat{P}_{W}(i)=e^{-\Delta F},$ (58)

where now $\Delta F=\frac{1}{2}\ln(k_{f}/k_{i})$ is the Helmholtz free energy
difference between the two equilibrium states with stiffness $k_{f}$ and
$k_{i}$. From the characteristic function, we can determine any moments of the
work distribution, including the average work and variance:

$\displaystyle\expectationvalue{W}$ $\displaystyle=$
$\displaystyle\frac{k_{m}}{2}\left(\frac{1}{k_{i}}-\frac{1}{k_{f}}\right),$
(59) $\displaystyle\sigma_{W}^{2}$ $\displaystyle=$ $\displaystyle
2\expectationvalue{W}^{2}+\frac{2(k_{m}-k_{f})\expectationvalue{W}}{k_{m}}.$
(60)

It is worth noting that $\sigma_{W}^{2}>2\expectationvalue{W}^{2}$ holds true
when $(k_{m}-k_{f})\expectationvalue{W}>0$. Fig. 6 shows the evolution of the
average work and its standard deviation.

The exact computation of the inverse Fourier transform in Eq. (44) can be
achieved for the characteristic function given in Eq. (53). To accomplish
this, it is imperative to use the vertex form of a quadratic polynomial.
Subsequently, an elementary change of variable results in an integral
representation of the modified Bessel function of the second kind of order
zero $K_{0}$. The final expression is as follows:

$\displaystyle P_{W}(W)$
$\displaystyle=\frac{1}{\pi\sqrt{\sigma_{W}^{2}-2\expectationvalue{W}^{2}}}\exp{\frac{\expectationvalue{W}}{\sigma_{W}^{2}-2\expectationvalue{W}^{2}}W}K_{0}\left(\frac{\sqrt{\sigma_{W}^{2}-\expectationvalue{W}^{2}}}{\sigma_{W}^{2}-2\expectationvalue{W}^{2}}\absolutevalue{W}\right).$
(61)

The tails of this work distribution can be obtained from the asymptotic
behavior of the Bessel function $K_{0}$, they are

$\displaystyle
P_{W}(W)\sim\left\\{\begin{array}[]{lcc}\frac{1}{\sqrt{2\pi\sigma_{r}W}}\exp{\frac{-1}{\sigma_{r}+\expectationvalue{W}}W}&\quad
W\rightarrow\infty,\\\
\frac{1}{\sqrt{-2\pi\sigma_{r}W}}\exp{\frac{1}{\sigma_{r}-\expectationvalue{W}}W}&\quad
W\rightarrow-\infty.\end{array}\right.$ (64)

where $\sigma_{r}^{2}=\sigma_{W}^{2}-\expectationvalue{W}^{2}$.

To test this analytical result and compare it with the PDF for more
complicated protocols, we performed simulations using a code developed by one
of us [32]. The comparison between this exact solution and the simulation is
depicted in figure 3.

Figure 3:  The curve is the theoretical prediction for $t>t_{f}$ given by (61)
and the histogram is built simulating the Langevin dynamic of the system under
study. This histogram is for $10^{6}$ simulations with parameters
$k_{i}=k_{f}/2=1/2$ and $t_{f}=1/(30k_{f})$.

The probability density function of the work exhibits an intriguing structure
characterized by the product of two components. The first component is a
symmetric function of $W$, specifically the Bessel $K_{0}$ function with an
argument proportional to the absolute value of $W$. The second component is an
exponential function of $W$, introducing an asymmetry in the probabilities of
obtaining positive or negative work values. Since the overall process is a
compression ($k_{f}>k_{i}$), the average work on the particle is positive.
Nevertheless there are rare events in which $W<0$ but those are less frequent
than the ones where $W>0$: the tail of $P_{W}(W)$ for $W>0$ is larger than the
one for $W<0$. Mathematically, this asymmetry factor can be expressed as
follows:

$\frac{P_{W}(W)}{P_{W}(-W)}=\exp{\frac{2\expectationvalue{W}}{\sigma_{W}^{2}-2\expectationvalue{W}^{2}}W}=\exp{\frac{k_{m}}{k_{m}-k_{f}}W}\,.$
(65)

Equivalently, this implies the symmetry relation for the characteristic
function

$\hat{P}_{W}(z)=\hat{P}_{W}\left(-z-i\frac{k_{m}}{k_{m}-k_{f}}\right).$ (66)

The other protocols follow a similar relation, at least for $t_{f}\ll 1$, as
observed in the plot shown in Figure 4. This suggests that the relation given
by Eq. (65) is universal for protocols that achieve fast thermal equilibration
in a short time, where $t_{f}\ll{\tau_{\text{relax}}}=1$. However, this
universality breaks down for longer times: when $t_{f}\sim 1$, the protocols
exhibit slight deviations from the predicted behavior of the TSP protocol, as
shown in Figure 5. Nonetheless, the universality of Eq. (65) remains a
prominent feature, valid within the range of interest for fast thermalization
protocols ($t_{f}\ll 1$).

Figure 4:  The protocol given by Martinez et al. [5], the protocol described
in Eq. (15), and the linear (22) and optimal protocol (20) exhibit a relation
predicted by the TSP in Eq. (65). We have performed $10^{6}$ simulations using
$k_{i}=1/2$ and $t_{f}=1/30$. Figure 5:  For $t_{f}\sim 1$, the protocols
shown in the figure, present a deviation from the behavior described by (65).
However, we can appreciate a linear tendency. This plot has been made for
$10^{6}$ simulations with parameters $t_{f}=0.7$ and $k_{i}=1/2$. Figure 6:
The standard deviation and work average in all time interval for the TSP
protocol whose parameters are $k_{i}=1/2$ and $t_{f}=1/30$.

####  III.1.3 Reverse protocol

In the previous section, we computed the probability distribution connecting
two equilibrium states, $k_{i}$ and $k_{f}$, through a suitable choice of
$k_{m}$ as given by Eq. (26). However, our results are more general than this
specific case. The probability density function of work retains the same
expressions as before for arbitrary $k_{m}$, except that the average work and
its variance will change. In this general setup, the expressions for $t\geq
t_{f}$ are as follows:

$\displaystyle\expectationvalue{W}$ $\displaystyle=$
$\displaystyle\frac{k_{f}-k_{i}}{2k_{i}}+\frac{(k_{m}-k_{f})(k_{m}-k_{i})}{2k_{i}}\sigma_{m}^{2}(t_{f}),$
(67) $\displaystyle\sigma_{W}^{2}$ $\displaystyle=$ $\displaystyle
2\expectationvalue{W}^{2}+\frac{(k_{m}-k_{f})(k_{m}-k_{i})}{k_{i}}\sigma_{m}^{2}(t_{f}).$
(68)

For $t<t_{f}$, the expressions remain unchanged, as is shown in Eq. (48).

Let us consider a time-reversed protocol defined as

$\displaystyle k^{R}(t)=\begin{cases}k_{f}&\text{if }\,t\leq 0,\\\
k_{m}&\text{if }\,0<t<t_{f},\\\ k_{i}&\text{if }\,t\geq t_{f}.\end{cases}$
(69)

In this time-reversed protocol, we do not change the value of $k_{m}$, it is
the same as in the forward protocol. As a result, the system will not be at
thermal equilibrium at $t_{f}$, since the value of $k_{m}$ has not been
adjusted properly (to obtain thermal equilibrium at $t_{f}$ in a reversed
protocol, the roles of $k_{i}$ and $k_{f}$ had to be interchanged in Eq.
(26)).

The Crooks relation can be verified by utilizing the general expression for
$\expectationvalue{W}$ and $\sigma_{W}^{2}$ [19]:

$\frac{P_{W}(W)}{P^{R}_{W}(-W)}=\exp{W-\Delta F},\quad\text{for any
{\hbox{t}}}.$ (70)

Here, $P^{R}_{W}(-W)$ represents the work probability distribution of the
reverse protocol. The cancellation of Bessel functions occurs due to the
symmetry of the argument under the interchange of $k_{i}$ and $k_{f}$. The
factor $\Delta F$ arises from the ratio of normalization factors, while the
appearance of $W$ stems from the exponential function (refer to Appendix A for
more details).

###  III.2 Heat distribution function for TSP

Now we move to the study of another relevant thermodynamics quantity: heat. In
the realm of stochastic thermodynamics, the concept of heat, as defined in
[26], is expressed as:

$Q=\int_{0}^{t}k(t)x\dot{x}dt{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}.}$
(71)

Using the first law of thermodynamics, the heat can also be computed from
$Q_{D}(t)=U(t)-U(0)-W_{D}(t)$. For the two-step protocol (TSP), the stochastic
heat is:

$\displaystyle
Q_{D}(t)=\left\\{\begin{array}[]{lcc}\frac{1}{2}k_{m}(x_{t}^{2}-x_{i}^{2})&\quad
t\leq t_{f},\\\ \\\
\frac{1}{2}k_{m}(x_{f}^{2}-x_{i}^{2})+\frac{1}{2}k_{f}(x_{t}^{2}-x_{f}^{2})&\quad
t>t_{f}.\end{array}\right.$ (75)

It is worth noting that the heat distribution is time-dependent because it
involves $x_{t}$. This complication leads to a non-stationary heat probability
density function (PDF). Similar to the work calculation, the analysis of heat
needs to be performed in two distinct cases:

####  III.2.1 For $0\leq t\leq t_{f}$

By following a similar methodology as with the work distribution, we can
obtain the characteristic function for the heat
$\hat{P}_{Q}(z)=(1+b_{1}z+b_{2}z^{2})^{-1/2}$ with the coefficients

$\displaystyle b_{1}$ $\displaystyle=$
$\displaystyle\frac{-i(k_{m}-k_{i})k_{m}}{k_{i}}\sigma_{m}^{2}(t),$
$\displaystyle b_{2}$ $\displaystyle=$
$\displaystyle\frac{k_{m}^{2}}{k_{i}}\sigma_{m}^{2}(t),$

where $\sigma_{m}(t)$ is given at 38. From this characteristic function, we
can derive the average $\expectationvalue{Q}_{t}$:

$\expectationvalue{Q}_{t}=-\frac{k_{m}-k_{i}}{2k_{i}}\left(1-e^{-2k_{m}t}\right),$
(76)

and the variance of the heat

$\sigma_{Q}(t)^{2}=2\expectationvalue{Q}_{t}^{2}-\frac{2k_{m}\expectationvalue{Q}_{t}}{k_{m}-k_{i}}.$
(77)

It is important to observe that the average heat and variance exhibit a
characteristic time $\tau_{m}=1/k_{m}$. The characteristic function has a
similar structure to that of the work probability density function for
$t>t_{f}$. Therefore, the inverse Fourier transform can be computed
explicitly, leading to the heat probability distribution for $0\leq t\leq
t_{f}$,

$\displaystyle P_{Q}(Q,t)$
$\displaystyle=\frac{1}{\pi\sqrt{\sigma_{Q}(t)^{2}-2\expectationvalue{Q}_{t}^{2}}}\exp{\frac{\expectationvalue{Q}_{t}}{\sigma_{Q}(t)^{2}-2\expectationvalue{Q}_{t}^{2}}Q}K_{0}\left(\frac{\sqrt{\sigma_{Q}(t)^{2}-\expectationvalue{Q}_{t}^{2}}}{\sigma_{Q}(t)^{2}-2\expectationvalue{Q}_{t}^{2}}\absolutevalue{Q}\right).$
(78)

This result can also be obtained using the methods described in [33].

It is worth noting that the exponential term in the heat distribution is time-
independent since
$\frac{\expectationvalue{Q}_{t}}{\sigma_{Q}^{2}-2\expectationvalue{Q}_{t}^{2}}=\frac{k_{i}-k_{m}}{2k_{m}}$
does not depend on $t$. On the other hand, the argument of the Bessel function
$K_{0}$ is time-dependent. The ratio of the heat probability distributions
$P(Q,t)/P(-Q,t)$ satisfies a similar relation to the one found for the work in
Eq. (65)

$\frac{P_{Q}(Q,t)}{P_{Q}(-Q,t)}=\exp{\frac{2\expectationvalue{Q}_{t}}{{\sigma_{Q}^{2}}_{t}-2\expectationvalue{Q}_{t}^{2}}Q}=\exp{\frac{k_{i}-k_{m}}{k_{m}}Q}.$
(79)

A numerical test of this relation for other protocols in the range
$t_{f}\ll{\tau_{\text{relax}}}=1$ is shown in Fig. (7) revealing again its
universality. For longer protocols $t_{f}\sim 1$, the relation given by Eq.
(79), which has been proven for the TSP protocol, remains valid at large
values of $|Q|$ but shows some deviations at small values as illustrated in
Fig. 8 for several protocols. A similar relation to (79) where
$P_{Q}(Q,t)/P_{Q}(-Q,t)$ is a simple exponential function of $Q$ appears also
in a different context when one considers the relaxation from the
microcanonical to the canonical ensembles of this system [34].

Figure 7:  The protocol described by Martinez et al. [5], the protocol $k(t)$
described in Eq. (15), the linear (22) and the optimal protocols (20) exhibit
the same relation as the TSP given in Eq. (79). This observation is based on
$10^{6}$ simulations performed with $k_{i}=1/2$ and $t_{f}=1/30$. Figure 8:
The relation predict by (79) continues to be valid even for $t_{f}\sim 1$ as
is shown for the plot. The figure is based on $10^{6}$ simulations with
parameters $t_{f}=0.7$ and $k_{i}=1/2$. Figure 9:  The presented plot provides
a comparative analysis of the probability density function for heat,
specifically focusing on $t=t_{f}$. The analysis includes both the theoretical
result (78) obtained for the TSP and the results from $10^{6}$ simulations.
The simulations were conducted using a specific parameter setting, where
$k_{i}$ is half of $k_{f}$. Notably, a remarkable level of agreement is
observed between the theoretical and simulation outcomes, indicating a high
degree of concurrence.

####  III.2.2 For $t>t_{f}$

After the second step of the TSP, the average heat reaches its equilibrium
value given by Eq. (76) evaluated at $t=t_{f}$. However, surprisingly, moments
of order $n\geq 2$ of the heat distribution take longer to reach their
equilibrium values. To investigate this phenomenon, we attempt to calculate
the heat distribution for $t>t_{f}$, which necessitates the joint probability
distribution of the three-points:

$P_{Q}(Q,t)=\int
dx_{i}dx_{f}dx,\delta(Q-Q_{D})P_{3}(x_{i},t_{i};x_{f},t_{f};x_{t},t).$ (80)

Here, $P_{3}$ is a product of two Ornstein-Uhlenbeck processes with stiffness
$k_{m}$ and $k_{f}$ and the initial probability distribution,

$\displaystyle P_{3}(x_{i},t_{i};x_{f},t_{f};x_{t},t)=$ $\displaystyle
P(x_{i},t_{i})P(x_{f},t_{f}|x_{i},t_{i})$ $\displaystyle\times
P(x_{t},t|x_{f},t_{f}),$ (81)

with

$P(x_{t},t|x_{f},t_{f})=\frac{1}{\sqrt{2\pi\sigma_{f}(t)^{2}}}\exp{-\frac{(x_{t}-x_{f}e^{-k_{f}(t-t_{f})})^{2}}{2\sigma_{f}(t)^{2}}},$
(82)

where

$\sigma_{f}(t)^{2}=\frac{1-e^{-2k_{f}(t-t_{f})}}{k_{f}},$ (83)

and $P(x_{f},t_{f}|x_{i},t_{i})$ is given by (37).

By combining these expressions, utilizing the Fourier representation of the
delta distribution, and performing the resulting Gaussian integrals, we can
obtain the characteristic function.

$\hat{P}_{Q}(z)=\frac{1}{\sqrt{1+c_{1}z+c_{2}z^{2}+c_{3}z^{3}}}$ (84)

whose coefficients are given by

$\displaystyle c_{1}$ $\displaystyle=$ $\displaystyle
A(1-e^{-2k_{f}(t-t_{f})})-\frac{ik_{m}(k_{m}-k_{i})}{k_{i}}\sigma_{m}^{2}(t_{f}),$
$\displaystyle c_{2}$ $\displaystyle=$ $\displaystyle
B(1-e^{-2k_{f}(t-t_{f})})+\frac{k_{m}^{2}}{k_{i}}\sigma_{m}^{2}(t_{f}),$
$\displaystyle c_{3}$ $\displaystyle=$ $\displaystyle
C(1-e^{-2k_{f}(t-t_{f})}),$

where $A,B$ and $C$ are constants given by

$\displaystyle A$ $\displaystyle=$
$\displaystyle\frac{i(k_{i}-k_{f})}{k_{i}}+\frac{ik_{f}(k_{m}-k_{i})}{k_{i}}\sigma_{m}^{2}(t_{f}),$
$\displaystyle B$ $\displaystyle=$
$\displaystyle\frac{k_{f}+(k_{f}k_{i}-2k_{m}k_{f}-k_{i}k_{m}+k_{m}^{2})\sigma_{m}^{2}(t_{f})}{k_{i}},$
$\displaystyle C$ $\displaystyle=$
$\displaystyle\frac{i(k_{m}-k_{f})k_{m}\sigma_{m}^{2}(t_{f})}{k_{i}}.$

These coefficients are completely general, meaning that $k_{m}$ does not
necessarily satisfy the consistency equation (26). If $k_{m}$ does satisfy the
consistency equation, the coefficients $B$ and $C$ are nonzero, while the
coefficient $A$ vanishes. As a result, the average heat for $t>t_{f}$ becomes
time-independent. However, the variance and higher-order moments of the heat
depend on time and have a characteristic time given by $1/(2k_{f})$. If
$k_{m}$ is a solution of the consistency equation (26), the average heat and
variance are given by:

$\displaystyle\expectationvalue{Q}$ $\displaystyle=$
$\displaystyle-\frac{k_{m}}{2}\left(\frac{1}{k_{i}}-\frac{1}{k_{f}}\right),$
(85) $\displaystyle\sigma_{Q}^{2}$ $\displaystyle=$ $\displaystyle
2\expectationvalue{Q}^{2}-\frac{2k_{m}\expectationvalue{Q}}{k_{m}-k_{i}}+D(1-e^{-2k_{f}(t-t_{f})}),$
(86)

here $D$ is

$D=\frac{(k_{f}-k_{m})(k_{i}^{2}+k_{f}k_{m}-k_{i}k_{m})}{k_{i}k_{f}(k_{i}-k_{m})}.$
(87)

Fig. 10 shows the evolution of the average heat and its standard deviation.
When $t>t_{f}$, the average heat has stabilized but its standard deviation
continues to change in time with a relaxation time of order $1/k_{f}$. This
finding holds considerable importance as it indicates that the Brownian
particle reaches a state of equilibrium in terms of its position and work
probability distribution by the final time $t_{f}$. However, it is worth
noting that the heat distribution requires additional time to reach its
ultimate state of equilibrium. The primary aim of the shortcut protocols is to
achieve equilibrium in the position distribution, while the attainment of
equilibrium in other relevant distributions is not necessarily simultaneous
with that of the position distribution.

Figure 10:  Standard deviation and average of heat for the TSP with parameters
$t_{f}=1/30$ and $k_{i}=1/2$. The dotted line is the asymptotic value of
$\sigma_{Q}$ when $t\rightarrow\infty$.

##  IV Entropy production

In this section, we specifically focus on the entropy production associated
with the two-step protocol. This quantity holds significant importance in
understanding the inherent irreversibility within nonequilibrium systems and
provides valuable insights into the fundamental thermodynamic principles that
govern the TSP.

We recall the definition of total entropy in the realm of stochastic
thermodynamics [19, 17]:

$\Delta S=-\frac{Q}{T}+k_{B}\ln{\frac{P_{i}}{P_{f}}},$ (88)

where $\Delta S$ is the entropy change between two states with the probability
density function (PDF) of the position given by $P_{i}$ and $P_{f}$,
respectively. The first term corresponds to the increase in entropy of the
environment, which is assumed to be in equilibrium at temperature $T$ and the
last term represents the increase in entropy of the system. In our notation,
we will use $\Delta S$ and $\Sigma$ interchangeably to represent entropy
production.

By employing a methodology analogous to that used for calculating the
probability density functions of heat and work, we can derive the PDF of
entropy production. Due to the similarities with the previous derivations, we
will omit several intermediate steps. Similar to before, the PDF of entropy
production can be expressed as:

$P_{\Sigma}(\Sigma)=\frac{1}{2\pi}\int_{-\infty}^{\infty}dz\hat{P}_{\Sigma}(z)\exp{(iz(\Sigma-\Sigma_{0}))},$
(89)

where $\hat{P}_{\Sigma}(z)e^{-iz\Sigma_{0}}$ is the characteristic function
associated with the entropy production random variable and

$\Sigma_{0}=\frac{1}{2}\ln{\frac{\sigma_{X}(t)^{2}}{\sigma_{i}^{2}}}.$ (90)

Once again, it is necessary to consider two distinct cases:

####  IV.0.1 For $0\leq t\leq t_{f}$

To calculate the second term in the entropy expression, we need the
probability density function of the position at the initial time as well as
the position PDF at $t\leq t_{f}$, as given by Eq. (25). By utilizing the
expression for the heat done, as shown in Eq. (75), and the two-point
probability distribution given in Eq. (36), we can determine the total entropy
(in the chosen units) for $t\leq t_{f}$ as follows:

$\Sigma=\frac{1}{2}k_{m}(x_{i}^{2}-x_{t}^{2})+\Sigma_{0}-\frac{x_{i}^{2}}{2\sigma_{i}^{2}}+\frac{x_{t}^{2}}{2\sigma_{X}(t)^{2}}.$
(91)

The function $\hat{P}_{\Sigma}(z)$ is

$\hat{P}_{\Sigma}(z)=\frac{1}{\sqrt{1+d_{1}z+d_{2}z^{2}}},$ (92)

where

$\displaystyle d_{1}$ $\displaystyle=$
$\displaystyle\frac{i\left(-k_{m}\sigma_{m}^{4}k_{i}+k_{m}^{2}\sigma_{m}^{4}-k_{m}\sigma_{m}^{2}+1\right)}{\sigma_{m}^{2}k_{i}},$
(93) $\displaystyle d_{2}$ $\displaystyle=$
$\displaystyle\frac{\left(k_{m}-k_{i}\right)\left(k_{m}\sigma_{m}^{2}-1\right)}{k_{i}}.$
(94)

Computing the integral (89), we arrive to the expression

$\displaystyle P_{\Sigma}(\Sigma,t)$
$\displaystyle=\frac{1}{\pi\sqrt{\sigma_{\Sigma}(t)^{2}-2\expectationvalue{\Sigma-\Sigma_{0}}_{t}^{2}}}\exp{\frac{\expectationvalue{\Sigma-\Sigma_{0}}_{t}}{\sigma_{\Sigma}(t)^{2}-2\expectationvalue{\Sigma-\Sigma_{0}}_{t}^{2}}(\Sigma-\Sigma_{0})}K_{0}\left(\frac{\sqrt{\sigma_{\Sigma}(t)^{2}-\expectationvalue{\Sigma-\Sigma_{0}}_{t}^{2}}}{\sigma_{\Sigma}(t)^{2}-2\expectationvalue{\Sigma-\Sigma_{0}}_{t}^{2}}\absolutevalue{\Sigma-\Sigma_{0}}\right).$
(95)

where

$\displaystyle\expectationvalue{\Sigma}_{t}$ $\displaystyle=$
$\displaystyle\Sigma_{0}+\frac{k_{m}-k_{i}}{2k_{i}}\left(1-e^{-2k_{m}t}\right),$
(96) $\displaystyle\sigma_{\Sigma}^{2}(t)$ $\displaystyle=$ $\displaystyle
2\expectationvalue{\Sigma-\Sigma_{0}}^{2}+\frac{(k_{i}-k_{m})^{2}\left(1-e^{-2k_{m}t}\right)}{k_{i}\left(k_{m}-k_{i}\left(1-e^{2k_{m}t}\right)\right)}.$
(97)

####  IV.0.2 For $t\geq t_{f}$

In this case, the entropy is given by

$\Sigma=\frac{1}{2}k_{m}(x_{i}^{2}-x_{f}^{2})+\frac{1}{2}k_{f}(x_{f}^{2}-x_{t}^{2})+\Sigma_{0}-\frac{x_{i}^{2}}{2\sigma_{i}^{2}}+\frac{x_{t}^{2}}{2\sigma_{X}(t)^{2}}.$
(98)

It is worth noting that if we ensure a final equilibrium situation at $t_{f}$
by imposing the consistency relation (26), there is a cancellation between
$-k_{f}x_{t}^{2}/2$ and $x_{t}^{2}/(2\sigma_{X}(t)^{2})$ making $\Sigma$
stationary.

In this time interval and assuming that $k_{m}$ is general, the function
$\hat{P}_{\Sigma}(z)$ takes the form

$\hat{P}_{\Sigma}(z)=\frac{1}{\sqrt{1+f_{1}z+f_{2}z^{2}+f_{3}z^{3}}},$ (99)

where the coefficients $f_{1},f_{2},$ and $f_{3}$ are complex expressions of
the system parameters (refer to Appendix B for details). Due to the cubic
polynomial inside the square root, the integral in Eq. (89) cannot be
evaluated in closed form.

Therefore, in this section, we will impose a constraint on the variable
$k_{m}$ to satisfy the consistency equation (26). Under this constraint, the
coefficient $f_{3}$ vanishes, allowing us to obtain a closed-form expression.
This is similar to the result in the previous section, Eq. (95), although the
average and variance will exhibit modifications. Specifically, they are time-
independent, given by:

$\displaystyle\expectationvalue{\Sigma}$ $\displaystyle=$
$\displaystyle\Sigma_{0}+\frac{k_{m}}{2}\left(\frac{1}{k_{i}}-\frac{1}{k_{f}}\right),$
(100) $\displaystyle\sigma_{\Sigma}^{2}$ $\displaystyle=$ $\displaystyle
2\expectationvalue{\Sigma-\Sigma_{0}}^{2}+\frac{2\expectationvalue{\Sigma-\Sigma_{0}}(k_{m}-k_{f})}{k_{m}}.$
(101)

with $\Sigma_{0}=\frac{1}{2}\ln(k_{i}/k_{f})$. These computed values coincide
with the final values obtained in the previous time interval, Eqs. (96) and
(97) for $t=t_{f}$. Furthermore, for $t>t_{f}$, these values remain constant
as expected for an equilibrium situation. Fig. 11 shows the evolution of the
average entropy production and its standard deviation.

The entropy production distribution satisfies a relation similar to the ones
for work and heat shown in Eqs. (65) and (79) if we define
$\tilde{\Sigma}=\Sigma-\Sigma_{0}$. This is given by

$\displaystyle\frac{P_{\tilde{\Sigma}}(\tilde{\Sigma})}{P_{\tilde{\Sigma}}(-\tilde{\Sigma})}$
$\displaystyle=$
$\displaystyle\exp{\frac{2\expectationvalue{\Sigma-\Sigma_{0}}_{t}}{\sigma_{\Sigma}^{2}-2\expectationvalue{\Sigma-\Sigma_{0}}_{t}^{2}}\tilde{\Sigma}}$
(102) $\displaystyle=$
$\displaystyle\exp{\frac{k_{m}}{k_{m}-k_{f}}\tilde{\Sigma}}\,.$

Fig. 12 represents this relation for $t=t_{f}\ll 1$ and the simulation results
for the other protocols considered in this paper. This figure shows the
universality of the relation (102) when $t_{f}\ll 1$. However this
universality breaks down for longer protocols $t_{f}\sim 1$ as can be
appreciated at Fig. 13.

Figure 11:  Standard deviation and average for the entropy production for the
TSP protocol with parameters $k_{i}=1/2$ and $t_{f}=1/30$. Figure 12:
Relation (102) for several protocols. $10^{6}$ simulations have been performed
with $k_{i}=1/2$ and $t_{f}=1/30$. Figure 13:  The prediction given by (102)
fails when $t_{f}\sim 1$; however, the behavior looks linear. This plot is for
$10^{6}$ simulations with parameters $k_{i}=1/2$ and $t_{f}=0.7$.

####  IV.0.3 Fluctuation theorem for the entropy production

The entropy distribution must satisfy a fluctuation theorem [20, 21], in the
same way that work distribution must satisfy the Jarzynski equality and the
Crocks relation. Hence, we will prove the integral fluctuation relation for
the entropy production. Using the definition of the expectation value and
inserting the integral representation given in Eq. (89), we can obtain the
following identity:

$\langle e^{-\Sigma}\rangle=\hat{P}_{\Sigma}(-i)e^{-\Sigma_{0}}.$ (103)

For arbitrary $k_{m}$ (not necessarily a solution of the consistency equation
(26)), we can compute $\hat{P}_{\Sigma}(-i)$ as:

$\hat{P}(-i)=\frac{\sigma_{X}(t)}{\sigma_{i}}.$ (104)

As a result, the expectation value is given by:

$\langle e^{-\Sigma}\rangle=1.$ (105)

This is the well-known integral fluctuation theorem for the entropy production
[20, 21, 22].

##  V Conclusion

We have reviewed how to construct ESE protocols, in particular, we have shown
that for a Brownian particle in a harmonic potential with time-dependent
stiffness and at equilibrium initially with a Gaussian distribution, the
position distribution remains Gaussian for later times and the variance of the
position distribution depends on the protocol. Moreover, we have thoroughly
described the two-step protocol and examined its energetic behaviors. Notably,
we obtained analytical expressions for the distribution functions of work,
heat, and entropy production. These analytical solutions comply with the
fluctuation theorems and identities well-established in the existing
literature [20, 21, 22, 18, 19, 16, 17].

The significance of the two-step protocol extends beyond analytical solutions,
as it also characterizes more general protocols within the time range of
interest of the shortcut protocols (i.e., when
$t_{f}\ll{\tau_{\text{relax}}}$). Our findings unveil intriguing insights into
the behavior of various protocols that are of common interest for theoretical
and experimental investigations. The average value of the stiffness $k_{m}$ is
inversely proportional to the protocol duration $t_{f}$, explicitly given by
Eq. (27). This relation can serve as a guide for designing fast thermalization
protocols. It also imposes some restrictions on them. If for some practical
application the stiffness has to be bounded, Eq. (27) imposes a limit on how
short the protocol can be [23]. Regarding the energetics statistics, we
obtained strong numerical evidence that suggest that the probability
distribution functions of work, heat and entropy production satisfy universal
simple relations given by Eqs. (65), (79) and (102) when we compare the
positive versus the negative production of these quantities.

In summary, our research significantly advances the comprehension of
stochastic thermodynamics in colloidal Brownian particles. We have gained
valuable insights into the dynamics of overdamped Brownian particles under
time-dependent potentials and external control, facilitating the acceleration
of equilibration times. These findings hold great potential for the design and
application of efficient protocols across diverse scientific and engineering
disciplines.

Despite the valuable insights gained, this study faced certain limitations,
particularly in the context of the overdamped regime. Future research may
involve finding analytical solutions in the underdamped regime and extending
the analysis to systems of interacting particles. By addressing these aspects,
future investigations can provide a more comprehensive understanding of the
fast thermalization protocols for Brownian particles in complex environments,
further enriching the field of stochastic thermodynamics.

##  VI Acknowledgments

This work was supported by Fondo de Investigaciones, Facultad de Ciencias,
Universidad de los Andes INV-2021-128-2267. We thank Emmanuel Trizac for
fruitful discussions.

##  VII Appendixes

## Appendix A

In order to prove the Crooks identity (70), it is necessary to demonstrate the
symmetry under the permutation $k_{i}\leftrightarrow k_{f}$ of the Bessel
function in (61). To achieve this goal, we use the followings definitions

$\displaystyle S$ $\displaystyle=$
$\displaystyle\frac{(k_{m}-k_{f})(k_{m}-k_{i})\sigma_{m}^{2}(t)}{2},$ (106)
$\displaystyle A$ $\displaystyle=$ $\displaystyle\frac{k_{f}-k_{i}}{2},$ (107)

where $S$ and $A$ stand for symmetric and anti-symmetric under the
permutation, respectively. The equations (67) takes the form

$\displaystyle\expectationvalue{W}$ $\displaystyle=$
$\displaystyle\frac{A+S}{k_{i}},$ (108) $\displaystyle\sigma_{W}^{2}$
$\displaystyle=$ $\displaystyle\frac{2(A+S)^{2}+2k_{i}S}{k_{i}^{2}}.$ (109)

Now, the combination in the argument of the Bessel function is

$\frac{\sqrt{\sigma_{W}^{2}-\expectationvalue{W}^{2}}}{\sigma_{W}^{2}-2\expectationvalue{W}^{2}}=\frac{\sqrt{A^{2}+S^{2}+2S(A+k_{i})}}{2S}.$
(110)

The only term that is not evident symmetric is the combination $A+k_{i}$,
however

$A+k_{i}=\frac{k_{i}+k_{f}}{2}$ (111)

that is symmetric under the permutation.

## Appendix B

In this appendix, where give the intricate expression of the coefficients
given intervening in Eq. (99):

$\displaystyle\sigma_{t}^{2}f_{3}$ $\displaystyle=$ $\displaystyle
i\sigma_{f}^{2}\sigma_{m}^{2}\left(k_{f}-k_{m}\right)\left(k_{f}\sigma_{t}^{2}-1\right)\left(\sigma_{i}^{2}k_{m}-1\right),$
$\displaystyle\sigma_{t}^{2}f_{2}$ $\displaystyle=$
$\displaystyle\sigma_{m}^{2}\left(\sigma_{i}^{2}k_{m}-1\right)\left(k_{m}\sigma_{t}^{2}-1\right)-\sigma_{f}^{2}\left(k_{f}\sigma_{t}^{2}-1\right)$
$\displaystyle\times\left(\sigma_{m}^{2}\left(2k_{f}-k_{m}\right)\left(\sigma_{i}^{2}k_{m}-1\right)-k_{f}\sigma_{i}^{2}+1\right),$
$\displaystyle\sigma_{t}^{2}f_{1}$ $\displaystyle=$
$\displaystyle-i\sigma_{f}^{2}\left(k_{f}\sigma_{t}^{2}-1\right)\left(k_{f}\left(\sigma_{m}^{2}\left(\sigma_{i}^{2}k_{m}-1\right)-\sigma_{i}^{2}\right)+1\right)+$
$\displaystyle+$ $\displaystyle
i\sigma_{i}^{2}\left(k_{m}\sigma_{m}^{2}\left(k_{m}\sigma_{t}^{2}-1\right)+1\right)-i\sigma_{t}^{2}\left(k_{m}\sigma_{m}^{2}+1\right)+$
$\displaystyle+$ $\displaystyle i\sigma_{m}^{2},$

where $\sigma_{t}$ is the solution of (11) in the interval $t\geq t_{f}$ given
by

$\sigma^{2}_{t}=\left(-\frac{1}{k_{f}}+\frac{1-k_{m}\sigma_{m}^{2}}{k_{i}}+\sigma_{m}^{2}\right)e^{-2k_{f}\left(t-t_{f}\right)}+\frac{1}{k_{f}}.$
(112)

Using the consistency equation, these coefficients collapse to

$\displaystyle f_{3}$ $\displaystyle=$ $\displaystyle 0,$ (113) $\displaystyle
f_{2}$ $\displaystyle=$
$\displaystyle\frac{\left(k_{i}-k_{f}\right)\left(k_{f}-k_{m}\right)}{k_{f}k_{i}},$
(114) $\displaystyle f_{1}$ $\displaystyle=$ $\displaystyle-
ik_{m}\left(\frac{1}{k_{f}}-\frac{1}{k_{i}}\right).$ (115)

Because $f_{3}$ is zero, the polynomial (99) is of second order and the
inverse Fourier transform can be computed in terms of Bessel functions (Eq.
(95)).

## References

  * Toda  _et al._ [2012] M. Toda, R. Kubo, R. Kubo, M. Toda, N. Saito, N. Hashitsume, and N. Hash tsume,  _Statistical Physics II: Nonequilibrium Statistical Mechanics_, Springer Series in Solid-State Sciences (Springer Berlin Heidelberg, 2012).
  * Unanyan  _et al._ [1997] R. Unanyan, L. Yatsenko, K. Bergmann, and B. Shore, Laser-induced adiabatic atomic reorientation with control of diabatic losses, Optics Communications 139, 48 (1997).
  * Chen  _et al._ [2010a] X. Chen, A. Ruschhaupt, S. Schmidt, A. del Campo, D. Guéry-Odelin, and J. G. Muga, Fast optimal frictionless atom cooling in harmonic traps: Shortcut to adiabaticity, Phys. Rev. Lett. 104, 063002 (2010a).
  * Guéry-Odelin  _et al._ [2019] D. Guéry-Odelin, A. Ruschhaupt, A. Kiely, E. Torrontegui, S. Martínez-Garaot, and J. Muga, Shortcuts to adiabaticity: Concepts, methods, and applications, Reviews of Modern Physics 91 (2019).
  * Martínez  _et al._ [2016] I. Martínez, A. Petrosyan, D. Guéry-Odelin, E. Trizac, and S. Ciliberto, Engineered swift equilibration of a brownian particle, Nature Physics 12, 843–846 (2016).
  * Li  _et al._ [2017] G. Li, H. T. Quan, and Z. C. Tu, Shortcuts to isothermality and nonequilibrium work relations, Phys. Rev. E 96, 012144 (2017).
  * Guéry-Odelin  _et al._ [2023] D. Guéry-Odelin, C. Jarzynski, C. A. Plata, A. Prados, and E. Trizac, Driving rapidly while remaining in control: Classical shortcuts from Hamiltonian to stochastic dynamics, Reports on Progress in Physics 86, 035902 (2023).
  * Chen  _et al._ [2010b] X. Chen, A. Ruschhaupt, S. Schmidt, A. del Campo, D. Guéry-Odelin, and J. G. Muga, Fast optimal frictionless atom cooling in harmonic traps: Shortcut to adiabaticity, Phys. Rev. Lett. 104, 063002 (2010b).
  * Chen  _et al._ [2010c] X. Chen, I. Lizuain, A. Ruschhaupt, D. Guéry-Odelin, and J. G. Muga, Shortcut to adiabatic passage in two- and three-level atoms, Phys. Rev. Lett. 105, 123003 (2010c).
  * Cunuder  _et al._ [2016] A. Cunuder, I. Martínez, A. Petrosyan, D. Guéry-Odelin, E. Trizac, and S. Ciliberto, Fast equilibrium switch of a micro mechanical oscillator, Applied Physics Letters 109 (2016).
  * Wang  _et al._ [2002] G. M. Wang, E. M. Sevick, E. Mittag, D. J. Searles, and D. J. Evans, Experimental Demonstration of Violations of the Second Law of Thermodynamics for Small Systems and Short Time Scales, Physical Review Letters 89, 050601 (2002).
  * Blickle  _et al._ [2006] V. Blickle, T. Speck, L. Helden, U. Seifert, and C. Bechinger, Thermodynamics of a Colloidal Particle in a Time-Dependent Nonharmonic Potential, Physical Review Letters 96, 070603 (2006).
  * Blickle and Bechinger [2012] V. Blickle and C. Bechinger, Realization of a micrometre-sized stochastic heat engine, Nature Physics 8, 143 (2012).
  * Martínez  _et al._ [2016] I. A. Martínez, E. Roldán, L. Dinis, D. Petrov, J. M. Parrondo, and R. A. Rica, Brownian Carnot engine, Nature Physics 12, 67 (2016), 1412.1282 .
  * Speck [2011] T. Speck, Work distribution for the driven harmonic oscillator with time-dependent strength: exact solution and slow driving, Journal of Physics A: Mathematical and Theoretical 44, 305001 (2011).
  * Jarzynski [1997a] C. Jarzynski, Nonequilibrium equality for free energy differences, Phys. Rev. Lett. 78, 2690 (1997a).
  * Jarzynski [1997b] C. Jarzynski, Equilibrium free-energy differences from nonequilibrium measurements: A master-equation approach, Physical Review E 56, 5018 (1997b).
  * Crooks [1998] G. E. Crooks, Nonequilibrium Measurements of Free Energy Differences for Microscopically Reversible Markovian Systems, Journal of Statistical Physics 90, 1481 (1998).
  * Crooks [1999] G. E. Crooks, Entropy production fluctuation theorem and the nonequilibrium work relation for free energy differences, Physical Review E 60, 2721 (1999).
  * Evans  _et al._ [1993] D. J. Evans, E. G. D. Cohen, and G. P. Morriss, Probability of second law violations in shearing steady states, Phys. Rev. Lett. 71, 2401 (1993).
  * Evans and Searles [1994] D. J. Evans and D. J. Searles, Equilibrium microstates which generate second law violating steady states, Phys. Rev. E 50, 1645 (1994).
  * Gallavotti and Cohen [1995] G. Gallavotti and E. G. D. Cohen, Dynamical ensembles in stationary states, Journal of Statistical Physics 80, 931 (1995).
  * Plata  _et al._ [2019] C. A. Plata, D. Guéry-Odelin, E. Trizac, and A. Prados, Optimal work in a harmonic trap with bounded stiffness, Phys. Rev. E 99, 012140 (2019).
  * Bechhoefer [2021] J. Bechhoefer,  _Control Theory for Physicists_ (Cambridge University Press, 2021).
  * Risken and Haken [1989] H. Risken and H. Haken,  _The Fokker-Planck Equation: Methods of Solution and Applications Second Edition_ (Springer, 1989).
  * Sekimoto [2010] K. Sekimoto,  _Stochastic Energetics_, Lecture Notes in Physics (Springer Berlin Heidelberg, 2010).
  * Peliti and Pigolotti [2021] L. Peliti and S. Pigolotti,  _Stochastic Thermodynamics: An Introduction_ (Princeton University Press, 2021).
  * Schmiedl and Seifert [2007] T. Schmiedl and U. Seifert, Efficiency at maximum power: An analytically solvable model for stochastic heat engines, Europhysics Letters 81, 20003 (2007).
  * S.  _et al._ [1964] P. L. S., B. V. G., G. R. V., and M. E. F.,  _The Mathematical Theory of Optimal Processes_ (New York: The Macmillian Co., 1964).
  * Bayati and Trizac [2021] P. Bayati and E. Trizac, Diffusiophoresis driven colloidal manipulation and shortcuts to adiabaticity, New Journal of Physics 23, 063028 (2021).
  * Gardiner [2004] C. W. Gardiner,  _Handbook of stochastic methods for physics, chemistry and the natural sciences_ , 3rd ed., Springer Series in Synergetics, Vol. 13 (Springer-Verlag, Berlin, 2004) pp. xviii+415.
  * Téllez [2022] G. Téllez, Langesim, https://github.com/GabrielTellez/langesim (2022).
  * Imparato  _et al._ [2007] A. Imparato, L. Peliti, G. Pesce, G. Rusciano, and A. Sasso, Work and heat probability distribution of an optically driven brownian particle: Theory and experiments, Physical Review. E, Statistical, nonlinear, and soft matter physics 76, 050101 (2007).
  * Ardila  _et al._ [2023] L. F. Ardila, N. Torres, J. D. Muñoz, and C. Viviescas, Heat exchange fluctuation relation for the transition from a micro-canonical to a canonical ensemble in a classical harmonic oscillator, STATPHYS28 Book of abstracts (2023).