PHYSICAL REVIEW E 100, 022607 (2019)

Nonequilibrium probability flux of a thermally driven micromachine
Isamu Sou, Yuto Hosaka, Kento Yasuda, and Shigeyuki Komura*

Department of Chemistry, Graduate School of Science, Tokyo Metropolitan University, Tokyo 192-0397, Japan(Received 15 May 2019; revised manuscript received 5 August 2019; published 28 August 2019)
We discuss the nonequilibrium statistical mechanics of a thermally driven micromachine consisting of threespheres and two harmonic springs [Y. Hosaka et al., J. Phys. Soc. Jpn. 86, 113801 (2017)]. We obtain thenonequilibrium steady state probability distribution function of such a micromachine and calculate its probability
flux in the corresponding configuration space. The resulting probability flux can be expressed in terms of afrequency matrix that is used to distinguish between a nonequilibrium steady state and a thermal equilibrium
state satisfying detailed balance. The frequency matrix is shown to be proportional to the temperature difference
between the spheres. We obtain a linear relation between the eigenvalue of the frequency matrix and the averagevelocity of a thermally driven micromachine that can undergo a directed motion in a viscous fluid. This relationis consistent with the scallop theorem for a deterministic three-sphere microswimmer.
DOI: 10.1103/PhysRevE.100.022607

I. INTRODUCTION

Microswimmers are tiny machines that swim in a fluid andthey are expected to be used in microfluidics and microsystems [1]. Over the length scale of microswimmers, the fluid
forces acting on them are dominated by the frictional viscous
forces. By transforming chemical energy into mechanical energy, however, microswimmers change their shape and moveefficiently in viscous environments. According to Purcell’sscallop theorem, time-reversible body motion cannot be used
for locomotion in a Newtonian fluid [2,3]. As one of thesimplest models exhibiting broken time-reversal symmetry,
Najafi and Golestanian proposed a three-sphere swimmer[4,5], in which three in-line spheres are linked by two armsof varying length. Recently, 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]. A similar model was also considered byother people [11–13]. Later, our model was further extended
to a thermally driven elastic microswimmer [14], suggestinga mechanism for locomotion that is purely induced by thermal fluctuations without any external forcing. As depicted
in Fig. 1, the key assumption is that the three spheres arein equilibrium with independent heat baths characterized bydifferent temperatures. In such a situation, heat transfer occurs
inside the micromachine from a hotter sphere to a colder one,
driving the whole system out of equilibrium. We have shownthat a combination of heat transfer and hydrodynamic interactions among the spheres can lead to directional locomotion
in a steady state [14]. Our model has a similarity to a classof thermal ratchet models [15–17], and the suggested new*

komura@tmu.ac.jp

2470-0045/2019/100(2)/022607(10)mechanism is relevant to nonequilibrium dynamics of proteinsand enzymes in biological systems [18,19].
For a thermally driven elastic micromachine, the averagevelocity was calculated to be [14]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 the viscosity of the surrounding fluid, and  is the natural length of the two springs. This result indicates that theswimming direction is from a colder sphere to a hotter one,and the velocity does not depend on the temperature of the
middle sphere. Moreover, we demonstrated that the averagevelocity is determined by the net heat flow between the first
and the third spheres (see Sec. VI later) [14]. This result
is consistent with the theoretical framework of “stochastic
energetics” [20–22]. However, a more detailed analysis basedon nonequilibrium statistical mechanics is required in orderto clarify the physical mechanism for the locomotion of athermally driven micromachine.
It is well known that systems in thermodynamic equilibrium obey detailed balance meaning that transition rates between any two microscopic states are pairwise balanced [23].
For nonequilibrium steady state situations, however, detailedbalance is broken and a probability flux loop should exist
in a configuration phase space [24–27]. Such a probability
flux has been experimentally measured in the periodic beatingof a flagellum from Chlamydomonas reinhardtii and in the
nonperiodic fluctuations of primary cilia of epithelial cells[28]. A similar analysis was performed for nonequilibrium
shape fluctuations of semiflexible filaments in a viscoelastic
environment [29,30]. Since the existence of a probability flux
loop is a direct verification of a nonequilibrium steady state, it
is a useful concept to characterize driven systems.
In this paper, we discuss the nonequilibrium statisticalmechanics of a thermally driven micromachine consisting of
three spheres and two harmonic springs [14]. We obtain the022607-1©2019 American Physical Society

SOU, HOSAKA, YASUDA, AND KOMURA

x1
T1
ζ1

K

x2
T2
ζ2

K

PHYSICAL REVIEW E 100, 022607 (2019)energy is given by

x3
T3

E=

ζ3

FIG. 1. Thermally driven elastic three-sphere micromachine.
Three spheres are connected by two identical harmonic springs
characterized by the elastic 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 the ith sphere. The three spheres are in equilibriumwith independent heat baths at temperatures Ti . In this paper, we
do not take into account hydrodynamic interactions acting between
different spheres.

steady state conformational distribution function of such a
micromachine and calculate its probability flux in the corresponding configuration space. The obtained probability flux
will be expressed in terms of a frequency matrix to discuss thenonequilibrium steady state of a micromachine. The main purpose of our work is to understand the physical mechanism that
underlies the locomotion of a thermally driven micromachine
within the nonequilibrium statistical mechanics. To this aim,
we shall obtain a relation connecting the eigenvalue of thefrequency matrix to the average velocity of a thermally driven
micromachine as shown in Eq. (1). With this relation, we show
explicitly that the concept of Purcell’s scallop theorem canbe generalized for thermally driven micromachines. Together
with the results in Ref. [14], the present study provides aunified description of the locomotion of a stochastic elastic
micromachine.
In the next section, we explain our model of a thermally
driven three-sphere micromachine by introducing the coupled
Langevin equations for the two spring extensions. In Sec. III,
we determine the steady state probability distribution functionof an elastic micromachine. By employing the Fokker-Planck
equation in Sec. IV, we obtain the steady state probability fluxin the configuration space. From the Gaussian probability flux,
we calculate the frequency matrix and the flux rotor. In Sec. V,
the average velocity of a micromachine will be expressed interms of the obtained quantities characterizing the scale ofnonequilibrium. Finally, a summary of our work and somediscussion are given in Sec. VI. In the Appendix, we give amatrix representation of linear stochastic dynamical systems,
which is useful to understand our results from a general pointof view.

KK(x2 − x1 − )2 + (x3 − x2 − )2 ,22

where xi (t ) (i = 1, 2, 3) are the positions of the three spheresin a one-dimensional coordinate system, and we also assume
x1 < x2 < x3 without loss of generality. In our previous model
for an elastic swimmer [9,10], the natural length of eachspring was assumed to undergo a prescribed cyclic motion intime, (t ), representing internal states of the micromachine.
However, for a thermally driven micromachine [14], as we
discuss in this paper,  is taken to be a constant.
We consider a situation in which the three spheres arein thermal equilibrium with independent heat baths at temperatures Ti (i = 1, 2, 3) [14]. When these temperatures aredifferent, the system is inevitably driven out of equilibrium
because heat flux from a hotter sphere to a colder one is
generated. The equations of motion of the three spheres are
written in the form of Langevin equations asK
2T1 1/2(x2 − x1 − ) +
ξ1 ,
(3)
ζ1
ζ1


2T2 1/2KK
ẋ2 = − (x2 − x1 − ) + (x3 − x2 − ) +
ξ2 ,
ζ2
ζ2
ζ2
(4)
1/2

K
2T3
ẋ3 = − (x3 − x2 − ) +
ξ3 ,
(5)
ζ3
ζ3
ẋ1 =

where the dot indicates the time derivative, ζi (i = 1, 2, 3) isthe friction coefficient for the ith sphere, and the Boltzmann
constant is set to unity hereafter [except in Eqs. (51) and(53)]. Furthermore, ξi (t ) (i = 1, 2, 3) is a zero mean andunit variance Gaussian white noise, independent for all thespheres:ξi (t ) = 0,(6)

ξi (t )ξ j (t  ) = δi j δ(t − t  ).

(7)

In contrast to Ref. [14], we do not consider hydrodynamic
interactions acting between different spheres in Eqs. (3)–(5).
Hence, the locomotion of a micromachine is not explicitly
taken into account in this paper, but such a treatment issufficient for the statistical analysis of the configurational
properties. Corrections due to hydrodynamic interactions will
briefly be discussed in Sec. VI.
In the following argument, it is convenient to introduce thetwo spring extensions with respect to :
r12 = x2 − x1 − ,

II. THERMALLY DRIVEN THREE-SPHERE
MICROMACHINE

We first explain the model of a thermally driven elastic micromachine that was introduced before by the present authors[14]. As schematically shown in Fig. 1, this model consists ofthree hard spheres connected by two harmonic springs. For
the sake of simplicity, we assume that the two springs are
identical, and the common spring constant and the naturallength are given by K and , respectively. Then the total elastic(2)

r23 = x3 − x2 − .

(8)

From Eqs. (3)–(5), we obtain the Langevin equations for r12 (t )and r23 (t ) as [31,32]

022607-2

KK
2T12 1/2
ξ12 ,ṙ12 = − r12 + r23 +
ζ12
ζ2
ζ12
1/2

K
K
2T23
ṙ23 = r12 −
r23 +
ξ23 .
ζ2
ζ23
ζ23

(9)
(10)

NONEQUILIBRIUM PROBABILITY FLUX OF A THERMALLY …

Here we have introduced the relevant effective friction coefficientζi ζ j,(11)
ζi j =
ζi + ζ jand the mobility-weighted average temperature
Ti j =

ζ j Ti + ζi Tj.
ζi + ζ j(12)

PHYSICAL REVIEW E 100, 022607 (2019)

Using Eqs. (18) and (19), one can calculate the threecorrelation functions r12 [ω]r12 [ω ], r23 [ω]r23 [ω ], andr12 [ω]r23 [ω ] in the frequency domain. Then the two vari22, σ23
and the covariance σ13 (or the equal-timeances σ12
correlation functions) are obtained as ∞
dω ∞ dω2
σ12 =r12 [ω]r12 [ω ]
−∞ 2π −∞ 2π
1(T12 + ζ12 ),K ∞
dω ∞ dω2
σ23 =r23 [ω]r23 [ω ]
−∞ 2π −∞ 2π=

The average temperatures T12 and T23 can be regarded as
two effective temperatures characterizing nonequilibrium behaviors of a thermally driven micromachine in the reduced
configuration space, i.e., in the (r12 , r23 ) space.
The definition of the effective temperature Ti j arises fromthe condition that the newly defined noises ξ12 (t ) and ξ23 (t ) in
Eqs. (9) and (10), respectively, satisfy the following statisticalproperties:ξ12 (t ) = ξ23 (t ) = 0,








ξ12 (t )ξ12 (t ) = δ(t − t ),
ξ23 (t )ξ23 (t ) = δ(t − t ),T2 ζ12 ζ23 1/2
ξ12 (t )ξ23 (t ) = −δ(t − t  ).
ζ2 T12 T23

(13)

=

(15)
(16)

III. STEADY STATE DISTRIBUTION FUNCTION

(22)
(23)

where we have used the notation(14)

It is important to note that the strength of the cross-correlationξ12 (t )ξ23 (t  ) in Eq. (16) is negative and its amplitude differsfrom unity. This noise property turns out to be important when
we discuss the Fokker-Planck equation for the probability
distribution function in Sec. IV.

ζ12 ζ23 (T12 + T23 − 2T2 ).
(ζ12 + ζ23 ) ζ22 − ζ12 ζ23(24)

With these results, we construct a symmetric covariancematrix C defined by 2σ12 σ13(25)
C=
2 .
σ13 σ23
Then the inverse of the covariance matrix can be simply given
by

 2
1
−σ13
σ23
C−1 = 2 2
,
(26)
2
σ12
σ12 σ23 (1 − ρ 2 ) −σ13
where the correlation factor ρ is defined by

A. Covariance matrix

In this section, we shall investigate the conformational
distribution of r12 and r23 that obey the coupled Langevinequations given by Eqs. (9) and (10). For this purpose, we
introduce a Fourier representation of any function f (t ) as ∞dωf (t ) =f [ω]eiωt ,
(17)
2π
−∞
where ω is the frequency and f [ω] is the Fourier component.
We rewrite Eqs. (9) and (10) in terms of r12 [ω] and r23 [ω] andsolve for them. After some calculations, they become
2T12 1/2ζ2
1 iωζ2
r12 [ω] = −+ξ12 [ω]
b
K
ζ23
ζ12



2T23 1/2+
ξ23 [ω] ,
(18)
ζ23



2T23 1/2ζ2
1 iωζ2+
r23 [ω] = −ξ23 [ω]
b
K
ζ12
ζ23



2T12 1/2+
ξ12 [ω] ,(19)
ζ12
where the common quantity b in the denominators is




1
K
ζ22
ζ2 2
1
ω+
1−
.
b = ω − iζ2
+
K
ζ12
ζ23
ζ2
ζ12 ζ23

1(T23 + ζ23 ),K ∞
dω ∞ dωζ2
σ13 =r12 [ω]r23 [ω ] =
,
2π
2π
K
−∞
−∞
=

(21)

(20)

ρ=

σ13
ζ2
=
.
σ12 σ23[(T12 + ζ12 )(T23 + ζ23 )]1/2(27)

Notice that the absolute value of the correlation factor should
satisfy the condition |ρ| < 1. As explained in Ref. [24] andsummarized in the Appendix, the steady state covariancematrix can generally be obtained by solving the corresponding
Lyapunov equation [see Eq. (A4)].
When the three friction coefficients are all identical, i.e.,ζ1 = ζ2 = ζ3 , the above correlation factor reduces toρ=2(T1 + T3 − 2T2 ).
[(7T1 + 4T2 + T3 )(T1 + 4T2 + 7T3 )]1/2(28)Here, we see that ρ generally vanishes when T1 + T3 = 2T2 ,i.e., the temperature of the middle sphere coincides with theaverage temperature between the first and the third spheres.
Obviously, ρ vanishes in thermal equilibrium, T1 = T2 = T3 ,although ρ = 0 does not mean that the micromachine is inthermal equilibrium.
B. Distribution function

Next we consider the steady state probability distributionfunction p(r), where r = (r12 , r23 )T and “T” indicates thetranspose. Owing to the reproductive property of Gaussian
distributions [23], p(r) should also be a Gaussian function for

022607-3

SOU, HOSAKA, YASUDA, AND KOMURAPHYSICAL REVIEW E 100, 022607 (2019)

FIG. 2. Dimensionless probability distribution function p(r)2
given by Eq. (30) as a function of r12 / and r23 /. The parametersare (a) ζ1 = ζ2 = ζ3 , τ1 = 1/900, τ2 = 41/900, and τ3 = 81/900(these temperatures satisfy τ1 + τ3 = 2τ2 ); (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 thesame 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)].

the present linear problem and is given by
1p(r) = N exp − rT C−1 r ,2(29)

where N is the normalization factor. Using Eq. (26) for C−1 ,we can write the explicit form of the distribution function asp(r) =

1
2π σ12 σ23 (1−ρ 2 )1/2
 2

2 
r12
1
r12 r23 r23, (30)
× exp −−2ρ+ 2
2
2(1 − ρ 2 ) σ12
σ12 σ23 σ23

together with Eqs. (21), (22), and (27). Despite the fact thatthe micromachine is out of equilibrium, the extensions of
the two springs obey the Boltzmann-type distribution thatis characterized by the effective temperatures T12 and T23[31]. Our further analysis crucially depends on the result of
Eq. (30).
Let us introduce the dimensionless temperature parameter
τi =

2Ti
,
K2

(31)

which is the ratio between the thermal energy of each sphereand the spring elastic energy (recall kB = 1). In Fig. 2, weplot the dimensionless steady state distribution function p(r)in Eq. (30) as a function of r12 / and r23 / for different
parameter combinations. The absolute value of the friction
coefficient is not required for these plots because we are focusing only on the stationary state distribution. The parametersin Fig. 2(a) are ζ1 = ζ2 = ζ3 , τ1 = 1/900, τ2 = 41/900, andτ3 = 81/900. Notice that these temperatures satisfy τ1 + τ3 =2τ2 and hence ρ = 0 according to Eq. (28). This means thatthe distribution function p(r) is simply a product of two
Gaussian functions, and there is no correlation between r12
and r23 for these temperatures. The distribution function iselongated in the r23 direction because τ1 < τ3 .
The parameters in Fig. 2(b) are ζ1 = ζ2 = ζ3 , τ1 = τ3 =25/900, and τ2 = 41/900. Here the temperatures of the first
and the third spheres are identical although the correlation factor obtained from Eq. (28) is negative, ρ < 0. The correlation
between r12 and r23 causes a tilt of the elongated distribution
and the tilt direction has a negative slope. Here the wholeplot is symmetric with respect to the straight line r12 = r23 .
In spite of the correlation between r12 and r23 , the systemis thermally balanced because the effective temperatures T12
and T23 are identical. As we mention in the next section, themicromachine behaves as if it were in thermal equilibrium inthe reduced configuration space.
Keeping the temperature parameters τi the same as inFigs. 2(a) and 2(b), we introduce asymmetries in the friction
coefficients, e.g., ζ1 /ζ2 = 0.5 and ζ3 /ζ2 = 5, in Figs. 2(c)and 2(d), respectively. Then the correlation factor shouldbe estimated according to Eq. (27). In Fig. 2(c), there is anegative correlation between r12 and r23 owing to the different
friction coefficients. On the other hand, Fig. 2(d) is no longer
symmetric with respect to the line r12 = r23 . Since T12 = T23
in this case [see Eq. (12)], the system is not thermally balanced
and exhibits nonequilibrium behaviors.
So far, we have calculated the steady state conformational
distribution function p of a thermally driven micromachine.
We have shown that it is given by a Gaussian function
characterized by the covariance matrix C. In the next section,
we shall calculate the probability flux by using the obtained
distribution function.
IV. PROBABILITY FLUX AND RELATED QUANTITIES
A. Probability flux

Next, we calculate the probability flux that can be used tocharacterize the nonequilibrium properties of a three-spheremicromachine. The Fokker-Planck equation corresponding toEqs. (9) and (10) can be written for the time-dependent conformational probability distribution function p(r, t ) as [33,34]ṗ = −∇ · j,(32)

expressing the fact that probability is neither created nor
annihilated. In the above, the two-dimensional nabla operator
indicates ∇ = (∂/∂r12 , ∂/∂r23 )T in the configuration space
and the vector j = ( j12 , j23 )T is the two-dimensional probability flux given by
j = Arp − D∇ p,
where the matrices A and D are
−K/ζ12A=K/ζ2and022607-4

T /ζ
D = 12 12−T2 /ζ2

K/ζ2−K/ζ23(33)
−T2 /ζ2,
T23 /ζ23

(34)

(35)

NONEQUILIBRIUM PROBABILITY FLUX OF A THERMALLY …

respectively. Here A describes the linear deterministic dynamics (without noise) in Eqs. (9) and (10), whereas D is called
the diffusion matrix (see also the Appendix). The nonzero offdiagonal components of D originate from the noise correlation
shown in Eq. (16).
In a steady state, i.e., ṗ = 0, the probability flux j(r) canbe obtained from Eq. (33) by using the Gaussian probability
distribution function p in Eq. (30): 
1
T12 r12r23
r12+
j12 (r) = K
−
2
2
ζ2
ζ12
(1 − ρ ) ζ12 σ12
T2 r12
T2 r23ρ
T12 r23−p, (36)+−2σ12 σ23ζ2
ζ12
ζ2 σ23

PHYSICAL REVIEW E 100, 022607 (2019) 1
T2 r12
r12
r23
+
−
−
j23 (r) = K
2
2
ζ2
ζ23
(1 − ρ )ζ2 σ12
T2 r23
T23 r23ρ
T23 r12+p.
+−
2
σ12 σ23ζ2
ζ23
ζ23 σ23

One can confirm that this probability flux is divergence-free,i.e., ∇ · j = 0, so that the steady state condition is satisfied.
B. Frequency matrix

The above probability flux j can be conveniently expressedin terms of a frequency matrix  as [24]j = rp,where the explicit expression of  is given by


⎛
T12
1
ρT2
K
+
+
−
⎜ ζ12
(1 − ρ 2 ) ζ12 σ12 2
σ12 σ23 ζ2


=⎜
⎝ K
T2
1
ρT23
−
+
ζ2
(1 − ρ 2 ) ζ2 σ12 2
σ12 σ23 ζ23
Since j is divergence-free,  should be a traceless matrix [25].
Using Eqs. (11) and (12), we indeed find that  is traceless.
In the Appendix, we generally show that the frequency matrix
 can be expressed as  = A + DC−1 within the matrix
formulation [see Eq. (A6)].
When the friction coefficients are identical, ζ1 = ζ2 =
ζ3 = ζ , the frequency matrix is simplified to=K (T1 − T3 )ζc
2(T1 + T3 − 2T2 )×
T1 + 4T2 + 7T3

−(7T1 + 4T2 + T3 ),−2(T1 + T3 − 2T2 )

(37)

(38)

 ⎞

K
T2
1
ρT12
−
+
ζ2
(1 − ρ 2 ) ζ2 σ23 2
σ12 σ23 ζ12 ⎟
⎟

⎠.
K
T23
1
ρT2
−
+
+
ζ23
(1 − ρ 2 ) ζ23 σ23 2
σ12 σ23 ζ2

by√i 3K (T1 − T3 ).
γ =±ζ c1/2(39)(42)

Since these eigenvalues are purely imaginary, the probability
current in the configuration space is rotational [24,27]. In theabove, the essential frequency scale is set by γ ∼ K/ζ whichcharacterizes the speed of the rotational motion.
In Fig. 3, we plot the dimensionless probability flux vectorj(r) in Eq. (38) as a function of r12 / and r23 / by using thesame parameter values as in Fig. 2. Note that the color scale(40)

where the quantity c in the denominators is
c = T12 + 16T1 T2 + 14T1 T3 + 16T2 T3 + T32 .

(41)

Here, we explicitly see that  is proportional to T1 − T3 ,and hence the probability flux j vanishes when T1 = T3 forany T2 of the middle sphere. Indeed, j = 0 or  = 0 is anecessary and sufficient condition for physical situations inwhich detailed balance is satisfied [24]. When the system
is in thermal equilibrium and detailed balance holds, theprobability of a transition between any two states is the sameas the probability of the reverse transition. As a result, a net
probability flux does not exist. In the Appendix, we generally
show that the necessary and sufficient condition for detailedbalance is given by the commutation relation AD − DAT = 0[see Eq. (A7)].
When T1 = T3 , on the other hand, we have nonzero frequency matrix . Then detailed balance is violated and themicromachine is in a nonequilibrium steady state. Hence thefrequency matrix is an important measure to quantify the outof-equilibrium properties of a thermally driven micromachine.
When T1 = T3 , the two eigenvalues of  in Eq. (40) are given

FIG. 3. Dimensionless probability flux vector j(r)ζ2 /K given
by Eq. (38) as a function of r12 / and r23 /. The parameter values
are the same as in Fig. 2. There is no probability flux in (b). The
color scale indicates the magnitude of the flux vector j.

022607-5

SOU, HOSAKA, YASUDA, AND KOMURAPHYSICAL REVIEW E 100, 022607 (2019)indicates the magnitude of the flux vector. In Figs. 3(a), 3(c)and 3(d), clockwise flux loops can be clearly seen. When thefriction coefficients are identical as in Fig. 3(a), the existenceof a flux loop is consistent with a directed motion of a threesphere micromachine [14], as we further discuss in Sec. V.
In Fig. 3(b) with τ1 = τ3 , however, the probability flux
does not exist simply because  = 0. Asymmetries both in
the temperatures and the friction coefficients lead to a flux
loop in Fig. 3(c). The presence of a flux loop in Fig. 3(d) evenfor τ1 = τ3 is due to the asymmetry in the friction coefficients.
Although the case with different friction coefficients was notstudied in Ref. [14], we expect that a finite flux loop leads toa directed motion of a three-sphere micromachine despite thebalanced temperature distribution.
C. Flux rotor

To further characterize the strength of the probability flux
loop in terms of a scalar quantity, we consider the following
flux rotor [35]:
∂
∂
j23 −
j12
∂r12
∂r23


2
2
+ w23 r23
+ w13 r12 r23
p w0 + w12 r12=,4 σ 4 (1 − ρ 2 )2
ζ2 ζ12 ζ23 σ1223s(r) =

where w0 , w12 , w23 , and w13 are
 2


2 2
2
w0 = σ12σ23 (1 − ρ 2 ) ζ12 ζ23 T2 σ12− σ23
+ ζ2 σ12 σ23 ρ(ζ23 T12 − ζ12 T23 ) ,
w12 = ζ2 σ12 σ23 ρ(ζ12 T23 − ζ23 T12 )2− Kζ23 σ12σ23 (1 − ρ 2 )(ζ12 σ23 − ζ2 σ12 ρ)
 2
 2
2
σ23 ,
− ρ 2 σ12+ ζ12 ζ23 T2 σ23
w23 = ζ2 σ12 σ23 ρ(ζ12 T23 − ζ23 T12 )2− Kζ12 σ12 σ23(1 − ρ 2 )(ζ2 σ23 ρ − ζ23 σ12 )
 2
 2
2
− ζ12 ζ23 T2 σ12σ12 ,− ρ 2 σ23 22
w13 = 2ζ12 ζ23 T2 ρ σ12 − σ2322+ Kζ2 σ12 σ23 (1 − ρ 2 ) ζ12 σ23− ζ23 σ12+ ζ2 σ12 σ23 (1 + ρ 2 )(ζ23 T12 − ζ12 T23 ) σ12 σ23 .

(43)

FIG. 4. Dimensionless flux rotor s(r)ζ2 2 /K given by Eq. (43) asa function of r12 / and r23 /. The parameter values are the same as
in Fig. 2. There is no flux rotor in (b).

s(r = 0). This quantity, denoted as s0 , is given bys0 =1
3 σ 3 (1 − ρ 2 )3/2
2π ζ2 ζ12 ζ23 σ12
23

 2



2
+ ζ2 σ12 σ23 ρ(ζ23 T12 − ζ12 T23 ) ,× ζ12 ζ23 T2 σ12
− σ23
(48)
(44)

(45)

(46)

because only the term proportional to w0 in Eq. (43) remains.
When the friction coefficients are identical, ζ1 = ζ2 = ζ3 = ζ ,Eq. (48) further reduces to√
16 3K 2 (T1 + T2 + T3 )(T1 − T3 )
s0 =
,
(49)
π ζ c3/2
where c is given by Eq. (41). Here the sign of s0 is purely
determined by the temperature difference T1 − T3 . In Fig. 5,
we plot the dimensionless s0 as a function of T1 /T2 and T3 /T2
by using a color representation.

(47)

In Fig. 4, we plot the dimensionless flux rotor s(r) in
Eq. (43) as a function of r12 / and r23 / by using the same
parameter values as in Fig. 2. In Figs. 4(a), 4(c) and 4(d),the flux rotor s takes negative values around the origin of theconfiguration space, r = 0, whereas it takes positive values in
the outer regions. In Fig. 4(a), the distribution of s is elongatedin the r23 direction, while there is no correlation between r12
and r23 . The flux rotor vanishes in Fig. 4(b) simply because
 = 0 for τ1 = τ3 . When the friction coefficients are different
as in Figs. 4(c) and 4(d), the flux rotor exhibits a negative
correlation.
As considered in Ref. [35], one can focus on the strengthof the flux rotor at the origin of the configuration space, i.e.,√
FIG. 5. Dimensionless flux rotor π ζ T2 s0 /(16 3K 2 ) given by
Eq. (49) as a function of T1 /T2 and T3 /T2 . The flux rotor is evaluatedat r = 0 when the friction coefficients are identical. Different curves
are the contour lines for the corresponding values.

022607-6

NONEQUILIBRIUM PROBABILITY FLUX OF A THERMALLY …

In this section, starting from the Fokker-Planck equation,
we have obtained the steady state probability flux j thatcan be expressed in terms of the frequency matrix . The
eigenvalues γ of the frequency matrix are purely imaginary
and proportional to the temperature difference. We have also
calculated the flux rotor s as a scalar quantity to characterize
the scale of nonequilibrium of a micromachine. In the next
section, we shall discuss the relation between these quantities
and the average velocity of a stochastic micromachine.
V. AVERAGE VELOCITY OF A MICROMACHINEAs mentioned before, the main purpose of this paper is tounderstand the physical mechanism that underlies the locomotion of a thermally driven micromachine within the nonequilibrium statistical mechanics. When the friction coefficientsare identical, the average velocity is given by Eq. (1) which is
proportional to T1 − T3 . This proportionality appears in several scalar quantities such as the eigenvalues of the frequency
matrix, γ [see Eq. (42)], or the flux rotor at the origin, s0 [seeEq. (49)]. Hence one can construct linear relations between
the average velocity and these scalar quantities.
For example, the absolute value of the average velocity canbe expressed in terms of γ as|V | =c1/2 
a|γ |,√16 32 K(50)

where we have used Stokes’ law, ζ = 6π ηa, for the friction
coefficient of a hard sphere, and γ  indicates the imaginary
part of γ . For the sake of completeness, we rewrite the above
expression by recovering the Boltzmann constant and using
Eq. (41) as follows:a|V | = √16 321/2
kB T12 + 16T1 T2 + 14T1 T3 + 16T2 T3 + T32|γ  |.
×
K
(51)


As mentioned before, γ has the dimension of frequency
and it determines the speed of the rotational motion of theprobability flux in the configuration space.
An interesting interpretation of the above expression can bemade by comparing it to the general form for the average velocity of a three-sphere swimmer. When the two arms undergo
a prescribed deterministic motion, as in the original NajafiGolestanian model [4], the average velocity is given by [5]V =G(r12 ṙ23 − ṙ12 r23 ),2(52)

where G ∼ a/2 is a geometric factor that depends only on
the structural parameters of a swimmer, and the averaging is
performed by time integration in a full cycle. In the above
form, the averaging part is proportional to the enclosed area(per unit time) that is swept in a full cycle in the configuration
space. Equation (52) can be regarded as a mathematicalrepresentation of Purcell’s scallop theorem for a three-sphereswimmer in low-Reynolds-number fluids [2,3].
For a thermally driven micromachine, V  is proportional
to a/2 as in Eq. (52). Moreover, c1/2 in Eq. (50) can be

PHYSICAL REVIEW E 100, 022607 (2019)

regarded as overall thermal energy, and c1/2 /K essentiallyrepresents the explored area in the configuration space by random motions of the spheres. This interpretation is reasonablebecause one can easily confirm the relation c1/2 /K ∼ σ12 σ23 ,where σ12 and σ23 are the variances in Eqs. (21) and (22),respectively.
While the timescale in the Najafi-Golestanian model isgiven by the frequency of the periodic arm motions [5], thecorresponding timescale for a thermally driven micromachineis set by the eigenvalue γ of the frequency matrix [seeEq. (42)]. Hence Eq. (50) or Eq. (51) provides us with anessential understanding concerning the locomotion of a thermally driven micromachine. It is interesting to note that themain concept of Purcell’s scallop theorem can be generalized
for thermally driven micromachines that undergo random
motions rather than deterministic cyclic motions.
In order to take into account the sign of the average velocity, namely, the direction of the locomotion, it is convenient
to use the flux rotor at the origin s0 [see Eq. (49)]. After
recovering the Boltzmann constant, the average velocity canalso be written asπaV  = −√256 323/2
kB2 T12 + 16T1 T2 + 14T1 T3 + 16T2 T3 + T32
s0 .
×
K 2 (T1 + T2 + T3 )(53)
In the above, the sign of s0 determines the direction of thelocomotion. For example, when T1 < T3 as in Fig. 4(a) orin Fig. 5, we have s0 < 0 and V  > 0. Although Eqs. (51)and (53) are essentially equivalent, we consider that Eq. (51)gives more physical insights to understand the locomotion ofa thermally driven micromachine.
VI. SUMMARY AND DISCUSSION

In this paper, we have discussed the nonequilibrium statistical mechanics of a thermally driven micromachine madeof three spheres and two harmonic springs as previously
proposed by the present authors [14]. First, we have calculated
the steady state conformational distribution function of sucha micromachine and showed that it is given by a Gaussianfunction characterized by the covariance matrix [see Eq. (30)].
Using this distribution function, we have obtained the steadystate probability flux of a micromachine.
The distribution function can be expressed in terms of a
traceless frequency matrix [see Eqs. (39) and (40)]. When thefriction coefficients are all identical, we have shown thatthe eigenvalues of the frequency matrix are proportional to thetemperature difference between the first and the third spheres[see Eq. (40)]. Moreover, the scale of nonequilibrium of a
micromachine can quantitatively be characterized by the flux
rotor [see Eq. (43) or Eq. (48)]. As one of the main results ofthis paper, the average velocity of a thermally driven machineis expressed in terms of the eigenvalue of the frequency matrix[see Eq. (50)]. This expression allows us to generalize theconcept of Purcell’s scallop theorem that is also applicablefor thermally driven stochastic micromachines.

022607-7

SOU, HOSAKA, YASUDA, AND KOMURAPHYSICAL REVIEW E 100, 022607 (2019)

An interesting situation to be discussed in more detail is thecase when ζ1 = ζ2 = ζ3 and T1 = T3 = T2 . Microscopically,such a micromachine is in a nonequilibrium steady statebecause the temperature of the second sphere is different
from those of the other two spheres. However, within theconfiguration space (r12 , r23 ), the two relevant average temperatures are equal, i.e., T12 = T23 . Hence the overall stateof a micromachine is effectively in thermal equilibrium anddetailed balance is satisfied, i.e.,  = 0. This is the reason
why the average velocity in Eq. (1) vanishes when T1 = T3
regardless of the value of T2 .
It is known that for some nonequilibrium systems, broken
detailed balance does not need to be apparent at larger scales,
and they can regain thermodynamic equilibrium when thesystem is coarse grained [36,37]. In our analysis, the reduction of the configuration space has been implicitly assumed
because we have focused only on the distribution of r12 andr23 , while the center of mass motion of a micromachine has
been neglected. Within this level of description, an elastic
three-sphere micromachine behaves as if it were in thermal
equilibrium when T1 = T3 .
In our previous paper [14], we argued that the averagevelocity of a micromachine can be related to the ensemble average of heat flows in a steady state. According to “stochastic
energetics,” the heat gained by the ith sphere per unit time isexpressed as [20–22]
Q̇i = ζi (−ẋi + ξi )ẋi ,(54)where ẋi (t ) and ξi (t ) are given in Eqs. (3)–(5). When thefriction coefficients are identical, we showed that the averagevelocity can also be expressed in terms of the lowest-order
average heat flows asa
V  =(Q̇3  − Q̇1 ).
(55)
8K2
This relation states that the average velocity is determined by
the net heat flow between the first and the third spheres. Our
result in this paper indicates that the net heat flow between thefirst and the third spheres is also proportional to the eigenvalueof the frequency matrix γ or the flux rotor s0 , as discussed inSec. V.
We mention here that our model of a three-sphere micromachine has a similarity to that of two overdamped, tethered
spheres coupled by a harmonic spring and also confined
between two walls [27,28]. Because these two spheres arein contact with heat baths having different temperatures, thesystem can be driven out of equilibrium. They numericallyshowed that displacements obey a Gaussian distribution andalso found probability flux loops that demonstrate the broken
detailed balance [28]. In Ref. [27], an analytical expression of
the frequency matrix for this two-sphere model was shown tobe proportional to the temperature difference.
Clearly, the two displacements r12 and r23 in Eq. (8) correspond to the sphere positions in their model. It is interesting
to note, for example, that our result of the frequency matrix in
Eq. (40) reduces to Eq. (42) in Ref. [27] when T2 = 0. WhenT2 = 0, however, the presence of the middle sphere changes
the structure of the frequency matrix as we have shown in
Eq. (39) or Eq. (40). Moreover, the important outcome ofthis work is the relation between the average velocity andthe eigenvalue of the frequency matrix for a three-spheremicromachine [see Eq. (50) or Eq. (51)]. Notice that a twosphere micromachine in a viscous fluid cannot have a directedmotion even if the temperatures are different [14].
As mentioned before, we have neglected long-ranged hydrodynamic interactions acting between different spheres. In
our previous paper [14], we explicitly took into account hydrodynamic interactions when the friction coefficients are allidentical. If hydrodynamic interactions are taken into accountin the present analysis, the variances and the covariance in
Eqs. (21)–(23) are modified in nonequilibrium situations.
Such hydrodynamic corrections should be proportional to a/
within the lowest-order expansion. Moreover, such corrections should vanish in thermal equilibrium, i.e., T1 = T2 = T3
because hydrodynamic interactions should not affect equilibrium statistical properties.
In the presence of hydrodynamic interactions, a thermally
driven elastic micromachine can undergo a directional motion[14]. The presence of the middle sphere is essential for a directional motion because the hydrodynamic interactions among
the three spheres are responsible for it. Such a locomotion
should be distinguished from the traditional thermophoresis(Ludwig-Soret effect) in which a temperature gradient in an
external fluid induces a directed motion of suspended particles[38]. In our model, the locomotion of a micromachine ispurely induced by nonequilibrium fluctuations of internal
degrees of freedom. Here, the three spheres are in thermal
equilibrium with independent heat baths having different temperatures, which is different from a situation where a temperature gradient is externally imposed in a surrounding fluid[38]. For example, we showed before, both analytically andnumerically, that a two-sphere elastic micromachine cannot
move even if the temperatures are different [14]. This clearly
indicates that the locomotion of a thermally driven micromachine cannot be explained within the standard thermophoresis,
although it is closely related to Purcell’s scallop theorem for
microswimmers under the force-free condition.
In the present work, we have shown that asymmetries in
the friction coefficients lead to a finite probability flux loop
[see Fig. 3(d)] or flux rotor [see Fig. 4(d)], anticipating a
locomotion of an asymmetric micromachine in a viscous fluid.
This prediction will be investigated in the future by performing numerical simulations of the coupled stochastic equations,as given by Eqs. (3)–(5), in the presence of hydrodynamic
interactions. Notice that analytical treatment of a case withdifferent sphere sizes is more difficult because one needs totake into account higher-order contributions in r12 and r23 aswell as in the sphere size a.
So far, various concepts have been proposed to quantitatively discuss whether a steady state is in thermal equilibrium
or not [37]. One of the most promising methods is to search
for the violation of the fluctuation-dissipation relation that is
guaranteed in thermal equilibrium situations [39,40]. However, it is not always easy to perform two separate measuresof the correlation function and of the response function forthe same system. Moreover, the measurement of the response
function can often be intrinsically invasive and is not suitable for biological systems. Moreover, the observation of
non-Gaussian distribution fluctuations is not a proof for thenonequilibrium steady state [28,41]. Since the emergence of022607-8

NONEQUILIBRIUM PROBABILITY FLUX OF A THERMALLY …

probability flux loops is a direct verification of a nonequilibrium steady state, it can be a powerful method to study
systems that are driven out of equilibrium. We consider thata thermally driven three-sphere micromachine is an excellent
example to show the usefulness of such an analysis.
ACKNOWLEDGMENTS

We thank T. Kato and M. Doi for useful discussions. We
also thank the anonymous referee for the useful suggestions.
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)
from the JSPS. S.K. acknowledges support by a Grant-in-Aid
for Scientific Research (C) (Grants No. 18K03567 and No.
19K03765) from the JSPS.
APPENDIX: MATRIX REPRESENTATION OF
STOCHASTIC DYNAMICAL SYSTEMS

Following closely the argument in Refs. [24,25], we shallbriefly review the general formulation of stochastic dynamicalsystems using a matrix representation. Let us start from thelinear stochastic Langevin modelṙ = Ar + ξ,

(A1)

where r is an N-dimensional vector characterizing the state of
the system, A is an (N × N )-dimensional matrix describing
the linear deterministic dynamics, and ξ is an N-dimensional
vector representing the noise forcing. The time-dependent
covariance matrix is introduced by C = r(t )rT (t  ). Then the
diffusion matrix D is defined by the relationξ(t )ξ T (t  ) = 2Dδ(t − t  ).

PHYSICAL REVIEW E 100, 022607 (2019)

Hereafter, we use the same notation C for the covariance matrix that satisfies the Lyapunov equation. It should be emphasized that the Lyapunov equation holds in both in-equilibriumand out-of-equilibrium situations. Hence it can be regarded asa generalized fluctuation-dissipation relation connecting thefluctuations (C) and the deterministic dissipation (A).
For linear systems with Gaussian noise, the steady state
probability distribution function p is a Gaussian function as in
Eq. (29). In this case, the probability flux defined by Eq. (33)
becomesj = Arp − D∇ p = (A + DC−1 )rp.

Hence the frequency matrix  introduced through the relationj = rp [see Eq. (38)] is given by = A + DC−1 .
(A6)When  = 0 and hence detailed balance holds, we have
A + DC−1 = 0. After the elimination of C in the Lyapunov
equation, we then haveAD − DAT = 0.

(A7)This commutation relation holds if and only if the system is inthermal equilibrium and satisfies detailed balance. Moreover,whether or not detailed balance is satisfied is coordinate
invariant [24].
When  = 0 and hence detailed balance is broken, thesystem is in nonequilibrium steady state situations. In thiscase, one can show that the following two relations hold:AC + D = C,(A2)

From Eq. (A1), one can show that the time evolution of thecovariance matrix is given by(A5)

CAT + D = −C.

(A8)

(A4)

Obviously, the sum of these two relations gives again the
Lyapunov equation in Eq. (A4). One can also show that theflow field, as expressed by r, is perpendicular to the gradientof the probability distribution, i.e., (r) · ∇ p = 0 [27].
Last but not least, the matrices C in Eq. (25), A in Eq. (34),
D in Eq. (35), and  in Eq. (39) obtained for a three-spheremicromachine satisfy all the relations in this Appendix.

[1] E. Lauga and T. R. Powers, Rep. Prog. Phys. 72, 096601 (2009).
[2] E. M. Purcell, Proc. Natl. Acad. Sci. USA 94, 11307 (1997).
[3] E. Lauga, Soft Matter 7, 3060 (2011).
[4] A. Najafi and R. Golestanian, Phys. Rev. E 69, 062901 (2004).
[5] R. Golestanian and A. Ajdari, Phys. Rev. E 77, 036308(2008).
[6] M. Leoni, J. Kotar, B. Bassetti, P. Cicuta, and M. C.
Lagomarsino, Soft Matter 5, 472 (2009).
[7] G. Grosjean, M. Hubert, G. Lagubeau, and N. Vandewalle,
Phys. Rev. E 94, 021101(R) (2016).
[8] G. Grosjean, M. Hubert, and N. Vandewalle, Adv. Colloid
Interface Sci. 255, 84 (2018).
[9] K. Yasuda, Y. Hosaka, M. Kuroda, R. Okamoto, and S. Komura,
J. Phys. Soc. Jpn. 86, 093801 (2017).

[10] M. Kuroda, K. Yasuda, and S. Komura, J. Phys. Soc. Jpn. 88,054804 (2019).
[11] J. Dunkel and I. M. Zaid, Phys. Rev. E 80, 021903 (2009).
[12] J. Pande and A.-S. Smith, Soft Matter 11, 2364 (2015).
[13] J. Pande, L. Merchant, T. Krüger, J. Harting, and A.-S. Smith,
New J. Phys. 19, 053024 (2017).
[14] Y. Hosaka, K. Yasuda, I. Sou, R. Okamoto, and S. Komura,
J. Phys. Soc. Jpn. 86, 113801 (2017).
[15] M. O. Magnasco, Phys. Rev. Lett. 71, 1477 (1993).
[16] R. D. Astumian and M. Bier, Phys. Rev. Lett. 72, 1766 (1994).
[17] J. Rousselet, L. Salome, A. Ajdari, and J. Prost, Nature(London) 370, 446 (1994).
[18] A.-Y. Jee, S. Dutta, Y.-K. Cho, T. Tlusty, and S. Granick, Proc.
Natl. Acad. Sci. USA 115, 14 (2018).

Ċ = AC + CAT + 2D.

(A3)

When the system is in a steady state, i.e., Ċ = 0, the covariance matrix must obey the following Lyapunov equation:
AC + CAT + 2D = 0.

022607-9

SOU, HOSAKA, YASUDA, AND KOMURAPHYSICAL REVIEW E 100, 022607 (2019)[19] K. K. Dey, Angew. Chem., Int. Ed. 58, 2208 (2019).
[20] K. Sekimoto, J. Phys. Soc. Jpn. 66, 1234 (1997).
[21] K. Sekimoto, Prog. Theor. Phys. Suppl. 130, 17 (1998).
[22] K. Sekimoto, Stochastic Energetics (Springer, Berlin/
Heidelberg, 2010).
[23] N. van Kampen, Stochastic Processes in Physics and Chemistry(North-Holland, Amsterdam, 1992).
[24] J. B. Weiss, Tellus A 55, 208 (2003).
[25] J. B. Weiss, Phys. Rev. E 76, 061128 (2007).
[26] A. Ghanta, J. C. Neu, and S. Teitsworth, Phys. Rev. E 95,
032128 (2017).
[27] F. S. Gnesotto, F. Mura, J. Gladrow, and C. P. Broedersz, Rep.
Prog. Phys. 81, 066601 (2018).
[28] C. Battle, C. P. Broedersz, N. Fakhri, V. F. Geyer, J. Howard,
C. F. Schmidt, and F. C. MacKintosh, Science 352, 604(2016).
[29] J. Gladrow, N. Fakhri, F. C. MacKintosh, C. F. Schmidt,
and C. P. Broedersz, Phys. Rev. Lett. 116, 248301(2016).
[30] J. Gladrow, C. P. Broedersz, and C. F. Schmidt, Phys. Rev. E
96, 022408 (2017).

[31] A. Y. Grosberg and J.-F. Joanny, Phys. Rev. E 92, 032118(2015).
[32] R. R. Netz, J. Chem. Phys. 148, 185101 (2018).
[33] H. Risken, The Fokker-Planck Equation: Methods of Solution
and Applications (Springer-Verlag, Berlin, 1984).
[34] R. Zwanzig, Nonequilibrium Statistical Mechanics (Oxford
University Press, New York, 2001).
[35] V. Dotsenko, A. Maciołek, O. Vasilyev, and G. Oshanin, Phys.
Rev. E 87, 062130 (2013).
[36] D. A. Egolf, Science 287, 101 (2000).
[37] J.-F. Rupprecht and J. Prost, Science 352, 514 (2016).
[38] H.-R. Jiang, N. Yoshinaga, and M. Sano, Phys. Rev. Lett. 105,
268302 (2010).
[39] M. Guo, A. J. Ehrlicher, M. H. Jensen, M. Renz, J. R. Moore,
R. D. Goldman, J. Lippincott-Schwartz, F. C. MacKintosh, and
D. A. Weitz, Cell 158, 822 (2014).
[40] H. Turlier, D. A. Fedosov, B. Audoly, T. Auth, N. S. Gov, C.
Sykes, J.-F. Joanny, G. Gompper, and T. Betz, Nat. Phys. 12,513 (2016).
[41] C. Bustamante, J. Liphardt, and F. Ritort, Phys. Today 58(7), 43(2005).

022607-10