PHYSICAL REVIEW RESEARCH
6
, 023195 (2024)
Theoretical limits of energy extraction in active fluids
Shahriar Shadkhoo
*
and Matt Thomson
California Institute of Technology, Pasadena, California 91125, USA
(Received 5 July 2023; revised 21 November 2023; accepted 26 April 2024; published 21 May 2024)
Active materials form a class of far-from-equilibrium systems that are driven internally and exhibit self-
organization, which can be harnessed to perform mechanical work. Inspired by experiments on synthetic active
polymer networks, in this paper we examine limits of work extraction from an active viscoelastic medium by
analyzing the transport of a particle. The active viscoelastic material possesses an equilibrium density where
the active and passive forces are balanced out. In a one-dimensional system, a gliding activation front (AF),
which converts a passive to an active medium, provides active energy at a constant rate, which is injected into
the system at one end and propagates to the other. We demonstrate that there exists a maximum velocity of the
AF, above which the activated region fails to deliver the transport power. We hypothesize and intuitively argue
based on the limit cases that the feasibility and the velocity of transport of the particle can be interpreted in
terms of the velocity of an equilibration domain wall of the field, which is set by two parameters: (i) a measure
of activity, and (ii) the viscoelastic timescale. The phase diagram is divided into “transport” and “no-transport”
sectors, namely, for any pair of the two parameters, there exists a threshold velocity of the AF above which
the particle transport becomes impossible. Constructing the phase diagram we find that (1) there are regions
of the phase diagram for which the threshold velocity of the AF diverges, and (2) larger viscoelastic timescale
makes the transport region more accessible and also increases the transport velocity therein. Furthermore, we
find that increasing the velocity of AF, results in larger extracted power but smaller transport coefficient; defined
as the ratio of the transport velocity and that of the AF. Our model provides a framework for understanding the
energetics of transport phenomena in biology and designing efficient mechanisms of transport in synthetic active
materials.
DOI:
10.1103/PhysRevResearch.6.023195
I. INTRODUCTION
For more than a century the problem of energy extraction
from nonequilibrium systems has been a subject of extensive
research, with significant implications in fundamental physics
and engineering. Thought experiments like Maxwell’s demon
and the generalizations thereof have established deep connec-
tions between information gain and work extraction [
1
]. Such
studies have also advanced our understanding of the optimal
protocols of energy extraction [
2
–
10
]. The role of feedback
loops in the process of information-to-energy conversion was
highlighted in an early instantiation of Maxwell’s demon sug-
gested by Szilárd [
11
]. This paradigm of energy extraction
at the microscopic scales leverages information to rectify the
fluctuations in the desired direction, hence decreasing entropy.
Energy can also be extracted at macroscopic scales from
systems that are externally supplied with energy—usually at
the boundaries. The injected energy excites the dynamical
modes of the system which may subsequently perform work.
*
shahriar@caltech.edu
Published by the American Physical Society under the terms of the
Creative Commons Attribution 4.0 International
license. Further
distribution of this work must maintain attribution to the author(s)
and the published article’s title, journal citation, and DOI.
Active materials constitute a class of far-from-equilibrium
systems that are driven internally via injection of energy
at microscopic scales and exhibit self-organization—a com-
monplace phenomenon in biology. In spite of its prevalence,
the mechanisms of energy propagation and conversion in
active systems is yet to be elucidated. Unraveling such
mechanisms would provide insights into the efficiency of
biological pathways and the barriers that impede the con-
version of active energy into mechanical work [
12
]. Besides
natural realizations, synthetic active systems allow us to pro-
gram novel states of matter through harnessing the power
of self-organization, and designing optimal control protocols
[
13
–
23
].
In this article we propose a model of an active system with
a built-in mechanical feedback which can be utilized to extract
energy. The feedback mechanism operates on the basis of
activity-induced contraction in viscoelastic active materials,
which in turn amplifies local activity and the resultant contrac-
tion; Fig.
1(a)
. Coupling a particle to the active medium and
studying the transport of the particle, we investigate the limits
of energy extraction from the
boundary
of scalar viscoelastic
active systems; see Fig.
1(b)
.
The limits of energy extraction are determined by the prop-
erties of the system at hand. Here we consider a continuum
model, unlike most studies where the stochastic thermody-
namics of the microscopic degrees of freedom is considered
[
24
–
29
]. A continuum description of viscoelasticity comprise
2643-1564/2024/6(2)/023195(16)
023195-1
Published by the American Physical Society
SHAHRIAR SHADKHOO AND MATT THOMSON
PHYSICAL REVIEW RESEARCH
6
, 023195 (2024)
FIG. 1. (a) The positive and negative feedback loops between
strain and the two competing contractile and extensile components
of stress, which originate from the active and passive components,
respectively. (b) Transport of a particle coupled to an active string.
The activation front generates stress gradient and pulls on the parti-
cle. The densities of the passive and active components are denoted
by
φ
p
and
φ
a
, respectively, the sum of which equals
, a constant.
an elastic and a viscous component expressed in terms of
the strain and stress fields. The elasticity and viscosity to-
gether define a timescale that characterizes the timescale over
which the dominant mode transitions from viscous to elastic,
or vice versa. In our model—inspired by the dynamics of
active biopolymer networks—activity, which introduces the
source of the to-be-extracted energy, is manifested in the form
of density-dependent contractile stresses (see Sec.
II
). Stress
gradients provide forces that can potentially perform work on
a second system that is coupled to the active one. In a system
with statically defined boundaries, i.e., no material exchange,
the energy extraction (and motion) cannot be sustained in the
steady-state sense. The active energy excites the dynamics,
which is subsequently damped out towards an equilibrium
state, in which the active energy is constantly converted into
the elastic energy of the equilibrium state. Should the sys-
tem be provided with constant input of fresh active energy,
however, the active energy can potentially provide the power
required for sustained dynamics.
Steady-state transport of an external particle in a viscous
medium exemplifies the extraction of energy. Transport phe-
nomena require net (thermodynamic) forces like pressure
gradients. In isolated passive systems, transport of the con-
stituents is driven by minimization of the free energy towards
the thermodynamic equilibrium [
30
]. Active materials utilized
bulk energy to perform work at large scales [
31
]. Important
examples of active transport in biology include cytoplasmic
transport of proteins and ions [
12
,
27
,
32
–
36
]. The study of
active transport has been mostly with a few exceptions, e.g., in
Refs. [
24
,
37
]—limited to self-propelled particles in a passive
medium [
38
–
43
]. An alternative scenario is for the particle to
be coupled to a reservoir of active constituents that provides
the transport energy.
Here we aim to understand the properties of the active
system that determine the limits of energy extraction. Our
TABLE I. Model parameters and variables with their definitions.
Parameter
Definition
α
1
,α
2
Coefficients of stress vs density
η
Viscosity of the fluid component
κ
Elasticity of the rigid component
τ
K
=
κ
−
1
Kelvin-Voigt viscoelastic time
γ
Drag coefficient of the field
φ
i
Initial density of activated region
V
AF
Velocity of the activation front
g
0
Field-particle coupling constant
Drag coefficient of the particle
Variable
Definition
,σ
Stress field of active region
E
,
ε
Strain field of active region
φ
a
Active density field
φ
e
Equilibrium density
φ
b
Boundary density
v
b
Boundary velocity
V
DW
Velocity of the domain wall
model is motivated by an experimental work in which motility
is induced in a mixture of motor proteins and microtubule
filaments. The motor proteins are activated upon shining light,
form cross-links, and walk along the filaments [
14
,
44
]. A
microtubule dense aster resembling the external rigid particle
is transported as the activation front sweeps across the system
at a fixed velocity. An upper bound on the transport velocity
is observed when the velocity of the activation front exceeds
a threshold value. Even though in the previous experiments
the mechanical properties of the active networks (such as
elasticity, viscosity, and the magnitude of active stress) vary
over narrow ranges of parameters, the synthetic nature of such
systems provides the possibility of exploring their behavior
beyond the currently available ones. Therefore, extrapolating
the dynamical behavior to broader ranges of parameters allow
us to design systems with higher efficiencies.
In our system, energy is provided by a gliding activation
front (AF) that activates an otherwise passive medium; like
a combine harvester machine. The progression of the AF
leaves behind a trail of activated material which evolves un-
der its active stress and induces a stress gradient that drives
the particle; see Fig.
1(b)
. Taking the experimental work as
our phenomenology, we formulate a model that allows us to
explore the phase diagram of the system. Our model involves
several parameters (see Table
I
), among which we choose the
following ones to study: (1) a measure of activity, (2) the
viscoelastic timescale, (3) the velocity of the activation front,
and (4) the coupling to the particle. Our theory provides in-
sight into the mechanism of transport within the active system
and also reveals how an upper bound on transport velocity
emerges within a contractile active system. Note that in our
setting, the particle is not a probe that extracts energy from
the active medium but instead interacts with both active and
passive components and “measures” the amount of extracted
energy expended by the active medium to transport the par-
ticle against the drag force. This is distinct from the viscous
023195-2
THEORETICAL LIMITS OF ENERGY EXTRACTION IN ...
PHYSICAL REVIEW RESEARCH
6
, 023195 (2024)
dissipation of the material that occurs even in the absence of
the particle. We show that the viscoelastic timescale which
characterizes the propagation velocity of the density wave is a
key determinant of the feasibility and the efficiency of trans-
port. In the context of active polymer networks, this timescale
is determined by the rate of cross-linking filaments and the
consequent increase of the rigidity over time.
A key difference between our setting and conventional
active systems is that here, the active and passive stresses can
cancel one another, where the system ceases to flow. There-
fore, even though active energy is required to be constantly
injected
locally
for maintaining the stress (force dipoles) in
the contracted state, a persistent flow requires nonzero diver-
gence of the stress, which occurs when fresh active material is
introduced to the system.
II. MODEL
The active system we adopt to study comprises a mixture
of active and passive components, the total density of which is
initially, and remains, uniform in space:
=
φ
a
+
φ
p
, where
φ
a
,
p
denote the corresponding densities of active and passive
components. To put our model in context we consider a mix-
ture where the passive component is a solution of free-floating
filaments and cross-linkers, that would assemble into an active
network (gelation) upon “activation” which entails recruit-
ing active cross-linkers [
45
]. Active networks of filaments
can be realized naturally in biological systems, or artificially
(e.g., using iLID technology [
14
,
44
,
46
]), and are driven out
of equilibrium by active cross-linkers that generate stress by
walking along the filaments. The magnitude of the active
stress is controlled by a parameter that depends on the density
of active cross-linkers [
47
]. We find it necessary to mention
that while we use active networks as a synthetic system, to
which our model is applicable, some aspects of the model
are independent of the details of such systems. Therefore we
define phenomenological parameters that would have various
microscopic interpretations depending on the system.
A. Active viscoelasticity
In our model, activation is defined as a process that
endows the passive component with the capacity of form-
ing an active gel of initial density
φ
i
. The density in
the case of polar filaments represents the density of fila-
ments’ plus (or minus) ends in the cross-linked network (see
Ref. [
47
]). While activation is instantaneous, gelation takes
place over a characteristic timescale of the system, suggesting
viscoelastic behavior [
45
,
48
–
50
]. The specific model of vis-
coelasticity to be adopted depends on the timescales in the
phenomenology of interest. There exist two minimal mod-
els of viscoelasticity, that can be represented by means of
a spring and a dashpot that are arranged either in series
or parallel configurations. The series configuration, called
Maxwell’s model describes high-frequency elastic followed
by low-frequency viscous behavior, where as the parallel
configuration, called Kelvin-Voigt model, represents the op-
posite behavior. The two models (Maxwell and Kelvin-Voigt)
both describe dynamics that transition between elasticity-
and viscosity-dominated modes. While the response of the
preformed polymer networks is primarily described in terms
of Maxwell’s viscoelasticity [
45
], the response of actively
assembling (as opposed to preformed) networks transitions
from viscous to elastic in time and can be captured by
Kelvin-Voigt model [
51
]. The elastic-to-viscous transition of
Maxwell’s model in preformed networks, originates from
breaking the cross-links under stress, which leads to relaxing
internal stresses. On the other hand, a network that is undergo-
ing cross-linking behaves like a fluid before the cross-linking
reaches the limit of stress percolation where the network
exhibits elastic behavior. As such, given that the activation
front in our system ignites the cross-linking in a previously
free-floating mixture of filaments, the short-term behavior is
dominated by viscosity and gradually transitions to elasticity;
hence the Kelvin-Voigt model.
In addition to the density field
φ
a
, the active component is
characterized by a displacement field
u
a
(
r
,
t
). Using Euler’s
notation for derivatives: D
t
f
≡
D
f
/
D
t
is the material time
derivative, and
∂
q
f
≡
∂
f
/∂
q
is the partial derivative with
respect to variable
q
, the velocity and strain fields equal
v
a
=
∂
t
u
a
, and
E
=
∇
u
a
. Together with compressible continuity
equation
∂
t
φ
a
+
∇
·
(
φ
a
v
a
)
=
0, momentum conservation and
the constitutive equation read:
D
t
(
φ
a
v
a
)
=
∇
·
+
f
,
(1a)
=
η
(
τ
−
1
K
+
D
t
)
E
=
κ
E
+
η
D
t
E
.
(1b)
Here
and
E
represent stress and strain fields;
η
and
τ
K
are the viscosity and the Kelvin timescale, and
κ
=
η/τ
K
denotes the rigidity of the elastic component. The body
force density
f
reads
f
=−
∇
·
int
−
γ
0
φ
a
φ
p
(
v
a
−
v
p
). The
second term is the drag force between passive and active
components that respects
{
a
}↔{
p
}
exchange symmetries.
Using the continuity equation
φ
a
v
a
+
φ
p
v
p
=
0, we can show
that the drag reduces to
−
γ
0
φ
a
v
a
≡−
γφ
a
v
a
(see Ap-
pendix
A
). The internal stress is the sum of contractile (active)
and an extensile (passive) components. The latter originates
from the collisions of the filaments, which guarantees stability
against the collapse of the network under contractile stress.
Using a virial-like expansion in terms of the density field we
get
int
=
σ
int
(
φ
a
)
I
=−
α
1
φ
a
I
+
α
2
φ
2
a
I
+
O
(
φ
3
a
)
=−
α
1
φ
a
(
1
−
φ
a
/φ
0
)
+
O
(
φ
3
a
)
,
(2)
where
α
1
,α
2
>
0 are constants, and
I
is the identity matrix.
To second order in the density field, the stress can be written
as
σ
int
=−
α
1
φ
a
(1
−
φ
a
/φ
0
), where
φ
0
=
α
1
/α
2
, is the den-
sity at which the active and passive densities cancel out one
another; but is in general distinct from equilibrium density for
finite Kelvin timescales
τ
K
; see below. Since we are interested
in one-dimensional transport, the tensors are reduced to scalar
fields
σ
and
ε
which represent
and
E
, respectively.
The signature of activity can be seen by noting that the
contractile stress and density form a positive
feedback loop
against the conventional equilibration in passive systems—
except for gravitational systems. The positive feedback holds
for
φ
a
φ
0
/
2, above which the effect of passive stress dom-
inates the contractile term. At equilibrium, the internal stress
023195-3
SHAHRIAR SHADKHOO AND MATT THOMSON
PHYSICAL REVIEW RESEARCH
6
, 023195 (2024)
is given by
σ
int
(
φ
e
)
=
κ
ε
e
. In Appendix
A
, we show that a
dimensionless number,
ν
=
α
2
1
τ
K
/ηα
2
=
α
2
1
/κα
2
, controls
the equilibrium density which varies between
φ
i
<φ
e
φ
0
for
α
1
α
2
φ
i
(the condition required for contraction in
the first place). In the limit of minuscule elasticity
κ
→
0, the asymptotic solution reads:
φ
e
=
φ
0
−
(
φ
0
−
φ
i
)
ν
−
1
+
O
(
ν
−
2
). In the opposite limit
ν
→
0, we get:
φ
e
=
φ
i
+
(
φ
2
i
/φ
0
)(1
−
φ
i
/φ
0
)
ν
+
O
(
ν
2
); also see Fig.
5
in the Ap-
pendix for numerical solutions.
Finally, the dissipation of energy occurs through the vis-
cosity of the active component and the drag against the
background medium. While the former is a characteristic of
the active material, the latter is medium-specific. Thus we
deliberately neglect the drag force to focus on the material
properties. In Sec.
V
below and in Appendix
C2
, we discuss
the effect of nonzero drag.
B. Field-particle interaction
The total force experienced by a particle embedded in the
solutions depends on passive and active pressures. To linear
order, the total pressure reads
P
=
P
a
+
P
p
, where
P
a
,
p
=
g
a
,
p
φ
a
,
p
, and
g
a
,
p
are the coupling constants. The driving force
equals:
F
=
∇
P
, where
is the particle’s linear size. Given
that
φ
a
+
φ
p
=
is constant, we write the pressure gradient
as
∇
P
=
(
g
a
−
g
p
)
∇
φ
a
=
g
0
∇
φ
a
, where
g
0
=
g
a
−
g
p
.Note
that
g
0
and
α
1
are of the same dimensions. In the experi-
ments of light-induced activity in the cross-linked networks,
high-density localized asters that are isolated from the passive
medium, can play the role of a rigid particle. With this analogy
in mind, we choose
g
0
>
0; meaning that the particle is pulled
by the newly activated regions. The case of
g
0
<
0 is briefly
discussed below and in Appendix
B
.
The motion of the particle is resisted by a drag force, which
can be shown to be equal to
−
V
p
, with
the effective
drag coefficient (see Appendix
B
). In the overdamped limit
the equation of motion reads:
V
p
=
g
0
∇
φ
a
. To focus on
the properties of the active system regardless of its specific
coupling to the particle, here we assume that the coupling to
the particle is negligible compared with the field’s internal
interactions:
g
0
α
1
; namely, we neglect the effect of the
particle on the field. The potential field to which the parti-
cle is exposed is obtained by
U
=−
g
0
∫
x
d
x
∇
φ
a
=−
g
0
φ
a
,
which is not influenced by the particle—similar to the system
studied in Ref. [
52
].
An interesting point is that the feasibility of transport de-
pends on the initial condition. If the
initial
positions of the
AF and the particle coincide, the transport is only possible
for
g
0
>
0. If the particle is introduced to an already activated
region (with density gradient), the transport can
also
occur for
g
0
<
0. In the latter case, the particle is trapped at the position
of the moving density domain wall which carries a negative
density gradient [see Fig.
2(a)
for the definition of the domain
wall]. The two scenarios introduce two potential wells that
could trap and drag the particle along, although with opposite
curvatures of the density profile. Here we are focusing on the
case of
g
0
>
0, which results in the particle being trapped at
the boundary of the system.
Throughout the paper we set
α
2
(as opposed to
φ
0
)
equal to unity; thus
φ
0
=
α
1
. The reason behind this choice
is that in real systems, the passive (excluded volume)
FIG. 2. (a) A schematic of the transport of a particle, as the
activation front [dark blue line in panels (b) and (c)] progresses in
time and activates the passive region. The newly activated matter
transitions from state 1 characterized by (
φ
1
,
̇
ε
1
), to the equilibrium
state 2 with (
φ
2
=
φ
e
,
̇
ε
2
=
0). The two states are separated by the
DW moving with velocity
V
DW
. The left boundary moves with
v
b
.
The density profiles and the position of the particle are shown for
two sets of parameters: (b) large enough coupling for the activation
front to drag the particle, hence steady transport; and (c) for small
coupling the particle is detached and the transport fails.
pressure, characterized by
α
2
, is independent of the activity.
Furthermore, even though the field’s viscosity
η
, the particle’s
drag coefficient
, and the initial density
φ
i
are kept around in
our calculations, the numerical results are obtained by setting
them equal to unity.
III. ACTIVE TRANSPORT
The hypothesis we follow is that, by continuously acti-
vating the material in a specific direction with velocity
V
AF
,
contractile stress provides the power required to move the
particle. Given the uniformity of the passive region,
V
AF
is
proportional to and can be interpreted as the
rate of the injec-
tion of active energy
into the system. We numerically explore
the phase diagrams of transport spanned by (
α
1
,
V
AF
,
g
0
);
see Fig.
3
. We elaborate on the numerical methods in Ap-
pendix
E
. In specific limits we also derive analytical solutions;
Sec.
III A 2
.
To gain intuition into the dynamics we use a discrete pic-
ture. Moving the activation front in
+
ˆ
x
direction, the activated
023195-4