# Geometric analysis of gaits and optimal control for three-link kinematic
swimmers

Oren Wiezel, Suresh Ramasamy, Nathan Justus, Yizhar Or and Ross Hatton The
work of OW and YO has been supported in part by Israel Ministry of Science and
Technology, under contracts no. 3-15622 and 3-17383The work of SR, NJ and RH
has been supported in part by the National Science Foundation awards 1653220
and 1935324OW and YO are with the Faculty of Mechanical Engineering at the
Technion - Israel Institute of Technology, Haifa 3200003, Israel (e-mail:
<EMAIL_ADDRESS>izi@technion.ac.il). SR, NJ and RH are Oregon
State University, Corvallis, OR 97333 USA. (e-mail<EMAIL_ADDRESS>,
<EMAIL_ADDRESS>Ross.Hatton@oregonstate.edu).

###### Abstract

Many robotic systems locomote using gaits – periodic changes of internal
shape, whose mechanical interaction with the robot’s environment generate
characteristic net displacements. Prominent examples with two shape variables
are the low Reynolds number 3-link “Purcell swimmer” with inputs of 2 joint
angles and the “ideal fluid” swimmer. Gait analysis of these systems allows
for intelligent decisions to be made about the swimmer’s locomotive
properties, increasing the potential for robotic autonomy. In this work, we
present comparative analysis of gait optimization using two different methods.
The first method is variational approach of “Pontryagin’s maximum principle”
(PMP ) from optimal control theory. We apply PMP for several variants of
3-link swimmers, with and without incorporation of bounds on joint angles. The
second method is differential-geometric analysis of the gaits based on
curvature (total Lie bracket) of the local connection for 3-link swimmers.
Using optimized body-motion coordinates, contour plots of the curvature in
shape space give visualization that enables identifying distance-optimal gaits
as zero level sets. Combining and comparing results of the two methods enables
better understanding of changes in existence, shape and topology of distance-
optimal gait trajectories, depending on the swimmers’ parameters.

##  I Introduction

Robotic swimmers are a promising avenue of research. Both small microswimmers
and large scale swimmers have many promising possible applications. Advances
in technology for manufacturing nano- and micro-systems have brought renewed
interest in simplified models of microswimmers and production of microscopic
robotic devices that could be applied in the medical field [1, 2, 3]. Such
micro-robots would be able to provide targeted drug delivery, tumor detection,
assisted sperm motility, and even perform minimally invasive surgical
procedures. At the other end of size scale, large, fast moving swimming robots
can be useful in search and rescue missions, maintenance operations within
pipe systems of complex infrastructures, and surveillance or protection in
marine environments [4].

The flow around microswimmers is governed by Stokes equations, which arise
from the Navier-Stokes equations at the limit of zero Reynolds number $Re\to
0$ [5, 6]. At very small Reynolds numbers, $Re\ll 1$, inertial forces become
negligible and the viscous forces are dominant. This leads to a linear
relationship between the body velocities and the internal shape velocities
[7]. These unique attributes call for drastically different swimming
strategies than the ones used in the familiar motion of large organisms.

Purcell suggested that the “simplest animal” that could swim is the three-link
swimmer [8], comprised of three thin rigid links connected by two rotary
joints (Fig. 1(a)). By alternately rotating the joint angles, this swimmer
would be able to propel itself in a low-$Re$ environment. The series of shape
changes, or “gait” proposed by Purcell appears as a square in the plane of
shape variables (the joint angles) as shown in Fig. 1(b). This gait results in
motion along the mean orientation of the central link. Becker et al. [9]
formulated the dynamics of Purcell’s swimmer using slender-body theory [10,
11] and found that the direction of net displacement of the swimmer depends on
the angular amplitude of the strokes. For small amplitude the swimmer will
move in one direction, but for larger amplitudes the translation will be in
the opposite direction (Fig. 2). This also implies the existence of an
amplitude of the gait that will result in maximal translation (marked with a
Purple X in Fig. 2). For very large amplitudes there is a second maximum-
displacement gait that will result in swimming in the opposite direction
(marked with a red X in Fig. 2).111We note that maximizing displacement per
cycle is not a true “optimal control” problem for the systems considered in
this paper, because per-cycle displacement does not account for the different
amounts of time or energy consumed by “large” or “small” gaits. Displacement-
maximizing gaits do, however, have a distinct mathematical structure and exist
in close proximity to speed and efficiency maximizing gaits [12], and so are
therefore an interesting topic of study. Importantly, in robotic realizations
of such swimmers, cf. [13], one commonly has to account for practical bounds
on joint angles which may limit the stroke amplitude.

(a)

(b)

Figure 1: (a) Purcell’s three-link swimmer. (b) The square, circular and
maximal displacement gaits. ++ Figure 2: Net-displacement over a period in the
$x$ direction as a function of the gait amplitude $\varepsilon$. Circular gait
in blue full line and square gait in dashed’ green line. X marks denote
optimal amplitudes for maximal (purple) and minimal (red) displacement for
circular gait.

In contrast to small microswimmers whose motion is governed by viscous drag
forces, the motion of large swimmers, such as fish or eels, is governed mainly
by inertial forces [14, 15]. A common simplified mathematical model for such
swimmers is that of “perfect fluid” that assumes potential flow in a inviscid
fluid and allows for an added mass effect [16]. Notably, using this model for
a swimmer in an unbounded fluid domain that dictates symmetry and conservation
laws leads to first order equations of motion that give a linear relation
between the body velocities and internal shape velocities. Therefore, motion
analysis of an inertial swimmer under this “perfect fluid” model can be done
using the same mathematical methods as the in the case of a microswimmer,
despite the differences in the physics of the mechanism creating the motion
[14, 17]. Previous works have examined the reduction of inertia-driven
“perfect fluid” dynamics into a first order system [14, 18, 19].

The common basis for the small-scale microswimmers and large-scale “perfect
fluid” swimmers is that their motion can be described by a “principal
kinematic” structure [19], where the body velocities are linearly related to
the shape velocities. Based on the Lie-group structure of rigid-body motion,
this locomotion principle is not unique to swimming and is shared by other
models of robotic locomotion induced by internal shape changes, such as the
terrestrial snake [20] and some wheeled robots [21, 22].

Optimization of design and actuation of low-$Re$ swimmers has been a subject
of several recent studies. One approach to gait optimization was taken in
[23], where geometric mechanics techniques based on Lie brackets were used to
provide an approximation of the net displacement for locomoting systems. These
approximations are then used to search for the optimal gait. Historically, the
approximation error has limited this analysis to small shape changes. Recent
works by Hatton and Choset [24] show that the choice of body coordinates can
significantly reduce the approximation error. This allowed them to find gaits
maximizing displacement per cycle of the three-link swimmer by examining zero-
contour lines of the total Lie bracket, as well as gaits for maximal energetic
efficiency which also maximize displacement per time at a given effort level
[25].

Another important approach for finding optimal gait trajectories is using
variational methods of optimal control, namely Pontryagin’s maximum principle
(PMP ) [26], [27], [28, 29]. This approach has been vastly applied in several
practical fields such as aerospace applications [30, 31, 32] and computer
vision [33], as well as biology [34, 35] and quantum dynamics [36, 37]. For
Purcell’s three-link microswimmer model, the PMP method has been applied in
[38], and analytically reproduced the maximum displacement-per-cycle gait
found in previous works [39, 19] using different approaches. The gait in [38]
was found using a polar representation of the shape variables and so is not
suitable for closed curves with a general shape. Another limitation of the PMP
formulation in [38] is that it did not allow for bounds on the shape
variables, which commonly exist in any practical application.

Extension of the maximum principle to account for inequality constraints on
state variables have been presented in several works, e.g. [30, 40, 41, 42,
43]. This extension has also been recently incorporated in optimal control for
specific applications such as autonomous transportation and vehicles [44, 45,
46], energy management systems [47, 48], and aerospace systems [49, 50].

A different variant of optimal control, which has significance for
computational aspects, is the discrete-time PMP [51, 52, 53]. In particular,
several works have focused on systems evolving on Lie groups, and development
of discrete integrators that preserve group invariants and thus avoid
cumulative integration error [54, 55, 56]. This method can incorporate bounds
on input and state variables due to its discrete nature. For Purcell’s three-
link microswimmer model, the recent work [57] considered the iso-holonomic
problem of finding gaits that minimize energy expenditure under given net
displacement per cycle under given state bounds, by employing discrete-time
PMP on a Lie group.

The goal of this work is to revisit the kinematic models of three-link
swimmers and explore optimal gait trajectories of joint angle inputs for
achieving maximal net displacement per cycle, with possible bounds on joint
angles. We make a comparative analysis of two different methods. The first
method is based on necessary conditions of variational calculus, using PMP
within the framework of optimal control. It requires solving a two-point
boundary-value problem which reduces to a simple “shooting” solution of a
scalar variable. The method captures locally-optimal smooth gait trajectories
for maximal positive displacement. However, for some swimmer’s model parameter
it fails to find a smooth trajectory of large-amplitude gait for maximizing
negative displacement, whereas a non-smooth locally-optimal gait can be found
when state bounds are incorporated.

The second method is based on differential geometry and on using an
approximate solution (formally a truncated Baker-Campbell-Hausdorff series) to
the net-displacement-per-cycle equations. This approximation allows the net
displacement over a gait cycle to be computed as the volume integral of a
scalar height function over the region enclosed by the gait trajectory.
Through a principled choice of coordinates for the system (in particular,
selecting a body frame at a nonlinearly weighted average of the link positions
instead of being attached to a link), we have previously shown that this
approximation can be made usefully accurate over large-amplitude motions of
the joints [23].

In particular, a very accurate approximation of gait trajectories that
maximize displacement per cycle is obtained by the loops of zero-level curves
of the scalar height function. Upon varying swimmer’s parameters, changes in
height function landscape and in topology of its zero-level curves give a
graphical explanation to the change in shape and nature of such gait
trajectories. In particular, junction points in the zero-level curves preclude
existence of locally-optimal trajectories, and state bounds become necessary
in order to obtain non-smooth optimal gaits that avoid the junction points.
Overall, the two methods reconcile together to combine a rigorous variational
formulation with graphical visualization of the changes in shape and nature of
optimal gaits.

The structure of the paper is as follows. In the next section, we present two
models of swimmers: Purcell’s swimmer and the “perfect fluid” swimmer. In
section III we offer the formulation of the optimal control problem using PMP,
and show our solutions for the models with and without joint angle bounds, as
well as cases where this formulation fails to find smooth optimal gaits. Next,
in section IV we review the geometric analysis of gaits introduced in [23]. In
section V we compare the results of the two methods. The closing section
contains discussion and conclusion.

##  II Mathematical models of three-link kinematic swimmers

We begin by formulating the dynamics of two models of 3-link swimmers. The
swimmer model consists of three thin rigid links with lengths
$l_{0},l_{1},l_{2}$, where $l_{1}=l_{2}$ . The links are connected by two
rotary joints whose angles are denoted by $\phi_{1}$ and $\phi_{2}$ (see Fig.
1(a)). The shape of the swimmer is denoted by the vector
$\boldsymbol{\upphi}=(\phi_{1},\phi_{2})^{T}$. It is assumed that the
swimmer’s motion is confined to a plane. The planar position of the middle
link’s center is $(x,y)$ and its orientation angle is $\theta$. The velocity
of the central link in an inertial frame is denoted by
$\dot{\mathbf{q}}_{b}=(\dot{x},\dot{y},\dot{\theta})$. The velocity of the
$i$th link is described by the linear velocity of its center and the links
angular velocity $\omega_{i}$, which are augmented in the vector
$\mathbf{v}_{i}=(\dot{x}_{i},\dot{y}_{i},\omega_{i})\in\mathbb{R}^{3}$. We
denote the body velocities by
$\mathring{\mathbf{q}}=(v_{x},v_{y},\dot{\theta})$. These are the velocities
of the central link expressed in a reference frame attached to the central
link. We show below that for each of the models, the relation between the
shape velocities and the body velocities can be written as:

$\accentset{\scriptstyle\circ}{\mathbf{q}}=\mathbf{A}(\boldsymbol{\upphi})\dot{\boldsymbol{\upphi}}$
(1)

also known as the kinematic reconstruction equation [58]. In order to obtain
the velocities in an inertial frame we multiply by a rotation matrix:

$\dot{\mathbf{q}}_{b}=\mathbf{R}(\theta)\accentset{\scriptstyle\circ}{\mathbf{q}}=\mathbf{R}(\theta)\mathbf{A}(\boldsymbol{\upphi})\dot{\boldsymbol{\upphi}}=\mathbf{G}(\theta,\boldsymbol{\upphi})\dot{\boldsymbol{\upphi}}$
(2)

While in this work we only consider planar three-link models, this form can
also be generalized to multi-link models and to spatial swimmers using three
rotation angles [59].

In the rest of this section, we give a short review on how to write the
relation (1) for both swimmer models.

###  II-A Purcell’s swimmer

The first swimmer model we consider is “Purcell’s swimmer”, in which the
system is submerged in an unbounded fluid domain whose motion is governed by
Stokes equations [5]. Here, we briefly summarize how to find the local
connection for Purcell’s swimmer. More in-depth derivations can be found using
inertial-frame and link-frame methods described in [25] and [60].

First, the velocity of each link through the fluid is written using Jacobian
matrices $\mathbf{J}_{i}$ relating overall swimmer body velocity and joint
velocities to individual link velocities:

$\mathbf{v}_{i}=\mathbf{J}_{i}\begin{bmatrix}\mathring{\mathbf{q}}\\\
\dot{\boldsymbol{\upphi}}\end{bmatrix}$ (3)

Resistive force theory is used to calculate the force on each link as being
linearly proportional to that link’s velocity $\mathbf{v}_{i}$ through a
resistance matrix $\mathbf{D}_{i}$ :

$\mathbf{F}_{i}=-\mathbf{D}_{i}\mathbf{v}_{i}$ (4)

Then, forces on each link can be mapped back into the body frame of the
swimmer using a dual adjoint-inverse mapping and summed together to obtain net
load on the body $\mathbf{F}_{b}$:

$\mathbf{F}_{b}=\sum_{i=0}^{2}{\mathbf{J}_{i}^{T}}\mathbf{F}_{i}=\left(\sum_{i=0}^{2}-{\mathbf{J}_{i}^{T}}\mathbf{D}_{i}\mathbf{J}_{i}\right)\begin{bmatrix}\mathring{\mathbf{q}}\\\
\dot{\boldsymbol{\upphi}}\end{bmatrix}$ (5)

At low Reynolds numbers, net forces on the swimmer are zero because it is at
quasistatic equilibrium. This allows us to write equation (5) as nonholonomic
constraints on velocities:

$\mathbf{F}_{b}=\mathbf{C}_{\mathbf{q}}{(\boldsymbol{\upphi})}\mathring{\mathbf{q}}+\mathbf{C}_{\boldsymbol{\upphi}}{(\boldsymbol{\upphi})}\dot{\boldsymbol{\upphi}}=0$
(6)

This allows us to solve for the viscous local connection:

$\mathring{\mathbf{q}}=-\mathbf{C}_{\mathbf{q}}^{-1}\mathbf{C}_{\boldsymbol{\upphi}}\dot{\boldsymbol{\upphi}}=\mathbf{A}(\boldsymbol{\upphi})\dot{\boldsymbol{\upphi}}$
(7)

###  II-B Perfect fluid swimmer

The second swimmer we consider is the “perfect fluid” swimmer model (Fig. 3),
whose dynamics are constructed by assuming inviscid and irrotational potential
flow. In [18], analysis for sinusoidal input was done and motion experiments
were performed with a controlled robotic swimmer in order to validate this
model. Each link is assumed to be an ellipse with principal radii of $a_{i}$,
$b_{i}$ and density $\rho$, which has mass $m_{i}$ and moment of inertia
$I_{i}$. The swimmer is submerged in an unbounded domain of ideal fluid with
equal density $\rho$.

The kinetic energy of the system can be written as a sum of the kinetic
energies of the individual links:

$T=\frac{1}{2}\sum_{i}\mathbf{v}_{i}^{T}\mathbf{M}_{i}\mathbf{v}_{i}$ (8)

where

$\mathbf{M}_{i}=\left[\begin{array}[]{c c c}m_{i}&0&0\\\ 0&m_{i}&0\\\
0&0&I_{i}\end{array}\right]+\pi\rho\left[\begin{array}[]{c c
c}b_{i}^{2}&0&0\\\ 0&a_{i}^{2}&0\\\
0&0&\frac{1}{8}(a_{i}^{2}-b_{i}^{2})\end{array}\right]$

The first term in $\mathbf{M}_{i}$ represents the inertia of the link while
the second term represents the added mass effect due to accelerating the
displaced fluid around the swimmer [14].

Using the link Jacobian to relate body-frame velocity and joint velocity to
individual link velocity as in (3), we can rewrite the kinetic energy in terms
of these two quantities:

$T=\frac{1}{2}\left[\begin{array}[]{c }\mathring{\mathbf{q}}\\\
\dot{\boldsymbol{\upphi}}\end{array}\right]^{T}\left[\begin{array}[]{c c
}\mathbf{M}_{bb}&\mathbf{M}_{bs}\\\
\mathbf{M}_{bs}^{T}&\mathbf{M}_{ss}\end{array}\right]\left[\begin{array}[]{c
}\mathring{\mathbf{q}}\\\ \dot{\boldsymbol{\upphi}}\end{array}\right]$ (9)

where all blocks $\mathbf{M}_{bb},\mathbf{M}_{bs},\mathbf{M}_{ss}$ are
functions of the shape variables $\boldsymbol{\upphi}$ only.

Figure 3: The three-link “perfect fluid” swimmer model [18].

Using generalized coordinates
$\mathbf{q}=(\mathbf{q}_{b},\boldsymbol{\upphi})$, the system’s dynamic
equations of motion can be written using Lagrange’s formulation [61] as

$\frac{d}{dt}\left(\frac{\partial
T}{\partial\dot{\mathbf{q}}}\right)-\frac{\partial
T}{\partial\mathbf{q}}=\mathbf{F}_{q},$ (10)

where $\mathbf{F}_{q}$ is a vector of generalized forces. In matrix form, this
equation reads as

$\mathbf{M}(\mathbf{q})\ddot{\mathbf{q}}+\mathbf{b}(\mathbf{q},\dot{\mathbf{q}})=\mathbf{F}_{q}.$
(11)

The system’s inertia matrix $\mathbf{M}(\mathbf{q})$ in (11) has tight
relation to the block matrices in (9) due to invariance of the dynamics with
respect to rigid body transformations, which induces conservation of
generalized momentum variables, see [18, 19, 22]. Assuming that the swimmer
starts from rest, $\mathbf{q}(0)=0$, this momentum conservation can be reduced
to a first-order relation between body velocity and shape changes as (see
[14]):

$\mathbf{M}_{bb}(\boldsymbol{\upphi})\mathring{\mathbf{q}}+\mathbf{M}_{bs}(\boldsymbol{\upphi})\dot{\boldsymbol{\upphi}}=0$
(12)

The body velocities can now be written in the form of equation (1) with
$\mathbf{A}(\boldsymbol{\upphi})=-\mathbf{M}_{bb}^{-1}\mathbf{M}_{bs}$.

##  III Optimal control using PMP

Applying periodic input of shape changes (gaits) to the swimmer models
presented above, results in net motion per cycle, which can be obtained by
integrating (2) over a period. Importantly, the kinematic relation (2) is
time-invariant, and corresponds to a geometric relation

$d\mathbf{q}_{b}=\mathbf{G}(\theta,\boldsymbol{\upphi})d\boldsymbol{\upphi}.$
(13)

This implies that for any geometric shape of gait trajectory
$\boldsymbol{\upphi}(s)$, the resulting net body motion $\Delta\mathbf{q}_{b}$
generated by this gait is independent of rescaling the time-rate $s(t)$ of the
input [62]. One may seek for optimal gaits that maximize per-cycle
displacement [19]. In order to find gaits that maximize per-cycle displacement
of the swimmer models presented, we apply methods of variational optimal
control, which lead to the PMP [27]. We first present the formulation of a
general solution via PMP and the formulation for our problems with dynamics
(1), with and without state bounds. Next, we present the optimal gaits found
using this method for both of the swimmer models and the influence of changes
in problem’s bounds and model parameters on the optimal solution.

###  III-A Formulation of OCP

We start with a short description of the optimal control problem and PMP
solution for a model with the dynamics in (1).

For an optimal control problem

$(OCP)\begin{cases}\max J=\int_{0}^{t_{f}}g(\mathbf{z},\mathbf{u})\ s.t.\\\
\dot{\mathbf{z}}=\mathbf{f}(\mathbf{z},\mathbf{u})&\forall t\in[0,t_{f}],\\\
\mathbf{u}\in\mathcal{U}&\forall t\in[0,t_{f}],\\\
\mathbf{z}(0)=\mathbf{z}_{0},&\end{cases}$ (14)

where $\mathbf{z}\in\mathcal{R}^{n}$ is the state of the system and
$\mathbf{u}$ is the control input, we define a costate vector
$\boldsymbol{\lambda}\in\mathcal{R}^{n}$ and the Hamiltonian will be

$H(\mathbf{z},\mathbf{u},\boldsymbol{\lambda})=\boldsymbol{\lambda}^{T}\mathbf{f}+g$
(15)

Pontryagin’s Maximum Principle (PMP) states that a necessary condition for
obtaining the optimal control input trajectory $\mathbf{u}^{*}(t)$ with
associated state trajectory $\mathbf{z}^{*}(t)$ is that it maximizes the
Hamiltonian:

$\mathbf{u}^{*}(\mathbf{z},\boldsymbol{\lambda})=\operatorname*{arg\,max}_{\mathbf{u}\in\mathcal{U}}H(\mathbf{z},\mathbf{u},\boldsymbol{\lambda}).$
(16)

Thus, an optimal input can be found by solving

$\mathbf{H}_{\mathbf{u}}=\frac{\partial H}{\partial\mathbf{u}}=\vec{0}$ (17)

.

The Maximum Principle requires solving an ODE with two-point boundary
conditions.

#### Singular arcs

In problems where the Hamiltonian is linear in the input, the optimal input
cannot be found from (17). In many cases, this implies a “bang-bang” solution
where the control switches between the upper and lower bounds. The switching
is determined by the sign of $\mathbf{H}_{\mathbf{u}}$. In some cases, when
$\mathbf{H}_{\mathbf{u}}$ can vanish for a finite time, the solution follows a
singular arc and can be determined by the time derivatives:

$\frac{d^{k}}{dt^{k}}\mathbf{H}_{\mathbf{u}}=\vec{0}$ (18)

After an even number of derivatives $k$, the input $\mathbf{u}$ appears and
can be extracted.

#### Bounded state

When formulating the maximum principle with inequality bounds on the state
variables, we follow the direct adjoining approach [42, 43]. For an inequality
bound in the form:

$w(\mathbf{z},t)\leq 0$ (19)

with time derivatives

$w^{k}(\mathbf{z},t)=\frac{d^{k}w}{dt^{k}}$ (20)

The Hamiltonian is defined as

$\tilde{H}(\mathbf{z},\mathbf{u},\boldsymbol{\lambda},\nu)=\boldsymbol{\lambda}^{T}\mathbf{f}+g+\nu
w$ (21)

where $\nu$ is an additional multiplier, and

$\nu w^{*}(\mathbf{z},t)=0$ (22)

so that $\nu=0$ whenever the bound is inactive ($w(\mathbf{z},t)<0$), while
$w(\mathbf{z},t)=0$ and $\nu>0$ whenever the bound is active.

A time interval $[\tau_{1},\tau_{2}]$ is called a bounded interval if
$w(\mathbf{z},t)=0$ for $t\in[\tau_{1},\tau_{2}]$ (the bound is active).
$\tau_{1}$ and $\tau_{2}$ are the entry time and exit time, respectively.

Whenever the bound is inactive ($w(\mathbf{z},t)<0$), the solution is found
through $\mathbf{H}_{\mathbf{u}}=\vec{0}$ (or its time derivatives (18) in the
case of a singular arc). When the bound is active, $w(\mathbf{z},t)=0$ and its
time derivatives determine the control input $u^{*}$. Generally, the costate
variables as well as the Hamiltonian do not necessarily change continuously at
entry time and exit time and additional necessary conditions can be written
that determine the discontinuity [42]. But, if the control input is
discontinuous over the entry or exit time, $u^{*}(\tau^{-})\neq
u^{*}(\tau^{+})$, then the costate variables and Hamiltonian must be
continuous, $\boldsymbol{\lambda}(\tau^{-})=\boldsymbol{\lambda}(\tau^{+})$,
$\tilde{H}^{*}(\tau^{-})=\tilde{H}^{*}(\tau^{+})$ [43] ( $\tau^{-}$ and
$\tau^{+}$ denote the left and right limits, respectively).

###  III-B OCP for three-link swimmer

We now formulate the displacement-per-cycle maximizing control problem for a
three-link kinematic swimmer and find the solution. Whereas the formulation in
[38] used polar representation of the shape variables, here we write a more
general formulation in order to allow for non-polar gaits. For the three-link
kinematic swimmer, we define the state as
$\mathbf{z}=[\phi_{1},\phi_{2},\theta]^{T}$, with the dynamics:

$\dot{\mathbf{z}}=\left[\begin{array}[]{c }\dot{\phi}_{1}\\\ \dot{\phi}_{2}\\\
\dot{\theta}\end{array}\right]=\left[\begin{array}[]{c }u_{1}\\\ u_{2}\\\
p(\phi_{1},\phi_{2})u_{1}+q(\phi_{1},\phi_{2})u_{2}\end{array}\right]$ (23)

with the input vector
$\mathbf{u}=[u_{1},u_{2}]^{T}=[\dot{\phi}_{1},\dot{\phi}_{2}]^{T}$. The cost
function is the net-displacement over a period in the $x$-direction:

$J=x(t_{f})=\int_{0}^{t_{f}}\dot{x}dt=\int_{0}^{t_{f}}g(\phi_{1},\phi_{2},\theta)u_{1}+h(\phi_{1},\phi_{2},\theta)u_{2}dt$

(24)

The Hamiltonian $H$ is given by:

$H(\mathbf{z},\mathbf{u},\boldsymbol{\lambda})=gu_{1}+hu_{2}+\lambda_{1}u_{1}+\lambda_{2}u_{2}+\lambda_{3}(pu_{1}+qu_{2})$
(25)

where $\boldsymbol{\lambda}$ is the vector of costate variables with dynamics:

$\dot{\boldsymbol{\lambda}}=-\frac{\partial
H}{\partial\mathbf{z}}=\leavevmode\resizebox{260.17464pt}{}{$\left[\begin{array}[]{c
}-g_{\phi_{1}}u_{1}-h_{\phi_{1}}u_{2}-\lambda_{3}(p_{\phi_{1}}u_{1}+q_{\phi_{1}}u_{2})\\\
-g_{\phi_{2}}u_{1}-h_{\phi_{2}}u_{2}-\lambda_{3}(p_{\phi_{2}}u_{1}+q_{\phi_{2}}u_{2})\\\
-g_{\theta}u_{1}-h_{\theta}u_{2}\end{array}\right]$}$ (26)

with subscripts denoting partial derivatives, for example:
$g_{\phi_{1}}=\frac{\partial g}{\partial\phi_{1}}$. The Hamiltonian is linear
in $u_{1},u_{2}$, and the control inputs do not appear in the derivative:

$\mathbf{H}_{\mathbf{u}}=\frac{\partial
H}{\partial\mathbf{u}}=\left[\begin{array}[]{c }g+\lambda_{1}+\lambda_{3}p\\\
h+\lambda_{2}+\lambda_{3}q\end{array}\right]$ (27)

Since we may have $H_{u}=0$ for a finite time period, the optimal control
following a singular arc may be found by requiring that the time derivatives
$\dot{\mathbf{H}}_{\mathbf{u}}$ and $\ddot{\mathbf{H}}_{\mathbf{u}}$ vanish as
well. Substituting (23) and (26) into the derivative of (27) we have the
first-order derivative as:

$\dot{\mathbf{H}}_{\mathbf{u}}=\leavevmode\resizebox{325.215pt}{}{$\left[\begin{array}[]{c
}u_{2}\\\
-u_{1}\end{array}\right]\overbrace{(g_{\phi_{2}}-h_{\phi_{1}}+g_{\theta}q-h_{\theta}p+\lambda_{3}p_{\phi_{1}}-\lambda_{3}q_{\phi_{1}})}^{\Psi(\mathbf{z},\boldsymbol{\lambda})}$}$
(28)

Dismissing the trivial solution of zero control input, the requirement of
$\ddot{\mathbf{H}}_{\mathbf{u}}=0$ reduces to the scalar equation
$\frac{d}{dt}\Psi(\mathbf{z},\boldsymbol{\lambda})=0$, which leads to an
equation in the form:

$A(\boldsymbol{\upphi},\theta,\lambda_{3})u_{1}+B(\boldsymbol{\upphi},\theta,\lambda_{3})u_{2}=0$
(29)

The last equations gives the angle of the local tangent to the gait’s curve
$\frac{d\phi_{2}}{d\phi_{1}}$. Due to the time invariance property of the
motion, eq. (13), one only needs this tangent in order to construct the
geometric shape of the optimal gait trajectory $\boldsymbol{\upphi}(s)$.
Therefore, we arbitrarily choose time-scaling $u_{1}^{2}+u_{2}^{2}=1$ and
write the controls as
$u_{1}=B/\sqrt{A^{2}+B^{2}},u_{2}=-A/\sqrt{A^{2}+B^{2}}$. The final time
$t_{f}$ is unspecified. Due to the symmetries of our system [13], we may
restrict our analysis to one quarter of the gait (from $\phi_{1}=\phi_{2}$ to
$\phi_{1}=-\phi_{2}$) and invoke symmetry. The total displacement will be four
times the displacement over the quarter gait. Because the rotation angle is
known only at the initial time $\theta(0)=0$ and $\theta(t_{f})$ is unknown,
using transversality conditions [27] we have $\lambda_{3}(t_{f})=0$. We also
have the boundary conditions on the shape variables $\phi_{1}(0)=\phi_{2}(0)$,
$\phi_{1}(t_{f})=-\phi_{2}(t_{f})$. The costate variables
$\lambda_{1},\lambda_{2}$ do not appear in the solution (29) and we do not
need to solve for them. Using the relation $\dot{\mathbf{H}}_{\mathbf{u}}=0$
we can find $\lambda_{3}(0)$ as a function of the state variables. We are left
with a system of four ODEs for the variables
$(\phi_{1},\phi_{2},\theta,\lambda_{3})$ where all but $\phi_{1}$ are known at
the initial time $t=0$. Using the shooting method we can find $\phi_{1}(0)$
that results in $\lambda_{3}(t_{f})=0$ with the relation
$\phi_{1}(t_{f})=-\phi_{2}(t_{f})$ defining the final time.

#### Solution with bounded joint angle

Assume a practical bound on the joint angles $\left|\phi_{i}\right|\leq b$.
Over a finite time interval with a non-zero control input, only the bound on
one joint angle can be active, while the other joint angle is varying. Using
the symmetry properties of the swimmer, we only consider one quadrant of the
gait where only one joint may reach the bound. Therefore, it is sufficient to
consider only a single, scalar state bound of the form (19). For our
demonstration, we will assume, without loss of generality, that the only bound
is $\phi_{2}\leq b$. We define

$w(\mathbf{z},t)=\phi_{2}-b\leq 0$ (30)

Assuming $\phi_{2}(0)<b$, the optimal gait starts on the singular arc and the
control is found from (29). On the bound $\phi_{2}=b$ we have
$\dot{w}=u_{2}=0$ which leads to $u_{1}=\pm 1$. The entry to the bound
$\tau_{1}$ is the time when the singular arc reaches the bound $\phi_{2}=b$.
We determine the exit time $\tau_{2}$ using the shooting method to satisfy the
end condition $\lambda_{2}(t_{f})=\lambda_{1}(t_{f})$. This means that we must
solve for $\lambda_{1},\lambda_{2}$ as well and ascertain the discontinuities
at the entry and exit times. As stated earlier, the discontinuity of the
control inputs at the entry and exit point implies continuity of the costate
variables and the Hamiltonian. The transversality conditions lead to
$\lambda_{1}(0)=-\lambda_{2}(0)$ and $\lambda_{1}(t_{f})=\lambda_{2}(t_{f})$.
We can find $\lambda_{1}(0)$ from $\mathbf{H}_{\mathbf{u}}(0)=\vec{0}$.

###  III-C Maximum-displacement gaits for three-link swimmers

We now present the solutions of the maximum displacement-per-cycle control
problem for the three-link swimmer models presented above, and the influence
of input bounds and swimmers’ parameters on existence and topology of optimal
solutions. The method presented in the previous section results in a set of
differential equations. These are solved using Matlab’s ode45 function. Event
functions detect crossing the bounds or reaching the final time at
$\phi_{1}=-\phi_{2}$. When a bound is reached, the simulation continues along
the bound until a exit time $\tau_{2}$. The initial value
$\phi_{1}(0)=\phi_{2}(0)$ and the exit time $\tau_{2}$ are found using fzero
function to satisfy $\lambda_{1}(t_{f})=\lambda_{2}(t_{f})$ and
$\lambda_{3}(t_{f})=0$

#### Purcell’s swimmer

For Purcell’s three-link swimmer model, a gait for maximal displacement was
presented in [38] using a polar representation of the shape variables. The
optimal gait (seen in purple in Fig. 4) is identical to that found numerically
in [39]. Restricting to circular gaits, one can see (Fig. 2) that there exists
a second optimal amplitude, that when followed in the same direction (counter-
clockwise) will give a net-displacement in the opposite direction with a
greater absolute value. Hence, it is reasonable to expect a second, general,
large amplitude, displacement-maximizing trajectory that would result in
translation in the opposite direction. Our attempts to find a second
displacement-maximizing gait for Purcell’s swimmer using PMP failed for the
unbounded problem. Nevertheless, when applying bounds of $b=3.2$ [rad] to the
problem, the gait shown in red in Fig. 4 is found. Gaits for some smaller
bounds are also shown in Fig. 4. Curiously, this method fails to find a
solution for larger bounds. In section IV we attempt to use geometric analysis
to explain this failure and why the unbounded problem does not have a
solution.

$X<0$$X>0$ Figure 4: Maximal displacement per cycle gait (purple) and second
optimal gait with bounds of $\phi_{i}\leq 3.6$ [rad] (red) and $\phi_{i}\leq
3.1$ [rad] (blue) for Purcell’s swimmer. The dashed lines are the
corresponding optimal circular gaits.

(a)

(b)

Figure 5: Maximal displacement gaits for perfect fluid swimmer with (a)
$\eta=1/3$ (red curve achieved only with bounded joints) and (b) $\eta=1/2$.

#### Perfect fluid swimmer

The maximum-displacement-per-cycle gait for a perfect fluid swimmer model with
elliptical links having radii ratio of $\alpha=a_{i}/b_{i}=0.2$ and links’
length ratio $\eta=1/3$ is presented in Fig. 5(a). The gait is qualitatively
similar to that found for Purcell’s swimmer and can also be found by the polar
formulation in [38]. Here too, we can expect to find a second, large-angle
gait that results in opposite displacement of the swimmer, $X<0$. However, for
a swimmer with this specific geometry, the unbounded OCP had no solution for
such gait. Once again, when imposing bounds on the joint angles, a gait with
$X<0$ is found for the perfect fluid swimmer as well. Fig. 5(a) shows the
second displacement-maximizing gait with joint angle bounds of $|\phi_{i}|\leq
3.2$. Unlike Purcell’s swimmer, for the perfect fluid model we found certain
parameter values that would yield a second displacement-maximizing solution
without imposing bounds. By changing the link length ratio to $\eta=1/2$, a
second solution emerged to the unbounded problem, that corresponds to a second
gait. Both gaits for the swimmer with $\eta=1/2$ are shown in Fig. 5(b).

Applying PMP reveals maximal displacement gaits for the swimmer models.
However, sometimes it fails to find an optimal solution, unless appropriate
bounds on joint angles are incorporated, which enable finding optimal gait
trajectories containing pieces along the bounds. In order to get a better
understanding of the influence of joint angle bounds on the existence of
displacement-maximizing gaits, in the following section we re-examine the
kinematic swimmer models using a geometric approach.

##  IV Geometric analysis of gaits

In order to have better insight into the behaviour of the kinematic models and
understand the changes in the optimal gaits, we now review a geometric
approach to the analysis of such swimmers [63]. This approach uses a
generalization of Stokes’ theorem to capture the relationship between the
geometry of a gait’s trajectory through shape space and the displacement it
produces. For the purposes of this paper, a key result of the geometric
approach is that the maximum-displacement-per-cycle gaits described in the
previous section approximately follow the zero-contour of the _curvature_
(augmented curl) of the local connection, such that the shape of the zero
contour can be used to predict the existence and form of the maximum-
displacement-per-cycle gaits.

++---- Figure 6: The height function plot for Purcell’s swimmer. Positive
regions are in shades of blue and negative regions in shades of yellow. The
zero level curve is the dashed black line.

###  IV-A Using Stokes’ Theorem for Measuring Net Displacement

As a preliminary to discussing the geometric locomotion formulation, we first
give a brief review of Stokes’ theorem. Let $S$ be an oriented smooth surface
that is bounded by a simple, closed, smooth boundary curve $\Gamma$ with
positive orientation. Also, let $\mathbf{F}$ be a vector field. Stokes’
theorem states that the line integral along the closed curve $\Gamma$ on the
vector field $\mathbf{F}$ is equal to the integral of the curl of that vector
field over a surface bounded by the curve,

$\ointctrclockwise_{\Gamma}\mathbf{F}\cdot
d\mathbf{r}=\iint_{S}\text{curl}\,\mathbf{F}\cdot d\mathbf{S}$ (31)

Intuitively, $\text{curl}\,\mathbf{F}$ measures how much $\mathbf{F}$ changes
“across” the gait, such that the contribution from flowing along one side of
the gait is not undone by flowing back along the opposite side of the gait.

Next, we show how the same logic as underlies Stokes’ theorem can be applied
in order to approximate the net displacement of the swimmer. The dynamic
equations of motion for the three-link swimmers, as given in (2), are

$\dot{\mathbf{q}}=\mathbf{G}\left(\boldsymbol{\upphi},\theta\right)\dot{\boldsymbol{\upphi}}=\mathbf{R}(\theta)\mathbf{A}\left(\boldsymbol{\upphi}\right)\dot{\boldsymbol{\upphi}}$
(32)

The net displacement of the swimmer over a period is equal to the line
integral over the gait,

$\Delta\mathbf{q}=\ointctrclockwise_{\boldsymbol{\upphi}}\mathbf{R}(\theta)\mathbf{A}\left(\boldsymbol{\upphi}\right)$
(33)

This line integral can be approximately converted to a surface integral in a
manner similar to Stokes’ theorem by evaluating the curvature of the local
connection over the surface $\boldsymbol{\upphi}_{a}$ as in [64], giving the
“corrected Body Velocity Integral” (cBVI)

$\Delta\mathbf{q}\approx\exp{\iint_{\boldsymbol{\upphi}_{a}}D\mathbf{A}}=\exp{\iint_{\boldsymbol{\upphi}_{a}}\mathrm{d}\mathbf{A}-[\mathbf{A}_{1},\mathbf{A}_{2}]}$
(34)

Here, $\mathrm{d}\mathbf{A}$ is the exterior derivative of the local
connection (the generalized row-wise curl), and
$[\mathbf{A}_{1},\mathbf{A}_{2}]$ is a local Lie bracket term that corrects
for noncommutativity as the swimmer translates and rotates through space (e.g.
if the system rotates counterclockwise before moving forward and/or clockwise
before moving backward, then the influence of the $\mathbf{R}(\theta)$ term
means that the “forward” and “backward” motions will be tilted slightly
“leftward”.) This bias is captured to first order by the Lie bracket, which
evaluates for planar translation and rotation as

$[\mathbf{A}_{1},\mathbf{A}_{2}]=\begin{bmatrix}\mathbf{A}_{1}^{y}\mathbf{A}_{2}^{\theta}-\mathbf{A}_{2}^{y}\mathbf{A}_{1}^{\theta}\\\
\mathbf{A}_{2}^{x}\mathbf{A}_{1}^{\theta}-\mathbf{A}_{1}^{x}\mathbf{A}_{2}^{\theta}\\\
0\end{bmatrix}$ (35)

The $\exp$ operator in (34) converts the cBVI from an “average body velocity”
to a displacement by flowing with the cBVI as body velocity for unit time. If
the rotational component of the cBVI is zero, the components of
$\exp{\text{cBVI}}$ are the same as those of the cBVI; if the rotational
component is non-zero, then the exponential flow is curved commensurately,
approximating the net rotation as having been evenly distributed along the
system’s motion.

###  IV-B Height Functions and Displacement-Per-Cycle Maximizing Gaits

Using the generalized Stokes’ theorem , we have a good approximation of the
net displacement of the swimmer over a gait. By plotting the x-integrand in
(34) as a height function
$H_{x}\left(\boldsymbol{\upphi}\right)=D\mathbf{A}^{x}(\boldsymbol{\upphi})$,
we can identify sign-definite areas of the shape plane. A positive area that
is encircled by the gait (in a counter-clockwise direction) will result in
positive net displacement of the swimmer in the $x$ direction. In Fig. 6 a
contour plot of the height function for the three-link swimmer is presented.
In order to maximize the displacement of the swimmer, the gait should enclose
a region of the height function that is as sign-definite as possible.
Obviously, this will be accomplished by following the zero level curve that
separates the positive and negative regions. The zero level curve for the
symmetric three-link swimmer, represented in Fig. 6 by a dashed closed curve,
is functionally identical to the optimal gait found using PMP in [38] and
numerically in [39].

###  IV-C Minimum Perturbation Body Frames

Note that our geometric approach involves two approximations for handling the
noncommutativity induced by the $\mathbf{R}(\theta)$ term: the Lie bracket
that makes a first-order approximation of the results of mixing _intermediate_
translation and rotation within a cycle by identifying whether
forward/backward motion was on-average preceded by clockwise/counterclockwise
motion; and the exponential map that takes the _net_ rotation over a cycle as
having been evenly distributed over the course of the motion. If these
approximations are not accurate, then the approximated net displacement will
have some residual error.

A core finding from our prior work [23, 24, 65] is that the choice of body
frame coordinates affects the magnitude of this residual error, and in
particular that the error correlates to the norm of the local connection
$\mathbf{A}$. For example, body frames attached to outer links have large
$\mathbf{A}$ values and so experience large motions in response to shape
changes. These motions accrue considerable error, both because they rotate
outside of the $\sin(\theta)\approx\theta$ linear approximation for rotation
and because they have significant translations while in these heavily-rotated
configurations. Conversely, frames such as those attached to the center of
mass and oriented to align with the mean of the link angles have small
$\mathbf{A}$ values, and this experience small motions that stay within the
valid domains of the approximations. (Note that the “large” and small” motions
under mentioned here are the intermediate motions; the net rigid body motion
is independent of the choice of body frame.)

For the examples in this paper, we placed the body frame at a location
$\bar{\mathbf{q}}_{b}=(\bar{x},\bar{y},\bar{\theta})$ found through our
algorithm in [23, 24], which identifies a weighted average of the link
positions and orientations that variationally minimizes the average norm of
$\mathbf{A}$ over the region of the shape space under consideration.

##  V Explaining changes in optimal gaits using geometric analysis

After reviewing the geometric method for gait analysis, we can now revisit the
cases where optimal gaits derived by PMP analysis depend on varying the
swimmer’s parameters and/or state bounds. By using the minimum perturbation
coordinates and plotting the height function for the models discussed here, we
get a better understanding of the way changes in the gait affect the
displacement. Moreover, we can see that topological changes in the zero-level
curves induce formation of junction points for which the variational equations
underlying PMP method become singular.

###  V-A “Reverse optimum” for Purcell’s swimmer.

First, we revisit the case of the symmetric three-link swimmer and show a
contour plot of the height function $H_{x}$ under the minimum-perturbation
coordinates. Fig. 7(a) shows the height function for the perfect fluid swimmer
at larger amplitudes of up to $6$ [rad]. In this plot, an additional zero
level curve can be seen in amplitudes of around $3$ [rad]. The junctions,
which can clearly be seen on the curve, explain the failure of the PMP method
to find an optimal gait. These junctions mean that there is no unique solution
to the unbounded optimal control problem presented.

###  V-B The perfect fluid swimmer.

For the perfect fluid model we found two gaits with opposite displacement
directions for a link length ratio of $\eta=1/2$, but were unable to find a
second, reverse direction, optimal gait for a ratio of $\eta=1/3$ without
applying bounds to the joint angles. Height functions for both cases are
presented in Fig. 7. In Fig. 7(c), showing the height function for $\eta=1/2$,
the two zero-level curves can be seen in black dashed lines. These curves
correspond to the two optimal gaits found using PMP (Fig. 5(b)). For
$\eta=1/3$, the height function in Fig. 7(a) shows the first zero level curve
representing the optimal gait which was found using PMP (Fig. 5(a)), while the
second zero level curve has junctions, once again explaining the failure of
PMP to find the second optimal gait for the unbounded problem.

Using the geometric method as presented here to plot the height function for
the displacement of a swimmer gives insight on the optimal gait and offers an
explanation for the cases where the PMP method fails to find an optimal gait,
as well as topological changes in the optimal gaits.

(a)

(b)

(c)

(d)

Figure 7: Height functions for “perfect fluid” swimmer with $\alpha=0.2$ (a)
$\eta=1/3$. Junctions can be seen in the larger zero level curves (dashed
black lines). (b) The corresponding PMP solutions. (c) $\eta=1/2$. Two loops
of zero level curves can be seen, corresponding to the two optimal gaits found
via PMP . (d) Corresponding PMP solutions.

##  VI conclusion

In this paper we examined and compared different methods of finding
displacement-maximizing gaits for two models of 3-link swimmers. We presented
the two swimmer models: Purcell’s three-link swimmer, and the ‘perfect fluid’
swimmer. We formulated the optimal control problem for a general case of a
three link swimmer with two joint angle inputs. We presented the solution
using PMP and also considered the case where there are bounds on the joint
angles. In some cases, we found that the PMP method fails to find the maximum
displacement gait for unbounded joint angles. We then reviewed the geometric
analysis approach to finding such gaits. Observing the results from this
method offered an explanation for the failures of PMP and to the unexplained
changes in the optimal gait. Junctions seen in the zero level curve mean there
is no unique solution to the optimal control problem, which causes PMP to
diverge. Adding bounds to the problem allowed us to avoid these junctions. The
geometric method gives insight into the changes in the distance-optimal gait
due to changes in the swimmer model and parameters. We now briefly discuss
limitations of our work and sketch some possible directions for future
extensions of the research. First, the methods used in this paper for
displacement-maximizing gaits can be applied to finding energy-optimal gaits
as well [25, 39]. Second, the work may be extended to time-dependent swimmer
models, such as a swimmer with elastic joints [66] and a magnetically actuated
swimmer [67], [68, 69]. Swimmers with multiple joints are another avenue of
research [70, 71], [72]. In [73] energy-optimal gaits were found for an
$N$-link swimmer assuming small amplitudes of the joint angles. This could
possibly be extended to large amplitudes by applying PMP method. Finally, the
results found here can be applied to real robotic micro swimmers in order to
verify the validity of the results [74, 59, 69].

Acknowledgement: The authors wish to thank the anonymous reviewers for their
useful and important suggestions. YO also thanks Prof. Joseph Ben-Asher from
Technion for guidance and fruitful discussions.

## References

  * [1] W. Gao, D. Kagan, O. S. Pak, C. Clawson, S. Campuzano, E. Chuluun-Erdene, E. Shipton, E. E. Fullerton, L. Zhang, E. Lauga, and J. Wang, “Cargo-towing fuel-free magnetic nanoswimmers for targeted drug delivery,” _Small_ , vol. 8, no. 3, pp. 460–467, 2012.
  * [2] JianFeng and S. Cho, “Mini and micro propulsion for medical swimmers,” _Micromachines_ , vol. 5, no. 1, pp. 97–113, 2014. [Online]. Available: http://dx.doi.org/10.3390/mi5010097
  * [3] G. Kósa, P. Jakab, N. Hata, F. Jólesz, Z. Neubach, M. Shoham, M. Zaaroor, and G. Székely, “Flagellar swimming for medical micro robots: theory, experiments and application,” in _2nd IEEE RAS & EMBS International Conference on Biomedical Robotics and Biomechatronics (BioRob)_. IEEE, 2008, pp. 258–263.
  * [4] J. M. Mirats Tur and W. Garthwaite, “Robotic devices for water main in-pipe inspection: A survey,” _Journal of Field Robotics_ , vol. 27, no. 4, pp. 491–508, 2010. [Online]. Available: https://onlinelibrary.wiley.com/doi/abs/10.1002/rob.20347
  * [5] J. Happel and H. Brenner, _Low Reynolds Number Hydrodynamics_. Prentice-Hall, 1965.
  * [6] N. Cohen and J. H. Boyle, “Swimming at low reynolds number: a beginners guide to undulatory locomotion,” _Contemporary Physics_ , vol. 51, no. 2, pp. 103–123, 2010. [Online]. Available: https://doi.org/10.1080/00107510903268381
  * [7] J. Koiller, K. Ehlers, and R. Montgomery, “Problems and progress in microswimming,” _Journal of Nonlinear Science_ , vol. 6, no. 6, pp. 507–541, 1996.
  * [8] E. M. Purcell, “Life at low Reynolds number,” _American journal of physics_ , vol. 45, no. 1, pp. 3–11, 1977.
  * [9] L. E. Becker, S. A. Koehler, and H. A. Stone, “On self-propulsion of micro-machines at low Reynolds number: Purcell’s three-link swimmer,” _Journal of Fluid Mechanics_ , vol. 490, pp. 15–35, 2003.
  * [10] R. G. Cox, “The motion of long slender bodies in a viscous fluid part 1. general theory,” _Journal of Fluid Mechanics_ , vol. 44, pp. 791–810, 1970\.
  * [11] G. Batchelor, “Slender-body theory for particles of arbitrary cross-section in Stokes flow,” _Journal of Fluid Mechanics_ , vol. 44, no. 3, pp. 419–440, 1970.
  * [12] R. L. Hatton, T. Dear, and H. Choset, “Kinematic cartography and the efficiency of viscous swimming,” _IEEE Transactions on Robotics_ , vol. 33, no. 3, pp. 523–535, Jun. 2017. [Online]. Available: https://doi.org/10.1109/tro.2017.2653810
  * [13] E. Gutman and Y. Or, “Symmetries and gaits for Purcell’s three-link microswimmer model,” _IEEE Transactions on Robotics_ , vol. 32, no. 1, pp. 53–69, 2016.
  * [14] E. Kanso, J. E. Marsden, C. W. Rowley, and J. B. Melli-Huber, “Locomotion of articulated bodies in a perfect fluid,” _Journal of Nonlinear Science_ , vol. 15, no. 4, pp. 255–289, 2005.
  * [15] J. B. Melli, C. W. Rowley, and D. S. Rufat, “Motion planning for an articulated body in a perfect planar fluid,” _SIAM Journal on applied dynamical systems_ , vol. 5, no. 4, pp. 650–669, 2006.
  * [16] H. Lamb, _Hydrodynamics_. Cambridge university press, 1993.
  * [17] K. A. Morgansen, B. I. Triplett, and D. J. Klein, “Geometric methods for modeling and control of free-swimming fin-actuated underwater vehicles,” _IEEE Transactions on Robotics_ , vol. 23, no. 6, pp. 1184–1199, 2007.
  * [18] E. Virozub, O. Wiezel, A. Wolf, and Y. Or, “Planar multi-link swimmers: Experiments and theoretical investigation using “perfect fluid” model,” _Robotica_ , vol. 37, no. 8, pp. 1289–1301, 2019.
  * [19] R. L. Hatton and H. Choset, “Geometric swimming at low and high reynolds numbers,” _IEEE Transactions on Robotics_ , vol. 29, no. 3, pp. 615–624, 2013.
  * [20] H. Marvi, C. Gong, N. Gravish, H. Astley, M. Travers, R. L. Hatton, J. R. Mendelson, H. Choset, D. L. Hu, and D. I. Goldman, “Sidewinding with minimal slip: Snake and robot ascent of sandy slopes,” _Science_ , vol. 346, no. 6206, pp. 224–229, 2014. [Online]. Available: https://science.sciencemag.org/content/346/6206/224
  * [21] T. Yona and Y. Or, “The wheeled three-link snake model: singularities in nonholonomic constraints and stick–slip hybrid dynamics induced by coulomb friction,” _Nonlinear Dynamics_ , vol. 95, no. 3, pp. 2307–2324, 2019.
  * [22] E. A. Shammas, H. Choset, and A. A. Rizzi, “Geometric motion planning analysis for two classes of underactuated mechanical systems,” _The International Journal of Robotics Research_ , vol. 26, no. 10, pp. 1043–1073, 2007\.
  * [23] R. L. Hatton and H. Choset, “Geometric motion planning: The local connection, stokes’ theorem, and the importance of coordinate choice,” _The International Journal of Robotics Research_ , vol. 30, no. 8, pp. 988–1014, 2011\.
  * [24] R. Hatton and H. Choset, “Nonconservativity and noncommutativity in locomotion,” _The European Physical Journal Special Topics_ , vol. 224, no. 17-18, pp. 3141–3174, Dec. 2015. [Online]. Available: https://doi.org/10.1140/epjst/e2015-50085-y
  * [25] S. Ramasamy and R. L. Hatton, “The geometry of optimal gaits for drag-dominated kinematic systems,” _IEEE Transactions on Robotics_ , vol. 35, no. 4, pp. 1014–1033, 2019.
  * [26] L. S. Pontryagin, _Mathematical theory of optimal processes_. CRC press, 1987.
  * [27] A. E. Bryson and Y. C. Ho, _Applied optimal control: optimization, estimation, and control_. Hemisphere, New York, 1975.
  * [28] D. Liberzon, _Calculus of Variations and Optimal Control Theory: A Concise Introduction_. Princeton University Press, 2011. [Online]. Available: https://doi.org/10.1515/9781400842643
  * [29] F. Clarke, _Functional analysis, calculus of variations and optimal control_. Springer, 2013, vol. 264.
  * [30] J. Z. Ben-Asher, _Optimal control theory with aerospace applications_. American institute of aeronautics and astronautics, 2010.
  * [31] R. Serra, C. Yanez, and C. Frueh, “Tracklet-to-orbit association for maneuvering space objects using optimal control theory,” _Acta Astronautica_ , vol. 181, pp. 271–281, 2021.
  * [32] M. R. Niknam, H. Kheiri, and N. Abdi Sobouhi, “Optimal control of satellite attitude and its stability based on quaternion parameters,” _Computational Methods for Differential Equations_ , vol. 10, no. 1, pp. 168–178, 2022.
  * [33] B. Kimia, A. Tannenbaum, and S. Zucker, “On optimal control methods in computer vision and image processing,” in _Geometry-Driven Diffusion in Computer Vision_. Springer, 1994, pp. 307–338.
  * [34] S. Rosa and D. F. Torres, “Optimal control of a fractional order epidemic model with application to human respiratory syncytial virus infection,” _Chaos, Solitons & Fractals_, vol. 117, pp. 142–149, 2018.
  * [35] J. A. Sharp, K. Burrage, and M. J. Simpson, “Implementation and acceleration of optimal control for systems biology,” _Journal of the Royal Society Interface_ , vol. 18, no. 181, p. 20210241, 2021.
  * [36] B. Bonnard and D. Sugny, _Optimal control with applications in space and quantum dynamics_. American Institute of Mathematical Sciences, 2012.
  * [37] U. Boscain, M. Sigalotti, and D. Sugny, “Introduction to the pontryagin maximum principle for quantum optimal control,” _PRX Quantum_ , vol. 2, no. 3, p. 030203, 2021.
  * [38] O. Wiezel and Y. Or, “Using optimal control to obtain maximum displacement gait for Purcell’s three-link swimmer,” in _Decision and Control (CDC), 2016 IEEE 55th Conference on_. IEEE, 2016, pp. 4463–4468.
  * [39] D. Tam and A. E. Hosoi, “Optimal stroke patterns for Purcell’s three-link swimmer,” _Physical Review Letters_ , vol. 98, no. 6, p. 068105, 2007.
  * [40] S. S. Chang, _Optimal control in bounded phase space_. New York University, College of Engineering, 1961, no. 1238.
  * [41] A. Y. Dubovitskii and A. A. Milyutin, “Necessary conditions for a weak extremum in optimal control problems with mixed constraints of the inequality type,” _USSR Computational Mathematics and Mathematical Physics_ , vol. 8, no. 4, pp. 24–98, 1968.
  * [42] J. L. Speyer and A. E. Bryson, “Optimal programming problems with a bounded state space.” _AIAA Journal_ , vol. 6, no. 8, pp. 1488–1491, 1968. [Online]. Available: https://doi.org/10.2514/3.4793
  * [43] R. F. Hartl, S. P. Sethi, and R. G. Vickson, “A survey of the maximum principles for optimal control problems with state constraints,” _SIAM Review_ , vol. 37, no. 2, pp. 181–218, 1995. [Online]. Available: https://doi.org/10.1137/1037043
  * [44] J. Ritzmann, A. Christon, M. Salazar, and C. Onder, “Fuel-optimal power split and gear selection strategies for a hybrid electric vehicle,” SAE Technical Paper, Tech. Rep., 2019.
  * [45] Y. Zhou, E. Chung, A. Bhaskar, and M. E. Cholette, “A state-constrained optimal control based trajectory planning strategy for cooperative freeway mainline facilitating and on-ramp merging maneuvers under congested traffic,” _Transportation Research Part C: Emerging Technologies_ , vol. 109, pp. 321–342, 2019.
  * [46] J. Z. Ben-Asher and E. D. Rimon, “Time optimal trajectories for a car-like mobile robot,” _IEEE Transactions on Robotics_ , 2021.
  * [47] B.-H. Nguyen, T. Vo-Duy, C. H. Antunes, and J. P. F. Trovão, “Multi-objective benchmark for energy management of dual-source electric vehicles: An optimal control approach,” _Energy_ , vol. 223, p. 119857, 2021\.
  * [48] S. Lee, Y. Kim, J. B. Siegel, and A. G. Stefanopoulou, “Optimal control for fast acquisition of equilibrium voltage for li-ion batteries,” _Journal of Energy Storage_ , vol. 40, p. 102814, 2021.
  * [49] O. Cots, J. Gergaud, and D. Goubinat, “Direct and indirect methods in optimal control with state constraints and the climbing trajectory of an aircraft,” _Optimal Control Applications and Methods_ , vol. 39, no. 1, pp. 281–301, 2018.
  * [50] K. Mall and E. Taheri, “Unified trigonometrization method for solving optimal control problems in atmospheric flight mechanics,” in _AIAA Scitech 2020 Forum_ , 2020, p. 0022.
  * [51] J. Holtzman, “Convexity and the maximum principle for discrete systems,” _IEEE Transactions on Automatic Control_ , vol. 11, no. 1, pp. 30–35, 1966\.
  * [52] V. G. Bolt’yanskii, “The method of tents in the theory of extremal problems,” _Russian Mathematical Surveys_ , vol. 30, no. 3, p. 1, 1975.
  * [53] A. Y. Dubovitskii, “The discrete principle of the maximum,” _Avtomatika i Telemekhanika_ , no. 10, pp. 55–71, 1978.
  * [54] M. B. Kobilarov and J. E. Marsden, “Discrete geometric optimal control on lie groups,” _IEEE Transactions on Robotics_ , vol. 27, no. 4, pp. 641–655, 2011\.
  * [55] A. Saccon, J. Hauser, and A. P. Aguiar, “Optimal control on lie groups: The projection operator approach,” _IEEE Transactions on Automatic Control_ , vol. 58, no. 9, pp. 2230–2245, 2013.
  * [56] K. S. Phogat, D. Chatterjee, and R. N. Banavar, “A discrete-time pontryagin maximum principle on matrix lie groups,” _Automatica_ , vol. 97, pp. 376–391, 2018.
  * [57] S. Kadam, K. S. Phogat, R. N. Banavar, and D. Chatterjee, “Exact isoholonomic motion of the planar purcell’s swimmer,” _IEEE Transactions on Automatic Control_ , vol. 67, no. 1, pp. 429–435, 2021.
  * [58] A. M. Bloch, “Nonholonomic mechanics,” in _Nonholonomic mechanics and control_. Springer, 2003, pp. 207–276.
  * [59] J. Wu, B. Jang, Y. Harduf, Z. Chapnik, O. B. Avci, X. Chen, J. Puigmartí-Luis, O. Ergeneman, B. J. Nelson, Y. Or, and S. Pané, “Helical klinotactic locomotion of two-link nanoswimmers with dual-function drug-loaded soft polysaccharide hinges,” _Advanced Science_ , vol. n/a, no. n/a, p. 2004458, 2021. [Online]. Available: https://onlinelibrary.wiley.com/doi/abs/10.1002/advs.202004458
  * [60] O. Wiezel and Y. Or, “Optimization and small-amplitude analysis of Purcell’s three-link microswimmer model,” _Proc. R. Soc. A_ , vol. 472, no. 2192, p. 20160425, 2016.
  * [61] L. N. Hand and J. D. Finch, _Analytical mechanics_. Cambridge University Press, 1998.
  * [62] S. D. Kelly and R. M. Murray, “Geometric phases and robotic locomotion,” _Journal of Robotic Systems_ , vol. 12, no. 6, pp. 417–431, 1995.
  * [63] R. L. Hatton and H. Choset, “Optimizing coordinate choice for locomoting systems,” in _2010 IEEE International Conference on Robotics and Automation_. IEEE, 2010, pp. 4493–4498.
  * [64] S. Ramasamy and R. L. Hatton, “Soap-bubble optimization of gaits,” in _2016 IEEE 55th Conference on Decision and Control (CDC)_. IEEE, 2016, pp. 1056–1062.
  * [65] C. Bass, S. Ramasamy, and R. L. Hatton, “Characterizing error in noncommutative geometric gait analysis,” in _2022 International Conference on Robotics and Automation (ICRA)_. IEEE, 2022, pp. 9845–9851.
  * [66] E. Passov (Gutman) and Y. Or, “Dynamics of Purcell’s three-link microswimmer with a passive elastic tail,” _The European Physical Journal E_ , vol. 35, no. 8, pp. 1–9, 2012.
  * [67] E. Gutman and Y. Or, “Simple model of a planar undulating magnetic microswimmer,” _Physical Review E_ , vol. 90, no. 1, p. 013012, 2014.
  * [68] M. Zoppello, A. De Simone, F. Alouges, and L. Giraldi, “Modeling and steering magneto-elastic micro-swimmers inspired by the motility of sperm cells,” _Atti della Accademia Peloritana dei Pericolanti-Classe di Scienze Fisiche, Matematiche e Naturali_ , vol. 96, no. S3, p. 12, 2018.
  * [69] Y. El Alaoui-Faris, J.-B. Pomet, S. Régnier, and L. Giraldi, “Optimal actuation of flagellar magnetic microswimmers,” _Physical Review E_ , vol. 101, no. 4, p. 042604, 2020.
  * [70] L. Giraldi, P. Martinon, and M. Zoppello, “Controllability and optimal strokes for n-link microswimmer,” in _52nd IEEE Conference on Decision and Control_. IEEE, 2013, pp. 3870–3875.
  * [71] R. Marchello, M. Morandotti, H. Shum, and M. Zoppello, “The N-link swimmer in three dimensions: Controllability and optimality results,” _Acta Applicandae Mathematicae_ , vol. 178, no. 1, pp. 1–25, 2022.
  * [72] S. Ramasamy and R. L. Hatton, “Geometric gait optimization beyond two dimensions,” in _2017 American Control Conference (ACC)_ , 2017, pp. 642–648.
  * [73] F. Alouges, A. DeSimone, L. Giraldi, Y. Or, and O. Wiezel, “Energy-optimal strokes for multi-link microswimmers: Purcell’s loops and taylor’s waves reconciled,” _New Journal of Physics_ , vol. 21, no. 4, p. 043050, 2019.
  * [74] B. Jang, E. Gutman, N. Stucki, B. F. Seitz, P. D. Wendel-García, T. Newton, J. Pokki, O. Ergeneman, S. Pané, Y. Or _et al._ , “Undulatory locomotion of magnetic multilink nanoswimmers,” _Nano Letters_ , vol. 15, no. 7, pp. 4829–4833, 2015.


# Heat transport in a flowing complex plasma in microgravity conditions

V. Nosenko<EMAIL_ADDRESS>Institut für Materialphysik im Weltraum,
Deutsches Zentrum für Luft- und Raumfahrt (DLR), D-82234 Weßling, Germany
Center for Astrophysics, Space Physics, and Engineering Research, Baylor
University, Waco, Texas 76798-7310, USA  S. Zhdanov  Institut für
Materialphysik im Weltraum, Deutsches Zentrum für Luft- und Raumfahrt (DLR),
D-82234 Weßling, Germany  M. Pustylnik  Institut für Materialphysik im
Weltraum, Deutsches Zentrum für Luft- und Raumfahrt (DLR), D-82234 Weßling,
Germany  H. M. Thomas  Institut für Materialphysik im Weltraum, Deutsches
Zentrum für Luft- und Raumfahrt (DLR), D-82234 Weßling, Germany  A. M. Lipaev
Institute for High Temperatures, Russian Academy of Sciences, Izhorskaya
13/19, 125412 Moscow, Russia  O. V. Novitskii  Gagarin Research and Test
Cosmonaut Training Center, 141160 Star City, Moscow Region, Russia

###### Abstract

Heat transport in a three-dimensional complex (dusty) plasma was
experimentally studied in microgravity conditions using Plasmakristall-4
(PK-4) instrument on board the International Space Station (ISS). An extended
suspension of microparticles was locally heated by a shear flow created by
applying the radiation pressure force of the manipulation-laser beam.
Individual particle trajectories in the flow were analysed and from these,
using a fluid heat transport equation that takes viscous heating and neutral
gas drag into account, the complex plasma’s thermal diffusivity and kinematic
viscosity were calculated. Their values are compared with previous results
reported in ground-based experiments with complex plasmas.

###### pacs:

52.27.Lw, 52.27.Gr, 66.60.+a

##  I Introduction

Transport phenomena in liquids are important in fundamental science and in
engineering applications. They are, however, challenging to study
experimentally. To perform experiments at the most fundamental kinetic level
(of individual particles), model systems are necessary. A particularly
convenient model system is complex plasma. This is a suspension of micron size
particles in ionized gas [1]. The particles become charged by collecting
electrons and ions from the plasma and interact with each other via a screened
Coulomb potential which in many cases can be approximated by the Yukawa
potential. The advantages of complex plasmas as model systems are the relative
ease of individually observing the constituent particles, stretched spatial
and temporal scales, and undamped dynamics. Complex plasmas were successfully
used to study diffusion [2], momentum [3, 4, 5, 6, 7], and heat transport [8,
9, 10]. Heat transport was experimentally studied in a crystalline phase of 2D
complex plasma [8], in a 2D complex plasma at melting conditions [9], and in a
3D liquid complex plasma [10]. So far, the heat transport in complex plasmas
was studied in ground-based laboratory experiments, where the particle
suspension was compressed by the force of gravity.

In this paper, we present an experiment with liquid complex plasma, where the
heat transport was studied by locally heating the particle suspension by a
laser-driven shear flow and analyzing the resultant particle kinetic
temperature profile. To achieve an extended unstressed 3D particle suspension,
the experiments were performed in microgravity conditions. One advantage of
studying transport phenomena in a 3D sample of complex plasma is that the
existence of valid transport coefficients is not questioned here, unlike in a
2D sample [11].

##  II Experimental method

The experiments described in this paper were performed using the PK-4
instrument on board the International Space Station (ISS) [12]. The
experimental procedure and parameters were similar to those used in Ref. [7].
We used DC discharge to produce plasma in Ne at the pressure of $15$ Pa. The
discharge current was $0.5$ mA and the maximum voltage was $2.5$ kV. Estimates
are available of the electron density $n_{e}=0.92\times 10^{8}$ cm-3 and
temperature $T_{e}=9.8$ eV (on the tube axis in the middle of the working
area) in a similar discharge in Ne at the pressure of $20$ Pa, see Ref. [12]
for more details. The screening length in these conditions was
$\lambda_{D}=122~{}\mu$m. A cloud of melamine formaldehyde (MF) microspheres
with a diameter of $3.38\pm 0.07$ was suspended in the working area of the
discharge tube. They were trapped using polarity switching with a frequency of
$500$ Hz and duty cycle $\simeq 50$%. The Wigner-Seitz radius of the particle
suspension was calculated as $r_{\rm WS}=r_{0}/1.79=100~{}\mu$m, where $r_{0}$
is the first peak position of the pair correlation function $g(r)$ measured in
2D cross sections of the particle suspension [13]. The particle number density
was calculated as $n=3(4\pi r_{\rm WS}^{3})^{-1}=2.4\times 10^{5}$ cm-3. For
the particle charge, we adopted the value $Q=2800e$ reported for our
experimental conditions in Ref. [14]. The screening parameter was
$\kappa=r_{\rm WS}/\lambda_{D}=0.8$.

A shear flow in the particle suspension was created by applying the radiation
pressure force from the focused beam of a $808$-nm manipulation laser. The
laser beam was aligned with the discharge tube axis and its output power was
set to $2.16$ W. As in Ref. [7], the main experimental run was supplemented by
a plasma-off phase, where the plasma was briefly (during $0.5$ s) switched off
while the manipulation laser remained on. During the plasma-off time, the
particle charge rapidly declined [15], as did the interparticle interactions
and external confinement; however, the particle suspension did not collapse
due to weightlessness. The particles were illuminated by a thin sheet of a
$532$-nm laser and imaged by two video cameras with slightly overlapping
fields of view. Here, we used the video recorded by camera 2. To better
resolve the particle motion, a higher than in Ref. [7] camera frame rate of
$100$ frames per second was used at the expense of narrower field of view
(FoV), which measured $1600\times 267~{}{\rm pixels}^{2}$ or $22.72\times
3.79~{}{\rm mm}^{2}$. The reduced FoV, however, included the full width of the
shear flow with boundary layers. The manipulation-laser beam was coming in the
positive $x$ direction, the particles were imaged in the $(x,z)$ plane, and in
the scanning direction the imaging system was fixed at $y=0$.

Figure 1:  (a) Trajectories of individual particles in the central cross
section ($y=0$) of the steady-state laser-induced shear flow during $0.3$ s.
The laser radiation force pushes the particles in the positive $x$ direction.
(b) Time-averaged particle velocity profile $v_{x}(z)$. (c) Map of the time-
averaged vorticity $(\nabla\times{\bf v})_{y}$, which shows the spatial
localisation of the shear in the particle flow. In (b) and (c), time averaging
was performed during $1.2$ s. Shown here is a part of the original FoV
delimited by $400~{}{\rm pixels}<x<1200~{}{\rm pixels}$.

Experimental data were analysed using the following method. A $1.2$-s section
of the experimental video recording of the steady-state shear flow was
identified and used for analysis. Since the shear flow has cylindrical
symmetry [7], it is sufficient to analyze its central cross section (defined
by $y=0$). Particles were identified and traced from frame to frame using a
moment method [16]. A part of the original FoV where the particle suspension
was rather uniform (delimited by $400~{}{\rm pixels}<x<800~{}{\rm pixels}$),
was divided into $33$ bins elongated along the $x$ axis and $8$-pixels wide in
the $z$ direction. The average number of particles was $\approx 50$ per bin
per frame, to the total of $\approx 6000$ per bin for the whole sequence. In
each bin, the distributions of the particle velocity $v_{x,z}$ were calculated
during $1.2$ s. In all bins, the distributions were close to Maxwellian
(shifted for the $v_{x}$ distribution). Therefore, the time-averaged profiles
of the particle speed $v_{x,z}(z)$ and kinetic temperature $T_{x,z}(z)$ were
calculated. The particle kinetic temperature was calculated as
$T_{x,z}=m\langle(v_{x,z}-\overline{v}_{x,z})^{2}\rangle/k_{B}$, where $m$ is
the particle mass and $k_{B}$ is the Boltzmann constant. The errors in
calculating $v_{x,z}(z)$ and $T_{x,z}(z)$ are estimated as $\pm 0.04$ mm/s and
up to $30\%$, respectively.

##  III Experimental results

The trajectories of individual particles in the central cross section of
steady-state shear flow are shown in Fig. 1(a). The flow appears to be laminar
on average with particles performing random motion on top of it. The time-
averaged particle velocity profile $v_{x}(z)$ is shown in Fig. 1(b). The map
of the time-averaged vorticity $(\nabla\times{\bf v})_{y}$, which shows the
spatial localisation of shear in the particle flow, is shown in Fig. 1(c). The
shear is localized in two rather narrow boundary layers adjacent to the flow.

In each narrow bin corresponding to one point in Fig. 1(b), the distributions
of the particle velocity $v_{x,z}$ were close to Maxwellian (shifted for
$v_{x}$). Two examples are shown in Fig. 2. The distributions shown in panel
(a) were measured in the area of maximum shear, $0.8~{}{\rm mm}<z<0.91~{}{\rm
mm}$ and those in panel (b) in the flow bulk, $1.7~{}{\rm mm}<z<1.81~{}{\rm
mm}$. The nearly Maxwellian velocity distributions validate our calculating
the mean velocities $\overline{v}_{x,z}$ and also allow us to introduce the
notion of the particle kinetic temperature
$T_{x,z}=m\langle(v_{x,z}-\overline{v}_{x,z})^{2}\rangle/k_{B}$, which
describes the random component of the particle motion.

Figure 2:  Distributions of the particle velocity $v_{x}$ (red circles) and
$v_{z}$ (blue diamonds) measured in (a) the area of maximum shear, $0.8~{}{\rm
mm}<z<0.91~{}{\rm mm}$ and (b) in the flow bulk, $1.7~{}{\rm mm}<z<1.81~{}{\rm
mm}$. The lines are (shifted) Maxwellian fits. Figure 3:  (a) Flow velocity
profiles $v_{x,z}(z)$ and (b) squared first derivative $(\partial
v_{x}(z)/\partial z)^{2}$, to which the heat source is proportional. The red
circles are for $v_{x}(z)$, blue diamonds for $v_{z}(z)$, solid symbols
correspond to the plasma-on phase, open symbols to the plasma-off phase. (c)
Longitudinal $T_{x}(z)$ (red circles), transverse $T_{z}(z)$ (blue diamonds),
and “thermalized” $T=(T_{x}+2T_{z})/3$ (pink squares) particle kinetic
temperature profiles for the plasma-on phase.

The time-averaged flow velocity profiles $v_{x,z}(z)$ are shown in Fig. 3(a).
In the plasma-on phase (shown by the solid lines), the $v_{x}(z)$ component is
approximately symmetric with respect to $z=1.65$ mm, while the $v_{z}(z)$
component is negligibly small, which supports the assumption of cylindrical
symmetry of the flow. In the plasma-off phase (shown by the dashed lines), the
$v_{x}(z)$ component does not change much (except in the central part).
However, the $v_{z}(z)$ component increases substantially due to the drift of
the particle cloud in the positive $z$ direction (presumably due to the
thermophoretic force). These observations corroborate the findings of Ref.
[7]. Since the source function of viscous heating in the flow geometry of our
experiment is proportional to $\eta(\partial v_{x}(z)/\partial z)^{2}$, where
$\eta$ is the shear viscosity of the particle suspension [17], it is
instructive to plot $(\partial v_{x}(z)/\partial z)^{2}$ as shown in Fig.
3(b). This figure serves to visualize the distribution of heat sources in the
shear flow.

The time-averaged particle kinetic temperature profiles in the plasma-on phase
$T_{x,z}(z)$ are shown in Fig. 3(c). The temperature is anisotropic with its
longitudinal component greater on average than transverse,
$T_{x}(z)>T_{z}(z)$. In what follows, we will use the “thermalized”
$T=(T_{x}+2T_{z})/3$. There are two prominent peaks in $T(z)$ corresponding to
the viscous heating sources [cf. Fig. 3(b)]. This is a direct observation of
viscous heating in a flowing liquid 3D complex plasma, which was previously
reported only in 2D complex plasmas [3, 18]. In the plasma-off phase, the
random part of the particle motion was very strong and it was not possible to
calculate their kinetic temperature.

In very simplified terms, the main idea of the subsequent analysis is to
deduce the shear viscosity $\eta$ of the particle suspension from the height
of the temperature peaks (since the heat source function is proportional to
$\eta$) and its thermal conductivity from the peak widths.

##  IV Heat transport model

The cylindrical symmetry of the flow significantly simplifies theoretical
consideration. In what follows, we will base our model on this suggestion
assuming that the main flow parameters are only dependent on the cylindrical
radius $r$, that is, $T=T(r),v_{x}=v_{x}(r)$, etc. Moreover, when plasma is on
and the radial confinement suppresses well the radial particle motion, the
radial flow component is small compared to the axial one, see Fig. 3(a),
therefore, it can be ignored in the model.

Generally, when considering heat transport in any system, one has to address
two major points: (i) the source of heat, (ii) the mechanism of heat
redistribution. In complex plasmas, the heat transport is traditionally
characterised studying the particle kinetic temperature distributions [8, 9,
10]. It is convenient to classify the sources of heat as local ones and those
globally distributed. We attribute the quite intense local heating observed in
the present experiment to the internal _viscous heating_ [17]. It is well
evidenced by the distribution of the squared derivative of the particle flow
velocity, to which this type of heat source is proportional in our case. It is
sharply peaked in the vicinity of the inflection point of the flow velocity
profile, see Fig. 3(b). Noticeably, the kinetic temperature also demonstrates
prominent maxima in the vicinity of these points, see Fig. 3(c). In contrast,
the distributed heat sources are spread approximately uniformly inside the
particle cloud. They can be attributed to external sources, e.g., those
stemming from the thermophoretic force caused by the non-uniformly heated
discharge tube walls [19] or the photophoretic force from the manipulation
laser [20].

In what follows, we demonstrate that the assumptions based on these
experimental observations significantly simplify the analysis of the heat
transport and allow us to come, with relatively simple means, to rather
important conclusions.

###  IV.1 The governing equations

The heat is transported within the particle cloud by mutual collisions which
are interpreted here as the thermal conductivity. Under steady-state
conditions of the particle flow with cylindrical symmetry, this process can be
described by the following simple model:

$b^{2}\Delta T\equiv b^{2}\frac{d}{rdr}(r\frac{dT}{dr})=T-T_{\rm bg}-S,$ (1)

where $T$ is the “thermalized” particle kinetic temperature, $T_{\rm bg}={\rm
const}$ is the background kinetic temperature due to the distributed heat
sources, while $b$ is the heat transport length defined by the cloud thermal
diffusivity $\chi$ via the relationship

$\chi=\frac{2\gamma b^{2}}{c_{p}},$ (2)

where $\gamma$ is the Epstein neutral gas drag rate [21] and $c_{p}=5/2$ is
the specific heat capacitance. Finally, $S=S(r)$ is the local heat source
proportional to the squared derivative of the flow velocity $v^{\prime 2}$,
the particle mass $m$, and the cloud kinematic viscosity $\nu=\eta(mn)^{-1}$:

$S=Cv^{\prime 2}(r)\equiv\frac{m\nu}{4\gamma}v^{\prime 2}(r).$ (3)

The values $b$, $T_{\rm bg}$, and $C$ are constant parameters that have to be
obtained by fitting the solution of the model (1),(2),(3) to experimental
data. An advantage of this model is that Eq. (1) can be solved analytically:

$\begin{split}T=T(r)=T_{\rm
bg}+I_{0}\left(\frac{r}{b}\right)\int_{r}^{\infty}\frac{d\xi\xi}{b^{2}}S(\xi)K_{0}\left(\frac{\xi}{b}\right)+\\\
K_{0}\left(\frac{r}{b}\right)\int_{0}^{r}\frac{d\xi\xi}{b^{2}}S(\xi)I_{0}\left(\frac{\xi}{b}\right),\end{split}$
(4)

where $I_{0}$ and $K_{0}$ are the modified Bessel functions of order zero.
Note an important consequence of this relationship:

$\int_{0}^{\infty}d\xi\xi\left[T(\xi)-T_{\rm
bg}\right]=C\int_{0}^{\infty}d\xi\xi v^{\prime 2}(\xi).$ (5)

Actually, the integrals in Eq. (5) can be measured experimentally (of course,
using some crude approximation for $T_{\rm bg}$), which allows one to estimate
$C$ and then viscosity by using Eq. (3). This may serve as a first step when
fitting the analytically computed kinetic temperature profiles to those
experimentally measured.

###  IV.2 Further simplifications

Since the exact analytical solution (4) is somewhat cumbersome, it is worth
considering further possible simplifications. Given quite narrow and sharply
peaked distribution of the heating term $S(r)$ [see Fig. 3(b)], Eq. (4) can be
reduced to the following form, which is well suited for fitting purposes:

$T=T(r)=T_{\rm
bg}+A\tau(r),\tau=\left\\{\begin{array}[]{ll}I_{0}\left(\frac{r}{b}\right)/I_{0}\left(\frac{R}{b}\right),&\hbox{r
$\leq$ R;}\\\
K_{0}\left(\frac{r}{b}\right)/K_{0}\left(\frac{R}{b}\right),&\hbox{r $>$
R,}\end{array}\right.$ (6)

where the parameter $R$, considered fixed in what follows, corresponds to the
radial coordinate of the inflection point well pronounced in the flow velocity
distribution, see Figs. 3(a),(b), while the parameter $A$ stands for

$A=\frac{C}{b^{2}}\int_{0}^{\infty}d\xi\xi v^{\prime 2}(\xi)$ (7)

and defines the excess of the maximum kinetic temperature $T_{\rm max}$ over
the background value $T_{\rm bg}$. The expression $T_{\rm max}=T(R)=T_{\rm
bg}+A$ becomes exact at $r=R$. Note that $T_{\rm bg}$, $A$, and
$\mathcal{R}=R/b$ represent an optimal set of fitting parameters.

###  IV.3 The scoring function

As a scoring function for performing fitting iterations, it is advantageous to
introduce

$f=f(T_{\rm bg},A,\mathcal{R})=\left<\left(T(r)-T_{\rm
bg}-A\tau(r)\right)^{2}\right>\Rightarrow\min,$ (8)

where $\left<...\right>$ denotes the average over available range of $r$. It
is convenient to introduce also the mean values $\overline{T}=\left<T\right>$,
$\overline{\tau}=\left<\tau\right>$ and for the sake of brevity the deviations
from the mean values $\widetilde{T}=T-\overline{T}$,
$\widetilde{\tau}=\tau-\overline{\tau}$; the scoring function then reduces to

$\begin{split}f=\left(\overline{T}-T_{\rm
bg}-A\overline{\tau}\right)^{2}-\frac{\left<\widetilde{\tau}\widetilde{T}\right>^{2}}{\left<\widetilde{\tau}^{2}\right>}+\left<\widetilde{T}^{2}\right>+\\\
\left<\widetilde{\tau}^{2}\right>\left(A-\frac{\left<\widetilde{T}\widetilde{\tau}\right>}{\left<\widetilde{\tau}^{2}\right>}\right)^{2},\end{split}$
(9)

with the evident minimum/optimal essentially positive value

$f=f^{\rm
opt}(\mathcal{R})=\left<\widetilde{T}^{2}\right>-\frac{\left<\widetilde{\tau}\widetilde{T}\right>^{2}}{\left<\widetilde{\tau}^{2}\right>},$
(10)

which is reached at

$A^{\rm
opt}=\frac{\left<\widetilde{T}\widetilde{\tau}\right>}{\left<\widetilde{\tau}^{2}\right>},~{}T_{\rm
bg}^{\rm opt}=\overline{T}-A^{\rm opt}\overline{\tau}.$ (11)

Figure 4:  Normalized scoring function $F=\langle(T_{\rm exp}-T_{\rm fit})^{2}\rangle/\langle(T_{\rm exp}-\langle T_{\rm exp}\rangle)^{2}$ with second-order polynomial fit. Figure 5:  Particle kinetic temperature profile $T(r)$ for the plasma-on phase with theoretical fit using Eq. 6. Blue down triangles are for $z<1.59$ mm, green up triangles for $z>1.59$ mm, and purple squares are their mean values.  Table 1: Fitting parameters, see the text for details. $R/b$ |  $A$ (K) |  $T_{\rm bg}$ (K) |  $T_{\rm max}^{\rm exp}$ (K) |  $T_{\rm max}^{\rm fit}$ (K) |  $R$ (mm)
---|---|---|---|---|---
$2.24\pm 0.11$ | $665\pm 22$ | $476\pm 16$ | $1144\pm 131$ | $1141\pm 38$ | 0.738

Finally, it is enough to further minimize the scoring function already
optimized over the $A$ and $T_{\rm bg}$ parameters, see Eq. (10), which is
quite a straightforward task, see Fig. 4. This procedure results in the
optimized set of parameters shown in Table I. Together with Eq. (6), they are
sufficient to successfully solve the fitting problem for our experimental
data, see Fig. 5.

##  V Thermal diffusivity and kinematic viscosity

An estimate of the heat transport length follows immediately from
$b=R/\mathcal{R}$ and then, with $b$ known, the thermal diffusivity can be
calculated using Eq. (2). Both estimates are collected in Table II. A fairly
high accuracy, despite quite large scatter of the experimental points, of
these estimates is worth mentioning. To estimate the kinematic viscosity, note
first that from the definitions Eqs. (4),(6) it follows

$A=Cb^{-2}I_{0}(\mathcal{R})K_{0}(\mathcal{R})\int_{0}^{\infty}d\xi\xi
v^{\prime 2}(\xi).$ (12)

Further, given that $C=\frac{m\nu}{4\gamma}$, the kinematic viscosity is easy
to estimate. The result is shown in Table II. The momentum transport scale
length $a$ is also shown there. Its value follows immediately from the simple
scaling relationship $\nu=a^{2}\gamma$. Note that, compared to the heat
transport parameters, the latter estimates are significantly less accurate.
The source of the inaccuracy is the quite shallow “scoring well”, see Fig. 4.

Table 2: Estimated complex-plasma’s thermal diffusivity $\chi$, kinematic viscosity $\nu$, momentum transport length $a$, heat transport length $b$, mean free path $\ell$, and energy exchange rate $\beta$. $\chi({\rm mm}^{2}$/s) |  $\nu({\rm mm}^{2}$/s) |  $a$ (mm) |  $b$ (mm) |  $\ell$ (mm) | $\beta({\rm s}^{-1})$
---|---|---|---|---|---
$3.3\pm 0.1$ | $0.8\pm 0.4$ | $0.15\pm 0.03$ | $0.33\pm 0.02$ | $0.01$ | $203$

The obtained value of thermal diffusivity $\chi=3.3\pm 0.1~{}{\rm mm}^{2}$/s
is between the values of $\approx 1~{}{\rm mm}^{2}$/s and $\approx 9$ mm2/s
reported for a 3D liquid complex plasma [10] and a 2D complex plasma near
melting transition [9], respectively, both in ground-based experiments. The
calculated value of kinematic viscosity $\nu=0.8\pm 0.4~{}{\rm mm}^{2}$/s is
consistent with the estimate reported in microgravity experiments [7]. On the
other hand, it is much lower than was previously reported in ground-based
experiments with 3D complex plasmas [4, 5, 6]. In Ref. [7], this difference
was attributed to the formation of strings in the particle cloud. The strings
tend to align themselves with the axis of discharge tube [22, 23, 24] and can
reduce the viscosity of complex plasma. It is interesting to note that our
result for kinematic viscosity of a 3D complex plasma is comparable to that of
a 2D complex plasma [3] $\nu_{\rm 2D}\simeq 2~{}{\rm mm^{2}/s}$ and also of
liquid water [25], $\nu_{w}\simeq 1.8~{}{\rm mm^{2}/s}$.

##  VI Energy exchange between the kinetic temperature components

An advantage of working with the “thermalized” particle kinetic temperature is
that the energy exchange terms are excluded from the heat transport equation.
Since the rate coefficient of the energy exchange is apriory unknown, this
approach significantly simplifies the consideration. Still, it would be
incorrect to ignore the exchange effects because the observed kinetic
temperature components are noticeably different, see Fig. 3(c).

The “thermalized” kinetic temperature is by definition
$T=\frac{1}{3}\left(T_{x}+T_{y}+T_{z}\right)$, where
$T_{x,y,z}=m\left<\left(v_{x,y,z}-\left<v_{x,y,z}\right>\right)^{2}\right>$
are actually the mean squared random velocity components. In our cylindrically
symmetric case, $T_{y}=T_{z}$ and this relationship simplifies to
$T=\frac{1}{3}\left(T_{\|}+2T_{\bot}\right)$, where $T_{\|}$ and
$T_{\bot}=\frac{1}{2}\left(T_{y}+T_{z}\right)$ are the longitudinal and
transverse kinetic temperature components with respect to the particle flow
direction $x$.

Let us consider the far-field region away from the intense viscous heating,
$r>1.4$ mm. In this region, the kinetic temperature components approach nearly
constant mean values, $\left<T\right>=730\pm 70$ K,
$\left<T_{\|}\right>=1100\pm 100$ K, $\left<T_{\bot}\right>=550\pm 70$ K. The
relaxed mean value, e.g., for $\left<T_{\bot}\right>$ should obey a balance
equation

$nc_{p}\frac{\partial}{\partial
t}\left<T_{\bot}\right>=0=-2n\gamma\left<T_{\bot}\right>-n\beta\left(\left<T_{\bot}\right>-\left<T\right>\right)+S_{\bot},$
(13)

where $\beta$ is the energy exchange rate defined as the inverse of the energy
exchange time, $\tau_{\varepsilon}=\beta^{-1}$, and $S_{\bot}$ symbolises the
(minor) external heating. Ignoring this factor, we get

$\beta=6\gamma\frac{\left<T_{\bot}\right>}{\left<T_{\|}\right>-\left<T_{\bot}\right>}\cong
203~{}\text{s}^{-1}.$ (14)

##  VII Concluding remarks

In conclusion, it is worth to discuss the validity of the simplifications made
in the heat transport model. First, the applicability of the fluid dynamics
approach is justified in the following way: Since the typical thermal velocity
of the particles is about $v_{T}\approx 1$ mm/s, for $\tau_{\varepsilon}\sim
1/200$ s the mean free path is of the order of
$\ell=2v_{T}\tau_{\varepsilon}\sim 10^{-2}$ mm, which is much smaller than all
other length scales in the problem. Second, the validity of the analysis using
simplified Eq. (6) follows from the fact that the momentum transport length
$a$ turns out to be smaller than the heat transport length $b$, see Table II.
It is also clear from Figs. 3(b),(c), that the temperature distribution is
wider than the heat source function (proportional to the squared derivative of
the flow velocity).

##  VIII Acknowledgments

The authors gratefully acknowledge the joint ESA-Roscosmos experiment
“Plasmakristall-4” on board the International Space Station. This work was
supported in part by DLR/BMWi Grant No. 50WM1441. We thank Ch. Knapek for
carefully reading the manuscript and helpful comments.

##  IX Author declarations

The authors have no conflicts of interest to disclose.

##  X Data availability

Raw data were generated at the PK-4 facility on the International Space
Station. Derived data supporting the findings of this study are available from
the corresponding author upon reasonable request.

## References

  * [1] A. Ivlev, H. Löwen, G. Morfill, C. P. Royall, Complex Plasmas and Colloidal Dispersions: Particle-resolved Studies of Classical Liquids and Solids, Series in Soft Condensed Matter Vol. 5 (World Scientific, Singapore, 2012).
  * [2] S. Nunomura, D. Samsonov, S. Zhdanov, and G. Morfill, Phys. Rev. Lett. 96, 015003 (2006).
  * [3] V. Nosenko and J. Goree, Phys. Rev. Lett. 93, 155004 (2004).
  * [4] A. Gavrikov, I. Shakhova, A. Ivanov, O. Petrov, N. Vorona, and V. Fortov, Phys. Lett. A 336, 378 (2005).
  * [5] N. A. Vorona, A. V. Gavrikov, A. S. Ivanov, O. F. Petrov, V. E. Fortov, and I. A. Shakhova, Journal of Experimental and Theoretical Physics 105, 824 (2007).
  * [6] A. V. Ivlev, V. Steinberg, R. Kompaneets, H. Höfner, I. Sidorenko, and G. E. Morfill, Phys. Rev. Lett. 98, 145003 (2007).
  * [7] V. Nosenko, M. Pustylnik, M. Rubin-Zuzic, A. M. Lipaev, A. V. Zobnin, A. D. Usachev, H. M. Thomas, M. H. Thoma, V. E. Fortov, O. Kononenko, and A. Ovchinin, Phys. Rev. Research 2, 033404 (2020).
  * [8] S. Nunomura, D. Samsonov, S. Zhdanov, and G. Morfill, Phys. Rev. Lett. 95, 025003 (2005).
  * [9] V. Nosenko, S. Zhdanov, A. V. Ivlev, G. Morfill, J. Goree, and A. Piel, Phys. Rev. Lett. 100, 025003 (2008).
  * [10] V. E. Fortov, O. S. Vaulina, O. F. Petrov, M. N. Vasiliev, A. V. Gavrikov, I. A. Shakova, N. A. Vorona, Yu. V. Khrustalyov, A. A. Manohin, and A. V. Chernyshev, Phys. Rev. E 75, 026403 (2007).
  * [11] M. H. Ernst, E. H. Hauge, and J. M. J. van Leeuwen, Phys. Rev. Lett. 25, 1254 (1970).
  * [12] M. Y. Pustylnik, M. A. Fink, V. Nosenko, T. Antonova, T. Hagl, H. M. Thomas, A. V. Zobnin, A. M. Lipaev, A. D. Usachev, V. I. Molotkov, O. F. Petrov, V. E. Fortov, C. Rau, C. Deysenroth, S. Albrecht, M. Kretschmer, M. H. Thoma, G. E. Morfill, R. Seurig, A. Stettner, V. A. Alyamovskaya, A. Orr, E. Kufner, E. G. Lavrenko, G. I. Padalka, E. O. Serova, A. M. Samokutyayev, and S. Christoforetti, Rev. Sci. Instrum. 87, 093505 (2016).
  * [13] B. Liu, J. Goree, and W. D. Suranga Ruhunusiri, Rev. Sci. Instrum. 86, 033703 (2015).
  * [14] T. Antonova, S. A. Khrapak, M. Y. Pustylnik, M. Rubin-Zuzic, H. M. Thomas, A. M. Lipaev, A. D. Usachev, V. I. Molotkov, and M. H. Thoma, Phys. Plasmas 26, 113703 (2019).
  * [15] A. V. Ivlev, M. Kretschmer, M. Zuzic, G. E. Morfill, H. Rothermel, H. M. Thomas, V. E. Fortov, V. I. Molotkov, A. P. Nefedov, A. M. Lipaev, O. F. Petrov, Yu. M. Baturin, A. I. Ivanov, and J. Goree, Phys. Rev. Lett. 90, 055003 (2003).
  * [16] U. Konopka “Super Particle Identification and Tracking” (unpublished).
  * [17] L. D. Landau and E. M. Lifshitz, Fluid Mechanics (Butterworth-Heinemann, Boston, 1997), Vol. 6.
  * [18] Y. Feng, J. Goree, and B. Liu, Phys. Rev. E 86, 056403 (2012).
  * [19] H. Rothermel, T. Hagl, G. E. Morfill, M. H. Thoma, and H. M. Thomas, Phys. Rev. Lett. 89, 175001 (2002).
  * [20] V. Nosenko, F. Luoni, A. Kaouk, M. Rubin-Zuzic, H. Thomas, Phys. Rev. Research 2, 033226 (2020).
  * [21] P. Epstein, Phys. Rev. 23, 710 (1924).
  * [22] M. Y. Pustylnik, B. Klumov, M. Rubin-Zuzic, A. M. Lipaev, V. Nosenko, A. D. Usachev, A. D. Zobnin, V. I. Molotkov, G. Joyce, H. M. Thomas, M. H. Thoma, O. F. Petrov, V. E. Fortov, and O. Kononenko, Phys. Rev. Res. 2, 033314 (2020).
  * [23] M. Schwabe, S. A. Khrapak, S. K. Zhdanov, M. Y. Pustylnik, C. Räth, M. Fink, M. Kretschmer, A. M. Lipaev, V. I. Molotkov, A. S. Schmitz, M. H. Thoma, A. D. Usachev, A. V. Zobnin, G. I. Padalka, V. E. Fortov, O. F. Petrov, and H. M. Thomas, New J. Phys. 22, 083079 (2020).
  * [24] S. Mitic, M. Y. Pustylnik, D. Erdle, A. M. Lipaev, A. D. Usachev, A. V. Zobnin, M. H. Thoma, H. M. Thomas, O. F. Petrov, V. E. Fortov, and O. Kononenko, Phys. Rev. E 103, 063212 (2021).
  * [25] G. E. Morfill, M. Rubin-Zuzic, H. Rothermel, A. V. Ivlev, B. A. Klumov, H. M. Thomas, U. Konopka, V. Steinberg, Phys. Rev. Lett. 92, 175004 (2004).