Physica A 562 (2021) 125277

Contents lists available at ScienceDirect

Physica A
journal homepage: www.elsevier.com/locate/physa

Irreversibility and entropy production of a thermally driven
micromachine
Isamu Sou, Yuto Hosaka, Kento Yasuda, Shigeyuki Komura∗

Department of Chemistry, Graduate School of Science, Tokyo Metropolitan University, Tokyo 192-0397, Japan

articleinfo

Article history:
Received 1 August 2020
Received in revised form 12 September 2020
Available online 18 September 2020a b s t r a c t
We discuss the non-equilibrium properties of a thermally driven micromachine consisting of three spheres which are in equilibrium with independent heat baths characterized
by different temperatures. Within the framework of a linear stochastic Langevin description, we calculate the time-dependent average irreversibility that takes a maximum
value for a finite time. This time scale is roughly set by the spring relaxation time. The
steady-state average entropy production rate is obtained in terms of the temperatures
and the friction coefficients of the spheres. The average entropy production rate depends
on thermal and/or mechanical asymmetry of a three-sphere micromachine. We also
obtain the center of mass diffusion coefficient of a thermally driven three-spheremicromachine as a function of different temperatures and friction coefficients. With
the results of the total entropy production rate and the diffusion coefficient, we finallydiscuss the efficiency of a thermally driven micromachine.
© 2020 Elsevier B.V. All rights reserved.

1. Introduction
Microswimmers are tiny machines that swim in a fluid and they are expected to be used in microfluidics andmicrosystems [1]. By transforming chemical energy into mechanical energy, microswimmers change their shape and moveefficiently in viscous environments. According to Purcell’s scallop theorem, reciprocal body motion cannot be used forlocomotion in a Newtonian fluid [2,3]. As one of the simplest models showing non-reciprocal body motion, Najafi andGolestanian proposed a three-sphere swimmer [4,5], in which three in-line spheres are linked by two arms of varying
length. Such a swimmer has been experimentally realized by using colloidal beads manipulated by optical tweezers [6]or by controlling ferromagnetic particles at an air–water interface [7,8].
Recently, the present authors have proposed a generalized three-sphere microswimmer model in which the spheresare connected by two harmonic springs, i.e., an elastic microswimmer [9,10]. Later, our model was further extended
to a thermally driven elastic microswimmer [11], suggesting a new mechanism for locomotion that is purely induced bythermal fluctuations without any external forcing. The key setting of the model is that the three spheres are in equilibriumwith independent heat baths characterized by different temperatures (as described later in Fig. 1). For this thermally driven
three-sphere micromachine, the average velocity was calculated to be [11]⟨V ⟩ =kB (T3 − T1 )96π ηℓ2,(1)

where kB is the Boltzmann constant, T1 and T3 are the temperatures of the first and the third spheres (see Fig. 1), η is theviscosity of the surrounding fluid, and ℓ is the natural length of the two springs. We have shown that a combination of∗ Corresponding author.
E-mail address: komura@tmu.ac.jp (S. Komura).
https://doi.org/10.1016/j.physa.2020.125277
0378-4371/© 2020 Elsevier B.V. All rights reserved.

I. Sou, Y. Hosaka, K. Yasuda et al.

Physica A 562 (2021) 125277

heat transfer and hydrodynamic interactions among the spheres leads to directional locomotion in a steady-state, which
can be described in terms of ‘‘stochastic energetics’’ [12–14].
Systems in thermodynamic equilibrium obey detailed balance meaning that transition rates between any two microscopic states are pairwise balanced [15,16]. For non-equilibrium steady-state situations, however, detailed balance isbroken and a probability flux loop exists in a configuration phase space [17–20]. Recently, the present authors discussed
the non-equilibrium steady-state probability distribution function of a thermally driven three-sphere micromachineand calculated its probability flux in the corresponding configuration space [21]. The resulting probability flux can be
expressed in terms of a frequency matrix to characterize a non-equilibrium steady-state [22,23]. Importantly, we haveobtained a linear relation between the eigenvalue of the frequency matrix and the average velocity of a thermally driven
micromachine [21].
Since our model micromachine offers a new type of thermal ratchet [11], a more detailed analysis based on nonequilibrium statistical mechanics is required in order to elucidate the physical mechanism for the locomotion of athermally driven three-sphere micromachine. One of the important quantities to measure the statistical irreversibilityof a non-equilibrium process is the entropy production rate [24–26]. In this paper, within the framework of a linear
stochastic Langevin description, we calculate the time-dependent average irreversibility of a thermally driven threesphere micromachine in the absence of hydrodynamic interactions. The average irreversibility is important because its
initial growth rate gives the average entropy production rate. Using this fact, we explicitly obtain the average entropy
production rate for a thermally driven three-sphere micromachine, and express it in terms of different temperatures andfriction coefficients. We examine in detail how the entropy production rate depends on the asymmetry of the temperatures
and/or the friction coefficients of a three-sphere micromachine.
Although a micromachine does not exhibit any directional motion in the absence of hydrodynamic interactions [21], it
undergoes a thermal Brownian motion. In addition to the above mentioned non-equilibrium quantities, we also calculate
the center of mass diffusion coefficient of a thermally driven three-sphere micromachine (even though we neglect theinertia of the spheres). Using the results of the diffusion coefficient for a dimer and trimer, we predict a simple anduseful expression for the center of mass diffusion coefficient of an elastic n-sphere micromachine. We consider that our
expression is useful because the center of mass diffusion is easier to measure than the center of friction diffusion [27].
Finally, we obtain the total entropy production rate including the center of mass motion, and further discuss the efficiency
of a thermally driven three-sphere micromachine. Our result shows that the efficiency becomes larger for specificcombinations of the temperatures.
In the next section, we briefly review the framework of a linear Langevin model and describe how the average
irreversibility and the entropy production rate are obtained in general. In Section 3, we explain our model of a thermally
driven three-sphere micromachine by introducing the coupled Langevin equations for the two spring lengths [21]. In
Sections 4 and 5, we explicitly calculate the average irreversibility and the average entropy production rate, respectively,
for a thermally driven three-sphere micromachine. The center of mass diffusion coefficient of a thermally driven
micromachine is given in Section 6. In Section 7, after calculating the total entropy production rate, we discuss theefficiency of a micromachine. Finally, a summary of our work and some discussion are given in Section 8.
2. Linear Langevin system
In this section, we briefly review Ref. [23] and pick up the important results for our calculation. Let us consider a linear
stochastic Langevin equation given by [28]
dr(t)
dt

= Ar(t) + Fξ (t),(2)

where r is the N-dimensional state vector of real numbers, A is an N × N real matrix representing the linear deterministic
dynamics, F is an N × N real matrix representing the noise forcing. We require that all eigenvalues of A have a negative realpart so that the system will eventually reach a steady-state. In Eq. (2), ξ is N-dimensional Gaussian white noise satisfying
the statistical properties⟨ξ (t)⟩ = 0,(3)

⟨ξ (t)ξ (t ) ⟩ = Bδ (t − t ),
′ T

′

(4)

where B is an N × N matrix representing the variance of noise, and superscript T represents the transpose. Then the
diffusion matrix D is obtained by
1

FBFT .
(5)
2
Using the above Langevin model, we discuss the probability of observing trajectory segments. For any two states r 0
and r 1 , we consider the trajectory probability p(r 0 , r 1 , t) as the probability of finding a trajectory segment within thelong trajectory which begins at r 0 and ends at r 1 a time t later. Such a trajectory probability can be expressed asD=p(r 0 , r 1 , t) = p(r 1 , t |r 0 )p0 (r 0 ),(6)2

I. Sou, Y. Hosaka, K. Yasuda et al.

Physica A 562 (2021) 125277

Fig. 1. Thermally driven elastic three-sphere micromachine. Three spheres are connected by two identical harmonic springs characterized by theelastic constant K and the natural length ℓ. The time-dependent positions of the spheres are denoted by xi (t) (i = 1, 2, 3) in a one-dimensional
coordinate system, and ζi is the friction coefficient for ith sphere. The three spheres are in equilibrium with independent heat baths having different
temperatures Ti .

where p(r 1 , t |r 0 ) is the transition probability of finding the system in state r 1 conditioned on the system being in stater 0 a time t earlier, and p0 (r 0 ) is the steady-state probability of finding the system in state r 0 . Similarly, the time-reversed
trajectory segment, one starting at r 1 and ending at r 0 , has a probabilityp(r 1 , r 0 , t) = p(r 0 , t |r 1 )p0 (r 1 ).

(7)

The irreversibility σ (r 0 , r 1 , t) of a trajectory segment with initial state r 0 and final state r 1 is defined by [23]σ (r 0 , r 1 , t) = lnp(r 0 , r 1 , t)p(r 1 , r 0 , t).

(8)

The system is reversible when σ = 0, while forward and reverse trajectories are distinguishable when σ ̸ = 0. By
introducing the probability P(σ ) of finding a trajectory segment with irreversibility σ , the fluctuation theorem can be
expressed as [23–26]P(σ )P(−σ )= eσ .

(9)

Although this fluctuation theorem gives a constraint on P(σ ), it does not completely fix its functional form. How to obtain
P(σ ) for a linear Langevin system is separately explained in Appendix A.
Because the model in Eq. (2) is linear with additive Gaussian white noise, the probabilities in Eqs. (6) and (7) are also
Gaussian, and they can be written in terms of the covariance of the dynamics. Weiss showed that the time-dependent
average irreversibility ⟨σ (t)⟩ is given by [23]

[
]
1
2At
⟨σ (t)⟩ = Tr C0 C−)−I ,t (I − e(10)where I is the N × N identity matrix. In the above, C0 is the steady-state covariance matrix satisfying the Lyapunov
equation [29]AC0 + C0 AT + 2D = 0,(11)which can be regarded as the fluctuation–dissipation relation [22]. Moreover, Ct is the time-dependent covariance matrixof the transition probability and is given by∫ t
Ct = 2TT

ds eA(t −s) DeA (t −s) = C0 − eAt C0 eA t .

(12)

0

Notice that C0 and Ct are related by C0 = limt →∞ Ct .
In Ref. [23], it was further shown that the steady-state average entropy production rate ⟨σ̇ ⟩ is given by the zero-time
growth rate of the average irreversibility in Eq. (10), i.e.,⟨σ̇ ⟩ =⏐⏐⏐d⟨σ (t)⟩ ⏐dt= Tr [AG] .

(13)

t =0

Here, G is the dimensionless gain matrix defined by
G = −(AC0 D−1 + I).

(14)

The gain matrix G is directly related to the violation of detailed balance [22], and the product AG in Eq. (13) measures
the noise amplification per unit time [23].
It is worth mentioning that the same average entropy production rate in Eq. (13) can also be obtained directly from thesteady-state probability flux [23,30]. In fact, the latter approach is more standard [31]. For a three-sphere micromachine,
the steady-state probability flux was calculated in our previous paper [21].
3

I. Sou, Y. Hosaka, K. Yasuda et al.

Physica A 562 (2021) 125277

3. Thermally driven three-sphere micromachine
In this section, we explain the model of a thermally driven elastic micromachine that has been introduced in ourprevious studies [11,21]. As schematically shown in Fig. 1, the model consists of three hard spheres connected by two
harmonic springs. We assume that the two springs are identical, and the common spring constant and the natural length
are given by K and ℓ, respectively. The positions of the three spheres in a one-dimensional coordinate system is defined
as xi (t) (i = 1, 2, 3).
Most importantly, we consider a situation where the three spheres are in thermal equilibrium with independent heatbaths at temperatures Ti [11,21]. When these temperatures are different, the system is driven out of equilibrium because
heat flux from a hotter sphere to a colder one is generated. The Langevin equations of motion of the three spheres are
given by
dx1
dt
dx2
dt
dx3
dt

=

K

ζ1

=−
=−

(

(x2 − x1 − ℓ) +
K

ζ2
K

ζ3

2T1

)1/2

ζ1

(x2 − x1 − ℓ) +(x3 − x2 − ℓ) +

K

ζ2
(

ξ1 ,

(x3 − x2 − ℓ) +

2T3

)1/2

ζ3

(15)

(

2T2

)1/2

ζ2

ξ2 ,

ξ3 ,

(16)
(17)

where ζi is the friction coefficient for ith sphere, and the Boltzmann constant kB is set to unity hereafter (except later in
Section 7). As we discuss later, the friction coefficient is generally proportional to the size of the sphere. Furthermore,ξi (t) is a zero mean and unit variance Gaussian white noise, independent for all the spheres:⟨ξi (t)⟩ = 0,(18)

⟨ξi (t)ξj (t )⟩ = δij δ (t − t ).
′

′

(19)

In contrast to Ref. [11], we do not consider hydrodynamic interactions acting between the different spheres. Although
the locomotion of a micromachine is not discussed in this paper (except in Section 7), we use a term ‘‘micromachine’’for the above three-sphere system because it can undergo a directional motion in the presence of hydrodynamic
interactions [11]. The effects of hydrodynamic interactions will be separately discussed in Section 8.
To describe the configuration of a micromachine, it is convenient to introduce the following two spring extensions
with respect to ℓ:
r12 = x2 − x1 − ℓ,r23 = x3 − x2 − ℓ.

(20)

From Eqs. (15)–(17), we obtain the reduced Langevin equations for r12 (t) and r23 (t) as [27]
dr12
dt
dr23
dt

=−
=

K

ζ12

K

ζ2

r12 +

r12 −

K

ζ2

K

ζ23

(
r23 +

(
r23 +

2T12

)1/2

ζ12
)1/2

2T23

ζ23

ξ12 ,

(21)

ξ23 .

(22)

Here we have introduced the relevant effective friction coefficient

ζij =

ζi ζj
,
ζi + ζj(23)

and the friction-weighted average temperature
Tij =ζj T i + ζi T j.
ζi + ζj(24)

The definition of the effective temperature Tij arises from the requirement that the newly introduced noises ξ12 (t) andξ23 (t) in Eqs. (21) and (22), respectively, satisfy the following statistical properties:⟨ξ12 (t)⟩ = ⟨ξ23 (t)⟩ = 0,(25)

⟨ξ12 (t)ξ12 (t )⟩ = δ (t − t ),(26)

⟨ξ23 (t)ξ23 (t )⟩ = δ (t − t ),()1/2
T2 ζ12 ζ23⟨ξ12 (t)ξ23 (t ′ )⟩ = −δ (t − t ′ ).
ζ2 T12 T23

(27)

′
′

′
′

(28)

The reduced Langevin equations in Eqs. (21) and (22) can be conveniently represented in the matrix form of Eq. (2).
Using the notations for the two-dimensional vectors r = (r12 , r23 )T and ξ = (ξ12 , ξ23 )T , the 2 × 2 matrices A, F, and B are
4

I. Sou, Y. Hosaka, K. Yasuda et al.

Physica A 562 (2021) 125277

now given by(A=−K /ζ12
K /ζ2

⎛(

2T12

K /ζ2
−K /ζ23

)

(29)

)1/2

⎞
0

ζ12

⎜

F=⎜
⎝

,

(
0

2T23

⎟
)1/2 ⎟
⎠,

(30)

ζ23

and

⎛
1

⎜

B=⎜
⎝

−

T2

(

ζ2

T2

−

ζ12 ζ23

ζ2

)1/2

(

ζ12 ζ23)1/2 ⎞

T12 T23⎟⎟,⎠

1

T12 T23(31)

respectively. Then, according to Eq. (5), the 2 × 2 diffusion matrix D becomes(D=

T12 /ζ12−T2 /ζ2)−T2 /ζ2.
T23 /ζ23(32)

In our previous work, we obtained the steady-state probability distribution function p0 (r) for a three-sphere micromachine [21]. Owing to the reproductive property of Gaussian distributions [15,16], p0 (r) should also be a Gaussian function
for the present linear problem and is given by

[

]

1
1
p0 (r) = N0 exp − r T C−
0 r .
2

(33)

In the above, N0 is the normalization factor, and the steady-state covariance matrix C0 is given by
C0 =

1

(

T12 + ζ12 ∆

ζ2 ∆

K

)
ζ2 ∆
,
T23 + ζ23 ∆

(34)

with

∆=

ζ12 ζ23 (T12 + T23 − 2T2 )(ζ12 + ζ23 )(ζ22 − ζ12 ζ23 )

.

(35)

In the following sections, the above matrices are used to calculate the average irreversibility and the average entropy
production rate.
4. Irreversibility
In this section, we calculate the average irreversibility ⟨σ (t)⟩ in Eq. (10) for a thermally driven three-sphere micromachine by using A and C0 in Eqs. (29) and (34), respectively (notice that Ct is also given by A and C0 according to Eq. (12)).
In Fig. 2, we plot ⟨σ (t)⟩ as a function of dimensionless time Kt /ζ2 . We also define the dimensionless temperature of thethree spheres by

τi =

2Ti
K ℓ2

,

(36)

which is the ratio between the thermal energy of each sphere and the spring elastic energy (recall kB = 1).
The chosen parameters in Fig. 2(a) are ζ1 = ζ2 = ζ3 , τ1 = 1/900, τ2 = 41/900, and τ3 = 81/900, and those in Fig. 2(b)are ζ1 = ζ2 = ζ3 , τ1 = τ3 = 25/900, and τ2 = 41/900. For these two cases, the friction coefficients are all identical.
In Fig. 2(a) for which τ1 ̸ = τ3 , the average irreversibility first increases from zero and it vanishes in the long time limitt → ∞. This is because p(x0 , x0 , 0) = p(x1 , x1 , 0) for t = 0, whereas the transition probability becomes a stationary one
for t → ∞. Since the average irreversibility is positive semidefinite for all time and goes to zero as t → 0 and t → ∞,there is a typical time Kt ∗ /ζ2 ≈ 0.39 for which ⟨σ (t)⟩ takes a global maximum value. Notice that ζ2 /K corresponds to thespring relaxation time, and t ∗ gives a characteristic time scale of the irreversible fluctuations. When τ1 = τ3 as in Fig. 2(b),on the other hand, ⟨σ (t)⟩ vanishes for all t. In this case, the system is in apparent equilibrium because the micromachineis thermally balanced [21].
Keeping the temperature parameters τi the same as in Fig. 2(a) and (b), we introduce asymmetry in the friction
coefficients in Fig. 2(c) and (d) for which we set ζ1 /ζ2 = 0.5 and ζ3 /ζ2 = 5. Although the time evolutions in Fig. 2(a)and (c) are similar, the maximum value of ⟨σ ⟩ in Fig. 2(c) is about 4.5 times larger than that in Fig. 2(a). Hence theasymmetry in the friction coefficients increases the average irreversibility. In Fig. 2(d), the irreversibility is non-zero and5

I. Sou, Y. Hosaka, K. Yasuda et al.

Physica A 562 (2021) 125277

Fig. 2. The average irreversibility ⟨σ (t)⟩ given by Eq. (10) as a function of dimensionless time Kt /ζ2 . The parameters are (a) ζ1 = ζ2 = ζ3 , τ1 = 1/900,
τ2 = 41/900, and τ3 = 81/900; (b) ζ1 = ζ2 = ζ3 , τ1 = τ3 = 25/900, and τ2 = 41/900; (c) ζ1 /ζ2 = 0.5, ζ3 /ζ2 = 5, τ1 = 1/900, τ2 = 41/900, andτ3 = 81/900 (these temperatures are the same as in (a)); (d) ζ1 /ζ2 = 0.5, ζ3 /ζ2 = 5, τ1 = τ3 = 25/900, and τ2 = 41/900 (these temperatures are
the same as in (b)).

Fig. 3. The dimensionless average entropy production rate ζ ⟨σ̇ ⟩/K given by Eq. (38) as a function of T1 /T2 and T3 /T2 . The average entropy production
rate vanishes along the diagonal white line T1 = T3 .

⟨σ (t)⟩ ≥ 0 even τ1 = τ3 . This is because the system is not thermally balanced, namely, T12 ̸= T23 , and the micromachineis in out-of-equilibrium. It is worth mentioning that the average irreversibility ⟨σ (t)⟩ depends only on the ratios of thetemperatures such as T1 /T2 and T3 /T2 .
In general, the irreversibility σ is a measure of non-equilibrium fluctuations, and one can distinguish forward fromreverse trajectories segments when σ is non-zero. The result in Fig. 2 clearly indicates that the asymmetry of either
the temperatures or the friction coefficients are necessary to drive a three-sphere micromachine out of equilibrium. The
maximum value of ⟨σ ⟩ in Fig. 2 directly reflects the magnitude of the irreversibility. It is useful to note that the asymmetryof both the temperatures and the friction coefficients can enhance the maximum value of ⟨σ ⟩ (compare Figs. 2(a) and (c)).
As we will discuss in the next section, the larger the maximum value of ⟨σ ⟩ is, the larger the average entropy production
rate of a thermally driven micromachine becomes.
6

I. Sou, Y. Hosaka, K. Yasuda et al.

Physica A 562 (2021) 125277

Fig. 4. The dimensionless entropy production rate ζ ⟨σ̇ ⟩/K given by Eq. (38) as a function of T2 /T1 for T3 /T1 = 1, 3, 5, and 7.

5. Entropy production rate
Next we calculate the steady-state average entropy production rate of a thermally driven micromachine. In thesteady-state, the entropy production rate balances with the entropy extraction rate, and both quantities become zeroat equilibrium [32]. Substituting A, D, and C0 in Eqs. (29), (32), and (34), respectively, to Eqs. (13) and (14), we obtainafter some calculation⟨σ̇ ⟩ =

K [ζ1 (T3 − T2 ) + ζ3 (T2 − T1 )]2(ζ1 T2 T3 + ζ2 T3 T1 + ζ3 T1 T2 )(ζ1 ζ2 + ζ2 ζ3 + 2ζ1 ζ3 )

.

(37)

This is an important result of this paper. Obviously, we have ⟨σ̇ ⟩ ≥ 0. When the system is in thermal equilibrium,i.e., T1 = T2 = T3 , the entropy production rate vanishes for any combination of the friction coefficients. We note thatEq. (37) does not depend on the spring natural length ℓ because we are considering only the fluctuations around ℓ as
defined in Eq. (20).
When the three friction coefficients are all identical, i.e., ζ1 = ζ2 = ζ3 = ζ , Eq. (37) reduces to⟨σ̇ ⟩ =

K (T1 − T3 )2
4ζ (T1 T2 + T2 T3 + T3 T1 )

.

(38)

Using Eq. (38), we give in Fig. 3 a color representation of the dimensionless average entropy production rate ζ ⟨σ̇ ⟩/K as afunction of T1 /T2 and T3 /T2 . Here ⟨σ̇ ⟩ vanishes when T1 = T3 and it increases as the difference between T1 and T3 becomeslarger. To see the role of the temperature T2 of the middle sphere in Eq. (38), we plot in Fig. 4 the steady-state average
entropy production rate ζ ⟨σ̇ ⟩/K as a function of T2 /T1 (not T1 /T2 ) for four different values of T3 /T1 . From this plot, one
can clearly see that ⟨σ̇ ⟩ becomes smaller as T2 /T1 is increased.
When the friction coefficients are different between the three spheres, Eq. (37) implies that the average entropy
production rate is non-zero, ⟨σ̇ ⟩ > 0, even when T1 = T3 . In general, ⟨σ̇ ⟩ becomes larger when either ζ1 /ζ2 or ζ3 /ζ2 isincreased. As an example of asymmetric situations, we give in Fig. 5 a color representation of the dimensionless average
entropy production rate ζ2 ⟨σ̇ ⟩/K as a function of ζ1 /ζ2 and ζ3 /ζ2 when T1 /T2 = 0.2 and T3 /T2 = 5. The value of ⟨σ̇ ⟩ isasymmetric with respect to the line ζ1 = ζ3 , and it becomes larger when ζ3 /ζ2 becomes smaller.
In multidimensional systems, it was generally discussed that enhanced fluctuation (or noise amplification) occurs not
through some additional forcing, but through violation of detailed balance which is measured by the gain matrix G in
Eq. (14) [22]. Hence, according to Eq. (13), the entropy production rate reflects the noise amplification that occurs when
detailed balance is not satisfied [23]. It is also known that the entropy production rate is related to the heat flow in asystem [30].
As we can see from Eq. (1) for the average velocity ⟨V ⟩ and Eq. (38) for the average entropy production rate ⟨σ̇ ⟩ of athermally driven three-sphere micromachine, the temperature difference between the first and the third spheres, T1 − T3 ,plays an essential role to characterize its non-equilibrium behaviors. According to their dependence on T1 − T3 , we see
a proportionality relation such that ⟨σ̇ ⟩ ∼ ⟨V ⟩2 . In our previous work [11], we showed that the average velocity ⟨V ⟩ of athree-sphere micromachine is proportional to the net heat flow between the first and the third spheres.
7

I. Sou, Y. Hosaka, K. Yasuda et al.

Physica A 562 (2021) 125277

Fig. 5. The dimensionless entropy production rate ζ2 ⟨σ̇ ⟩/K given by Eq. (37) as a function of ζ1 /ζ2 and ζ3 /ζ2 . Here the temperature ratios are fixed
to T1 /T2 = 0.2 and T3 /T2 = 5.

6. Diffusion coefficient
In this section, we discuss the Brownian motion of a thermally driven three-sphere micromachine. We introduce thecenter of mass position of a micromachine by1

X (t) =

3[x1 (t) + x2 (t) + x3 (t)],(39)

even though we neglect the inertia of the spheres. From Eqs. (15)–(17), the Langevin equation for X can be written interms of r12 and r23 as
dX
dt

=

K

(

3

ζ2 − ζ1
ζ1 ζ2

)
r12 +

K

(

3

ζ3 − ζ2
ζ2 ζ3

√ (

)
r23 +

2

T1

3

ζ1

+

T2

ζ2

+

T3

)1/2

ζ3

ξX ,

(40)

where ξX (t) is a zero mean and unit variance Gaussian white noise defined by

ξX =

(

T1

+

ζ1

T2

ζ2

+

T3

)−1/2 [(

ζ3

T1

ζ1

)1/2

ξ1 +

(

T2

ζ2

)1/2

ξ2 +

(

T3

ζ3

)1/2

]
ξ3 ,

(41)

and satisfies the following statistical properties⟨ξX (t)⟩ = 0,

(42)

⟨ξX (t)ξX (t )⟩ = δ (t − t ).
′

′

(43)

Using Eqs. (21) and (22) for the dynamics of r12 and r23 , respectively, the mean squared displacement of the center ofmass position becomes⟨X 2 (t)⟩ = 2Dt ,(44)

where the diffusion coefficient is obtained asD=ζ1 T1 + ζ2 T2 + ζ3 T3.
(ζ1 + ζ2 + ζ3 )2(45)

See Appendix B for the detailed derivation. When the temperatures are all identical, T1 = T2 = T3 , Eq. (45) becomesD=T2

ζ1 + ζ2 + ζ3

,

(46)

as expected for the equilibrium case.
On the other hand, when the three friction coefficients are all identical, ζ1 = ζ2 = ζ3 = ζ , Eq. (45) reduces toD=

T1 + T2 + T3
9ζ

,

(47)
8

I. Sou, Y. Hosaka, K. Yasuda et al.

Physica A 562 (2021) 125277

which is proportional to the sum of the three temperatures. Introducing an average temperature by
TX =T1 + T2 + T3
3

,

(48)

we can rewrite Eq. (47) as D = TX /(3ζ ).
In general, the diffusion coefficient of a thermally driven n-sphere swimmer is predicted to be∑n

Dn = (∑i=1

n
i=1

ζi Ti
)2 .
ζi

(49)

This expression can be explicitly confirmed also for n = 2. For n = 2, the diffusion coefficient for the center of friction
was obtained in Ref. [27], and it is different from that of the center of mass diffusion. This difference is not physicallyessential because it only depends on the choice of the coordinate system. Nevertheless, we consider that the center ofmass diffusion is much easier to be measured in the experiments.
7. Efficiency
Finally, we shall estimate the efficiency of a thermally driven micromachine. When the three friction coefficients are
all identical and given ζ , the average velocity is given by Eq. (1) which can be rewritten as [11]⟨V ⟩ =kB (T3 − T1 )a16ζ ℓ2.

(50)

Here a is the radius of the spheres, and we employ the Stokes relation ζ = 6πηa for the friction coefficient. Notice that
we recover the Boltzmann constant kB in this section for the sake of clarity. The above result indicates that the swimming
direction is from a colder sphere to a hotter one, and the velocity does not depend on the temperature of the middle
sphere [11].
Following Ref. [5], we define the efficiency of a thermally driven micromachine byε=3ζ ⟨V ⟩2,kB ⟨Σ̇ ⟩TX(51)

where ⟨Σ̇ ⟩ is the total entropy production rate and the average temperature TX is given by Eq. (48). Notice that theentropy production rate ⟨σ̇ ⟩ in Eq. (38) takes into account only the internal motions (r12 and r23 ) of a micromachine andit also vanishes when T1 = T3 .
In order to obtain ⟨Σ̇ ⟩, one needs to take into account the center of mass motion X given by Eq. (39), and solve
Eqs. (21), (22), and (40) simultaneously. For the three-dimensional vectors r = (r12 , r23 , X )T and ξ = (ξ12 , ξ23 , ξX )T , thecorresponding 3 × 3 matrices A, F, and D are shown in Appendix C. Repeating the same calculation as in Section 5, we
obtain the following total entropy production rate⟨Σ̇ ⟩ =

K
12ζ T1 T2 T3(T12 T2 + 4T1 T22 + 4T12 T3 + T32 T2 + 4T3 T22 + 4T1 T32 − 18T1 T2 T3 ),(52)

when the friction coefficients are identical. Unlike ⟨σ̇ ⟩ in Eq. (38), ⟨Σ̇ ⟩ in Eq. (52) vanishes only when T1 = T2 = T3 ,i.e., thermal equilibrium. It should be noted here that, although the above total entropy production rate ⟨Σ̇ ⟩ takes into
account all the three positional degrees of freedom of the spheres, it still does not include hydrodynamic interactions
acting between different spheres, which have been neglected throughout this paper. Hence, the maximum of ε in Eq. (51)
is not necessarily unity. Nevertheless, ε gives a useful measure of the ratio between the hydrodynamic dissipation of a
micromachine and its total entropy production rate due to the sphere motions.
With the above result, the efficiency ε in Eq. (51) can be obtained asε=27a2 kB T1 T2 T3 (T1 − T3 )2
64K ℓ

4 (T

1

+ T2 +T3 )(T12 T2+ 4T1 T22 + 4T12 T3 + T32 T2 + 4T3 T22 + 4T1 T32 − 18T1 T2 T3 )

.

(53)

Clearly, ε vanishes when T1 = T3 as it should. The important outcome of Eq. (53) is that the efficiency scales as (a/ℓ)2
and is proportional to the temperature T2 of the middle sphere. By using the dimensionless temperature τ2 = 2T2 /(K ℓ2 ),we give in Fig. 6 a color representation of the scaled efficiency 128ℓ2 ε/(27a2 τ2 ) as a function of T1 /T2 and T3 /T2 . The
efficiency also becomes smaller when the temperatures are highly asymmetric.
On the other hand, the efficiency ε becomes larger along a certain characteristic curve. Such a dependence on thetemperatures can be explained as follows. The numerator of ε (hydrodynamic dissipation) always vanishes for T1 = T3 ,and all the contour lines are parallel to the line T1 = T3 . On the other hand, the denominator of ε (total entropy productionrate) increases when T1 and T3 are either small or large even for T1 = T3 (the minimum occurs at T1 = T2 = T3 ), andthe contour lines are perpendicular to the line T1 = T3 . These different temperature dependences between the numerator
and the denominator give rise to the nontrivial increase of ε as shown in Fig. 6.
9I. Sou, Y. Hosaka, K. Yasuda et al.

Physica A 562 (2021) 125277

Fig. 6. The scaled efficiency 128ℓ2 ε/(27a2 τ2 ) given by Eq. (53) as a function of T1 /T2 and T3 /T2 . The efficiency vanishes along the diagonal whiteline T1 = T3 .

8. Summary and discussion
In this paper, we have discussed the non-equilibrium behaviors of a thermally driven elastic three-sphere micromachine. In our model, the three spheres are in contact with independent heat baths having different temperatures [11,21].
Using the formulation of a linear stochastic Langevin dynamics [22,23], we have calculated the time-dependent average
irreversibility ⟨σ (t)⟩ in Eq. (10). When the temperatures and/or the friction coefficients are asymmetric, the average
irreversibility is non-zero and takes a global maximum value for a finite time. The corresponding characteristic timescale is roughly set by the spring relaxation time ζ2 /K . We have further obtained the average entropy production rate ⟨σ̇ ⟩in Eq. (37) which is the zero-time growth rate of the average irreversibility. This quantity decreases as the temperatureof the middle phase increases.
We have also discussed the Brownian motion of a thermally driven three-sphere micromachine and calculated its
center of mass diffusion coefficient D as in Eq. (45). The obtained expression can be generalized for a many-spheremicromachine. When the friction coefficients are identical, an average temperature TX can be introduced as in Eq. (48). Our
result is different from the diffusion coefficient for the center of friction obtained for a non-equilibrium dimer model [27].
Finally, with the results of the total entropy production rate and the average temperature, we have estimated the efficiency
of a micromachine in Eq. (53).
Our model of a three-sphere micromachine has a similarity to that of two over-damped, tethered spheres coupled bya harmonic spring and also confined between two walls [17,20]. In these works, the authors numerically showed that
displacements obey a Gaussian distribution and also found probability flux loops that demonstrate the broken detailedbalance [17,20]. The two displacements r12 and r23 in Eq. (20) correspond to the sphere positions in their model. However,the presence of the middle sphere changes the structure of the frequency matrix for a three-sphere micromachine
when T2 ̸ = 0 [21]. Moreover, a two-sphere micromachine in a viscous fluid cannot have a directed motion even if thetemperatures are different [11]. Recently, Li et al. used the two-sphere model to calculate the entropy production rate [30].
We note that our result in Eq. (38) reduces to their expression when T2 = 0.
In this work, we have neglected long-ranged hydrodynamic interactions acting between different spheres and wehave not considered the locomotion of a micromachine [11]. If hydrodynamic interactions are taken into account in
the present analysis, the covariance matrix in Eq. (34) is modified in non-equilibrium situations. Such hydrodynamic
corrections should be proportional to a/ℓ within the lowest-order expansion. Moreover, these corrections should vanishin thermal equilibrium, i.e., T1 = T2 = T3 because hydrodynamic interactions should not affect equilibrium statisticalproperties.
CRediT authorship contribution statement
Isamu Sou: Investigation, Formal analysis, Writing - original draft. Yuto Hosaka: Conceptualization, Methodology. Kento Yasuda: Conceptualization, Methodology. Shigeyuki Komura: Conceptualization, Methodology, Supervision,
Writing - review & editing.
10

I. Sou, Y. Hosaka, K. Yasuda et al.

Physica A 562 (2021) 125277

Declaration of competing interest
The authors declare that they have no known competing financial interests or personal relationships that could haveappeared to influence the work reported in this paper.
Acknowledgments
Y.H. acknowledges support by a Grant-in-Aid for JSPS Fellows (Grant No. 19J20271) from the Japan Society for the
Promotion of Science (JSPS). K.Y. acknowledges support by a Grant-in-Aid for JSPS Fellows (Grant No. 18J21231) fromthe JSPS. S.K. acknowledges support by a Grant-in-Aid for Scientific Research (C) (Grant No. 18K03567 and Grant No.
19K03765) from the JSPS, and support by a Grant-in-Aid for Scientific Research on Innovative Areas ‘‘Information physics
of living matters’’ (Grant No. 20H05538) from the Ministry of Education, Culture, Sports, Science and Technology of Japan.
Appendix A. Irreversibility distribution P(σ )Following Ref. [23], we show how to calculate the probability density function of irreversibility P(σ ) for a systemdescribed by the linear Langevin dynamics in Eq. (2). We remind that P(σ ) satisfies the fluctuation theorem in Eq. (9).
We first introduce a 2N-dimensional state space

( )
r0
r1

z=

,

(A.1)

where r 0 and r 1 are the initial and final states. The probability density function of irreversibility P(σ ) is given by∫P(σ ) =

d2N z δ (σ − z T Rz /2)p(r 0 , r 1 , t),(A.2)

where R = R10 − R01 is the 2N × 2N matrix and we have used the matrices(
R 10 =1
C−t
AT t −1
−e Ct

)
1 At
−C−
t e
T
1 At
1 ,
eA t C−
+ C−
t e
0

(A.3)

and

(
R 01 =

T

T

1 At
1
eA t C−
+ C−
t e
0
−1 At
−Ct e

1
−eA t C−
t

)

1
C−
t

.

(A.4)

The characteristic function of the probability density function is defined by∫dσ P(σ )eikσ .

P [k] =(A.5)

For linear Langevin systems, it is shown that the characteristic function can be expressed as [23]1
P [k] = ∏2N √,1 − ikλm
m=1

(A.6)

1
where λm are the 2N eigenvalues of the matrix R −01 R. Then the probability density function of irreversibility in Eq. (A.2)
can be obtained by the inverse transform of Eq. (A.7):

1

P(σ ) =∫

2π

dk P [k]e−ikσ .

(A.7)

We discuss here the four eigenvalues for a thermally driven micromachine. When T1 = T3 , we find that all the
eigenvalues vanish, i.e., λm = 0. In this thermally balanced situation, the characteristic function is simply P [k] = 1 andthe probability density function of irreversibility is P(σ ) = δ (σ ).
As the simplest non-equilibrium situation, we consider the case when T1 = T2 = 0 but T3 ̸ = 0. Then the four
eigenvalues can be obtained in the short time limit as

√

λ1 ≈

3

√

(2t̄ − 1)(t̄ + 2) t̄(23t̄ + 8),

(A.8)

8√

λ2 ≈ −λ3 ≈ −
λ4 ≈ −

3√(2t̄ − 1)(t̄ + 2) t̄(23t̄ + 8),8[(2t̄ − 1)t̄
(2t̄ − 1)
t̄

√
3−

3
8

(t̄ + 2)2√√](t̄ + 2)(7t̄ + 6) ,

√

[
3+

3
8

(t̄ + 2)2

(A.9)

]

(A.10)(t̄ + 2)(7t̄ + 6) ,(A.11)

where t̄ = Kt /ζ2 is the dimensionless time. The above expressions are valid when t̄ ≪ 1.
11

I. Sou, Y. Hosaka, K. Yasuda et al.

Physica A 562 (2021) 125277

Appendix B. Derivation of diffusion coefficient
In this Appendix, we show the derivation of the diffusion coefficient in Eq. (45). We first integrate Eq. (40) over time,and obtain the mean squared displacement of a three-sphere micromachine as
K2

2⟨X (t)⟩ =

(

9
K

(
2

ζ2 − ζ1
ζ1 ζ2)2 ∫ t00

ζ3 − ζ2
ζ2 ζ3)2 ∫ t

∫ t

∫ t
dt1

dt2 ⟨r12 (t1 )r12 (t2 )⟩dt1
dt2 ⟨r23 (t1 )r23 (t2 )⟩)∫ t 0 ∫ t
T2
T3
2 T1
dt1
dt2 ⟨ξX (t1 )ξX (t2 )⟩
+
+
+
9 ζ1
ζ2
ζ3
0
0
(
)(
)∫ t
∫ t
2K 2 ζ2 − ζ1
ζ3 − ζ2+dt1
dt2 ⟨r12 (t1 )r23 (t2 )⟩
9
ζ1 ζ2
ζ2 ζ3
0
0
√ (
)(
)1/2 ∫ t∫ t
T2
T3
2 2K ζ2 − ζ1
T1
+
+
+
dt1
dt2 ⟨r12 (t1 )ξX (t2 )⟩
9
ζ1 ζ2
ζ1
ζ2
ζ3
0
0
√ (
)(
)1/2 ∫ t
∫ t
T1
T2
T3
2 2K ζ3 − ζ2dt2 ⟨r23 (t1 )ξX (t2 )⟩.
dt1
+
+
+
9
ζ2 ζ3
ζ1
ζ2
ζ3
0
0

+

9

0

(

(B.1)

Our task is to calculate the various noise correlation functions in the above expression.
Let us introduce the Fourier transform of a function f (t) by∫ ∞f [ω] =

dt f (t)eiωt ,f (t) =

−∞

1

∫ ∞

2π

−∞

dω f [ω]e−iωt .

(B.2)

Then we can solve the Langevin equations in Eqs. (21) and (22) in the Fourier domain as

)(
)1/2
(
)1/2
ζ2
iζ2 ω
2T12
2T23
+
ξ12 [ω] +ξ23 [ω]
ζ23
K
ζ12
ζ23
r12 [ω] = −
(
)
(
) ,
ζ22
1
K
1
ζ2 2
ω+
1−
ω − iζ2
+
K
ζ12
ζ23
ζ2
ζ12 ζ23
(
)1/2
(
)1/2
)(
ζ2
2T12
iζ2 ω
2T23
ξ23 [ω] +ξ12 [ω]
+
ζ12
K
ζ23
ζ12
r23 [ω] = −
(
)
(
) .
ζ22
1
1
K
ζ2 2
ω − iζ2
+
ω+
1−
K
ζ12
ζ23
ζ2
ζ12 ζ23
(

(B.3)

(B.4)

Calculating the products of the noise and taking the average, we obtain for example[⟨r12 (t1 )r12 (t2 )⟩ =(ζ2 ζ12 + ζ2 ζ23 + ζ12 ζ23 H)T12

2ζ12 (ζ22 T12 + ζ12 ζ23 T23 − 2ζ12 ζ23 T2 )]−2KH ζ12 (ζ12 + ζ23 )KH(ζ12 + ζ23 )(ζ2 ζ12 + ζ2 ζ23 + ζ12 ζ23 H)[]
K (ζ2 ζ12 + ζ2 ζ23 + ζ12 ζ23 H)
× exp −
|t1 − t2 |2ζ2 ζ12 ζ23[]
2ζ12 (ζ22 T12 + ζ12 ζ23 T23 − 2ζ12 ζ23 T2 )(ζ2 ζ12 + ζ2 ζ23 − ζ12 ζ23 H)T12+ −+2KH ζ12 (ζ12 + ζ23 )KH(ζ12 + ζ23 )(ζ2 ζ12 + ζ2 ζ23 − ζ12 ζ23 H)[]
K (ζ2 ζ12 + ζ2 ζ23 − ζ12 ζ23 H)
× exp −
|t1 − t2 | ,2ζ2 ζ12 ζ23

(B.5)

where

(
H=

4+

ζ22
ζ2
2ζ22
+ 22 −
2ζ12 ζ23
ζ12
ζ23

)1/2

.

(B.6)

The other noise correlation functions can be obtained in a similar way. However, it should be noted that, for nonequilibrium situations, the time-reversal invariance is not generally satisfied for the cross correlation functions, i.e.,⟨r12 (t1 )r23 (t2 )⟩ ̸= ⟨r12 (t2 )r23 (t1 )⟩. On the other hand, the time-translational invariance of the noise correlation functions isalways satisfied because we are dealing with steady-states.
Collecting all the noise correlation functions and taking the limit of t → ∞, we finally obtain
D = limt →∞⟨X 2 (t)⟩2t=ζ12 (ζ1 − ζ2 )2 (ζ22 T12 + ζ12 ζ23 T23 − 2ζ12 ζ23 T2 )9ζ12 (ζ22 − ζ12 ζ23 )2
12

I. Sou, Y. Hosaka, K. Yasuda et al.

Physica A 562 (2021) 125277

ζ23 (ζ3 − ζ2 )2 (ζ22 T23 + ζ12 ζ23 T12 − 2ζ12 ζ23 T2 )9ζ32 (ζ22 − ζ12 ζ23 )2
(
)
T2
T3
1 T1
+
+
+
9 ζ1
ζ2
ζ3
2ζ12 ζ23 (ζ1 − ζ2 )(ζ2 − ζ3 )(ζ22 T12 + ζ22 T23 − ζ22 T2 − ζ12 ζ23 T2 )+9ζ1 ζ2 ζ3 (ζ22 − ζ12 ζ23 )2
2ζ12 (ζ2 − ζ1 )(ζ1 ζ2 ζ23 T3 − ζ1 ζ3 ζ23 T2 + ζ1 ζ2 ζ3 T2 − ζ22 ζ3 T1 )+9ζ12 ζ2 ζ3 (ζ22 − ζ12 ζ23 )2ζ23 (ζ2 − ζ3 )(ζ2 ζ3 ζ12 T1 − ζ1 ζ3 ζ12 T2 + ζ1 ζ2 ζ3 T2 − ζ1 ζ22 T3 ).
+9ζ1 ζ2 ζ32 (ζ22 − ζ12 ζ23 )
+

(B.7)

Notice that the different lines in Eq. (B.1) correspond to the different lines in the above equation. The right hand side ofEq. (B.7) reduces to Eq. (45).
Appendix C. 3 × 3 matrices
Let us consider the three coupled Langevin equations Eqs. (21), (22), and (40). By introducing the three-dimensional
vectors r = (r12 , r23 , X )T and ξ = (ξ12 , ξ23 , ξX )T , the corresponding 3 × 3 matrices A, F, and D are given by−K /ζ12⎜ ( K /ζ2 )A = ⎝K ζ − ζ
2
1
3
ζ1 ζ2

K /ζ2
−
( K /ζ23 )

⎛

⎛(
⎜
⎜
⎜
⎜
F=⎜
⎜
⎜
⎝

2T12

ζ3 − ζ2
ζ2 ζ3

K
3

0
0⎟

⎞
⎠,

(C.1)

0

)1/2

⎞
0

ζ12

(
0
0

2T23

0

⎟
⎟
⎟
⎟
⎟,
0
⎟
√ (
)1/2 ⎟
⎠
2 T1
T2
T3
+
+
3
ζ1
ζ2
ζ3

)1/2

ζ23
0

(C.2)

and

⎛

⎜ T12 /ζ12
⎜
⎜
⎜
D=⎜
⎜ −T2 /ζ2
⎜ (
)
⎝1 T
T1
2
−
3 ζ2
ζ1

1

−T2 /ζ2

3

(

T3

ζ3

−

T2

ζ2

−

T1

)

⎞

ζ1

⎟
⎟
⎟
⎟
1 T3
T2
⎟,
−
⎟
3 ζ3
ζ2
(
)⎟
1 T1
T2
T3 ⎠
+
+
9 ζ1
ζ2
ζ3
(

)

T2

ζ2

3

T23 /ζ231()

respectively. The above matrices are the generalization of Eqs. (29), (30), and (32) to a higher dimension.
References
[1]
[2]
[3]
[4]
[5]
[6]
[7]
[8]
[9]
[10]
[11]
[12]
[13]
[14]
[15]
[16]
[17]
[18]

E. Lauga, T.R. Powers, Rep. Progr. Phys. 72 (2009) 096601.
E.M. Purcell, Proc. Natl. Acad. Sci. USA 94 (1997) 11307.
E. Lauga, Soft Matter 7 (2011) 3060.
A. Najafi, R. Golestanian, Phys. Rev. E 69 (2004) 062901.
R. Golestanian, A. Ajdari, Phys. Rev. E 77 (2008) 036308.
M. Leoni, J. Kotar, B. Bassetti, P. Cicuta, M.C. Lagomarsino, Soft Matter 5 (2009) 472.
G. Grosjean, M. Hubert, G. Lagubeau, N. Vandewalle, Phys. Rev. E 94 (2016) 021101(R).
G. Grosjean, M. Hubert, N. Vandewalle, Adv. Colloid Interface Sci. 255 (2018) 84.
K. Yasuda, Y. Hosaka, M. Kuroda, R. Okamoto, S. Komura, J. Phys. Soc. Japan 86 (2017) 093801.
M. Kuroda, K. Yasuda, S. Komura, J. Phys. Soc. Japan 88 (2019) 054804.
Y. Hosaka, K. Yasuda, I. Sou, R. Okamoto, S. Komura, J. Phys. Soc. Japan 86 (2017) 113801.
K. Sekimoto, J. Phys. Soc. Japan 66 (1997) 1234.
K. Sekimoto, Progr. Theoret. Phys. Suppl. 130 (1998) 17.
K. Sekimoto, Stochastic Energetics, Springer, Berlin Heidelberg, 2010.
N. van Kampen, Stochastic Processes in Physics and Chemistry, North-Holland, Amsterdam, 1992.
H. Risken, The Fokker–Planck Equation: Methods of Solution and Applications, Springer-Verlag, Berlin, 1984.
C. Battle, C.P. Broedersz, N. Fakhri, V.F. Geyer, J. Howard, C.F. Schmidt, F.C. MacKintosh, Science 352 (2016) 604.
J. Gladrow, N. Fakhri, F.C. MacKintosh, C.F. Schmidt, C.P. Broedersz, Phys. Rev. Lett. 116 (2016) 248301.
13

(C.3)

I. Sou, Y. Hosaka, K. Yasuda et al.
[19]
[20]
[21]
[22]
[23]
[24]
[25]
[26]
[27]
[28]
[29]
[30]
[31]
[32]

Physica A 562 (2021) 125277

J. Gladrow, C.P. Broedersz, C.F. Schmidt, Phys. Rev. E 96 (2017) 022408.
F.S. Gnesotto, F. Mura, J. Gladrow, C.P. Broedersz, Rep. Progr. Phys. 81 (2018) 066601.
I. Sou, Y. Hosaka, K. Yasuda, S. Komura, Phys. Rev. E 100 (2019) 022607.
J.B. Weiss, Tellus A 55 (2003) 208.
J.B. Weiss, Phys. Rev. E 76 (2007) 061128.
U. Seifert, Phys. Rev. Lett. 95 (2005) 040602.
C. Jarzynski, Annu. Rev. Condens. Matter Phys. 2 (2011) 329.
U. Seifert, Rep. Progr. Phys. 75 (2012) 126001.
A.Y. Grosberg, J.-F. Joanny, Phys. Rev. E 92 (2015) 032118.
R. Zwanzig, Nonequilibrium Statistical Mechanics, Oxford University Press, New York, 2001.
R.R. Netz, J. Chem. Phys. 148 (2018) 185101.
J. Li, J.M. Horowitz, T.R. Gingrich, N. Fakhri, Nat. Commun. 10 (2019) 1666.
V.Y. Chernyak, M. Chertkov, C. Jarzynski, J. Stat. Mech. Theory Exp. (2006) P08001.
M.A. Taye, Phys. Rev. E 92 (2015) 032126.

14