# Geometric methods for efficient planar swimming of copepod nauplii

Jonas Balisacan  University of Hawai‘i at Mānoa, 2565 McCarthy Mall, Honolulu,
HI 96822, USA<EMAIL_ADDRESS>,  Monique Chyba  University of Hawai‘i at
Mānoa, 2565 McCarthy Mall, Honolulu, HI 96822, USA<EMAIL_ADDRESS>,  Corey
Shanbrom  California State University, Sacramento, 6000 J St., Sacramento, CA
95819, USA<EMAIL_ADDRESS>and  George Wilkens  University of Hawai‘i
at Mānoa, 2565 McCarthy Mall, Honolulu, HI 96822, USA<EMAIL_ADDRESS>
###### Abstract.

Copepod nauplii are larval crustaceans with important ecological functions.
Due to their small size, they experience an environment of low Reynolds number
within their aquatic habitat. Here we provide a mathematical model of a
swimming copepod nauplius with two legs moving in a plane. This model allows
for both rotation and two-dimensional displacement by periodic deformation of
the swimmer’s body. The system is studied from the framework of optimal
control theory, with a simple cost function designed to approximate the
mechanical energy expended by the copepod. We find that this model is
sufficiently realistic to recreate behavior similar to those of observed
copepod nauplii, yet much of the mathematical analysis is tractable. In
particular, we show that the system is controllable, but there exist singular
configurations where the degree of non-holonomy is non-generic. We also
partially characterize the abnormal extremals and provide explicit examples of
families of abnormal curves. Finally, we numerically simulate normal extremals
and observe some interesting and surprising phenomena.

###### Key words and phrases:

Microswimmer; Planar Motion; Maximum Principle; Abnormal Extremals, Elastica

M.C. is partially supported by the Simons Foundation grant number 359510.

##  1\. Introduction

Microcrustaceans known as copepods are one of the most abundant animals on
Earth. They are a type of zooplankton that serve as an important link in the
marine food web. As they are prey for many larger aquatic creatures, they must
adapt strategies that help them maximize their survivability. Observations
have shown that different copepods have adapted different types of movement to
efficiently forage for food and evade predators [9, 10, 20, 28]. Figure 1
shows a nauplius of the copepod species Bestiolina similis.

In the world of microorganisms, water becomes a very viscous fluid in which
movement produces negligible inertia. This is known as a low Reynolds number
environment, as the Reynolds number $\mathcal{R}$ represents the ratio between
the inertial force due to momentum and the viscous force experienced from the
resistance of the liquid. An object swimming in some fluid experiences the
Reynolds number $\mathcal{R}=\frac{av\rho}{\eta}$, where $a$ is the
characteristic dimension of the object, $v$ is the velocity of the object,
$\rho$ is the density of the fluid, and $\eta$ is the fluid viscosity. In
water, organisms as small as bacteria have a Reynolds number of approximately
$10^{-6}$ to $10^{-4}$, whereas humans have a Reynolds number of approximately
$10^{6}$. Nauplii of the paracalanid copepod Bestiolina similis, as shown in
Figure 1, have lengths 70–200 $\mu$m and swim at Reynolds numbers of
$10^{-1}-10^{1}$ (see [16] and the references therein). To put swimming at low
Reynolds number in perspective, humans have a Reynolds number around $10^{2}$
when swimming in molasses.

Figure 1. Scanning electron micrograph of a larval copepod, showing three
pairs of appendages: first antenna (A1), second antenna (A2) and mandible
(Md). Image reproduced from [8] under Creative Commons Attribution 4.0
International License.

Edward Mills Purcell’s talk from 1976, Life at Low Reynolds Number ([22]),
first popularized the concept of swimming in an environment with low Reynolds
number and was foundational in the study of microswimmers. He stated the
scallop theorem, which says that complete reciprocal motion cannot produce any
displacement when swimming at low Reynolds number. Here a reciprocal motion
involves two sequences of motion where the second motion, called the recovery
stroke, is the reverse of the first motion, which is called the power stroke.
For humans swimming in water, the time it takes to do a stroke plays a role in
the induced displacement due to the inertia terms in the Navier-Stokes
equations. So humans are able to swim forward by completing the second
recovery stroke faster than the first. But for microswimmers, this inertia is
negligible, so any displacement induced by the first stroke is reversed by the
second recovery stroke. As a consequence, microorganisms must move in other
ways, like utilizing a flagella or moving pairs of legs in an asynchronous
manner. For example, for copepods that are bilaterally symmetrical, the
symmetric pairs of legs move simultaneously, but adjacent legs move out of
sync in order to produce a net displacement. These observations are important
as they could potentially be used to design microscopic robots that move in
similar ways. One application involves using bacteria based nanoswimmers to
transport drugs from a loading point to a destination such as cancer cells
([13]).

The main microorganisms of focus in this paper are copepods. Most types of
copepod are only able to move in a motion called swimming-by-jumping ([2]).
This motion is similar to the one described in the previous paragraph: it
involves moving symmetric pairs of legs in a way such that adjacent pairs move
asynchronously. In other words, the power stroke and recovery stroke alternate
among each of the pairs of legs. One way to model this is by moving each
symmetric pair of legs in a reciprocal motion while introducing a phase lag
between each pair of legs ([26]). As in such models, here we restrict the
copepod’s motion to a plane for simplicity, despite the fact that the actual
animals live and move in three dimensions.

Here we model the copepod as a slender body in Stokes flow as in [26, 4, 7, 8,
11, 6]. Other models of micro-swimmers capable of rotation appear in [12, 23,
15]. It is important to note that real living copepods do indeed perform
rotations to both evade predators ([24]) and to capture prey ([10]). In [19],
the authors analyze such rotational maneuvers (yaw, pitch and roll) via high-
speed video observations of copepod larvae.

One dimensional translational motion for the copepod model has been well
studied. It is possible to achieve positive displacement along an axis using
as few as two pairs of legs moving in a reciprocal motion ([26]). Methods from
sub-Riemannian geometry and Hamiltonian dynamics have been used to find
efficient optimal strokes in the translational case; efficiency was defined as
the ratio between the displacement resulting from a stroke and the length of
the stroke. Numerical methods were used to determine the optimal strokes
maximizing this efficiency ([4, 7, 8]). Here the term stroke refers to a
periodic motion of the legs.

Here we generalize this prior work by analyzing planar motions. To produce
orientation changes we need to break the symmetry of the pair of legs. A first
attempt was made in [11] by looking at three independent legs oscillating
sinusoidally; here we generalize that approach to include all strokes but
using two legs. We find that rotation by strokes is indeed possible with only
two legs. We also show that the two-legged system is controllable, although
the difficulty in steering locally depends on the initial state; in other
words, the system possesses singularities, which we classify. Taking the
mechanical energy expended as our cost function, we develop the two-legged
copepod movement as an optimal control problem and apply the Pontryagin
maximum principle ([21]) to study both the normal and abnormal extremals. We
partially characterize the abnormal extremals, and provide some explicit
examples. Finally, we utilize the optimal control software Bocop to simulate
normal extremals ([27]). Among our simulations we find copepod motions which
produce net rotation without net displacement, we characterize the optimal
motions which produce rotation with no conditions on displacement, and we
discover paths in the $xy$-plane which appear to be Euler elastica ([3]).

##  2\. Methods

We consider a simplified copepod microswimmer in a low Reynolds number
environment. The idealized copepod consists of stiff slender legs and a body
of negligible radius in comparison to the length of its legs. In this section
we will develop our mathematical model, derive the equations of motion,
describe the copepod’s motion as an optimal control problem, and develop the
appropriate version of the maximum principle.

###  2.1. Model

We assume the copepod moves in a plane and possesses $6$ independently moving
legs, three on each side of the body. The position of the copepod at time $t$
can be described by the vector $\left(x(t),y(t),\phi(t)\right)^{T}$, where $x$
and $y$ represent the usual Cartesian coordinates on the plane and $\phi$
represents the orientation of the copepod with respect to the positive
$x$-axis. Let $\theta_{i}$ denote the angle between the copepod’s orientation
and the $i^{th}$ leg, and let $\alpha_{i}=\theta_{i}+\phi$ denote the angle
between the $i^{th}$ leg and the positive $x$-axis. See Figure 2 for an
illustration.

Figure 2. The rotating copepod with $6$ legs. Note that angles $\theta_{i}$
are associated to the body frame while the angles $\alpha_{i}$ are associated
to the inertial frame .

We denote the state of the copepod at time $t$ by the vector

${q}(t)=\left(x(t),y(t),\phi(t),\theta_{1}(t),\dots\theta_{6}(t)\right)^{T},$

while the position and orientation coordinates alone will be written as
$\hat{{q}}=\left(x(t),y(t),\phi(t)\right)^{T}.$ Thus our configuration space
is ostensibly $\mathbb{R}^{9}$, however in order to prevent the legs from
passing each other, we impose the constraint

$0\leq\theta_{1}\leq\theta_{2}\leq\theta_{3}\leq\pi\leq\theta_{4}\leq\theta_{5}\leq\theta_{6}\leq
2\pi.$ (1)

For the rest of this paper we will focus on a simplified copepod with two
independent legs, one on each side of the body. See Figure 3. This
simplification allows us to conduct a mathematical analysis and is justified
by assuming the three legs on each side of the body are collapsed into one
stronger leg. As will be seen in Section 4, even with this simplification we
obtain swimming motions reflecting actual laboratory observations.

Figure 3. The rotating copepod with $2$ legs.

In most of Section 3, including all of Section 3.2, we use standard techniques
from optimal control ([21]) and sub-Riemannian geometry ([18]). In Section
3.3, however, we utilize the optimal control software Bocop. As stated in
[27], this software approximates our optimal control problem by a finite
dimensional optimization problem using the direct transcription approach to
time discretization. The resulting nonlinear programming problem is solved
using the software package Ipopt, using sparse exact derivatives computed by
ADOL-C.

###  2.2. Equations of Motion

We first develop the equations of motion for $n$ legs, then specify to the
case $n=2$. Our equations consist of a system of differential equations of the
form ${M}\dot{\hat{{q}}}=K$. The equations of motion for a copepod moving in
two dimensions are derived in [11], which focuses on legs moving in an
oscillatory motion: $\theta_{i}(t)=a\cos(t+k_{i})+\beta_{i}$. Parameters are
constrained to ensure that adjacent legs never overlap but possess a phase
lag. The author shows that no net rotation is possible with such a motion for
two legs, thus most of the analysis concerns the case of three legs. For
numerical simulations, the values of $a,k_{1},\beta_{1},\beta_{2},$ and
$\beta_{3}$ are fixed and the total change in orientation and displacement is
computed for varying values of $k_{2}$ and $k_{3}$. The change in displacement
and orientation is maximized when
$(k_{2},k_{3})=(\frac{2\pi}{3},\frac{4\pi}{3})$ and
$(k_{2},k_{3})=(\frac{4\pi}{3},\frac{2\pi}{3})$. In addition, the total work
done by the microswimmer is calculated and a notion of turning efficiency is
introduced.

The system is derived from slender-body Stokes flow, using the fact that, at
low Reynolds number, inertial forces are negligible and the Navier-Stokes
equations can be linearized. Here

$K=\sum_{i=1}^{n}\dot{\theta}_{i}\left(\begin{array}[]{ccc}\sin\alpha_{i}\\\
-\cos\alpha_{i}\\\ -\frac{2}{3}\end{array}\right)$ (2)

and $M$ is the resistance matrix given by

$M=\left(\begin{array}[]{ccc}\sum(1+\sin^{2}\alpha_{i})&-\sum(\sin\alpha_{i}\cos\alpha_{i})&-\sum\sin\alpha_{i}\\\
-\sum(\sin\alpha_{i}\cos\alpha_{i})&\sum(1+\cos^{2}\alpha_{i})&\sum\cos\alpha_{i}\\\
-\sum\sin\alpha_{i}&\sum\cos\alpha_{i}&2\end{array}\right).$ (3)

By computing the mobility matrix $M^{-1}$ (well defined since $M$ is symmetric
and positive definite) we obtain the equations of motion

$\dot{\hat{{q}}}=M^{-1}K.$ (4)

###  2.3. Optimal Control Framework

We now consider this system from the control theoretic point of view, where
the angular velocities of the leg are taken as controls. That is, we set
$u_{i}=\dot{\theta}_{i}$, and assume these are measurable functions of time.
Now let

${K}_{i}=\left(\begin{array}[]{ccc}\sin\alpha_{i}\\\ -\cos\alpha_{i}\\\
-\frac{2}{3}\end{array}\right).$ (5)

Then our control vector fields are

${F}_{i}=\left({M}^{-1}{K}_{i},0,...,0,1,0,...,0\right)^{T},$ (6)

where the $1$ appears in the $i^{th}$ entry after the ${M}^{-1}{K}_{i}$
entries. Then the copepod motion is described by the driftless affine control
system

${\dot{q}}(t)=\sum_{i=1}^{n}u_{i}(t){F}_{i}({q}(t))$ (7)

where $n$ refers to the number of legs. In this work we assume no bounds on
the control. In reality, of course, there are limits to how quickly an actual
copepod can move its legs. But some of this issue is mitigated by the fact
that we will be minimizing some form of energy; see Equation (11) below. We
will see in the next section that, when $n=2$, this is a controllable system;
that is, it is possible to find controls steering the copepod from any given
initial state $q_{\text{initial}}$ to any given terminal state
$q_{\text{final}}$.

he rest of this paper will concern the special case of the copepod with two
legs, so here we record these equations of motion explicitly. When $n=2$, we
have

${M}=\left(\begin{array}[]{ccc}2+\sin^{2}\alpha_{1}+\sin^{2}\alpha_{2}&-\sin\alpha_{1}\cos\alpha_{1}-\sin\alpha_{2}\cos\alpha_{2}&-\sin\alpha_{1}-\sin\alpha_{2}\\\
-\sin\alpha_{1}\cos\alpha_{1}-\sin\alpha_{2}\cos\alpha_{2}&2+\cos^{2}\alpha_{1}+\cos^{2}\alpha_{2}&\cos\alpha_{1}+\cos\alpha_{2}\\\
-\sin\alpha_{1}-\sin\alpha_{2}&\cos\alpha_{1}+\cos\alpha_{2}&2\end{array}\right)$

and

${K}=\dot{\theta}_{1}\left(\begin{array}[]{ccc}\sin\alpha_{1}\\\
-\cos\alpha_{1}\\\
-\frac{2}{3}\end{array}\right)+\dot{\theta}_{2}\left(\begin{array}[]{ccc}\sin\alpha_{2}\\\
-\cos\alpha_{2}\\\
-\frac{2}{3}\end{array}\right)=\dot{\theta}_{1}{K}_{1}+\dot{\theta}_{2}{K}_{2}.$

Our control system then takes the form

${\dot{q}}(t)=u_{1}(t){F}_{1}(q(t))+u_{2}(t){F}_{2}(q(t).$ (8)

where ${q}=(x,y,\phi,\theta_{1},\theta_{2})^{t},$ our controls are
$u_{1}=\dot{\theta}_{1}$ and $u_{2}=\dot{\theta}_{2}$, and the control vector
fields are

${F}_{1}=\left(\begin{array}[]{ccc}{M}^{-1}{K}_{1}\\\ 1\\\
0\end{array}\right)\qquad\text{and}\qquad{F}_{2}=\left(\begin{array}[]{ccc}{M}^{-1}{K}_{2}\\\
0\\\ 1\end{array}\right).$

Straightforward calculations give

${F}_{1}=\left(\begin{array}[]{ccc}-\frac{\sin(\theta_{1}-2\theta_{2}-\phi)-\sin(2\theta_{1}-\theta_{2}+\phi)+17\sin(\theta_{1}+\phi)-7\sin(\theta_{2}+\phi)}{24(\cos(\theta_{1}-\theta_{2})-3)}\\\
-\frac{\cos(\theta_{1}-2\theta_{2}-\phi)+\cos(2\theta_{1}-\theta_{2}+\phi)-17\cos(\theta_{1}+\phi)+7\cos(\theta_{2}+\phi)}{24(\cos(\theta_{1}-\theta_{2})-3)}\\\
\frac{1}{12}(\cos(\theta_{1}-\theta_{2})-3)\\\ 1\\\ 0\end{array}\right),$ (9)

${F}_{2}=\left(\begin{array}[]{ccc}-\frac{\sin(\theta_{1}-2\theta_{2}-\phi)-\sin(2\theta_{1}-\theta_{2}+\phi)+17\sin(\theta_{2}+\phi)-7\sin(\theta_{1}+\phi)}{24(\cos(\theta_{1}-\theta_{2})-3)}\\\
-\frac{\cos(\theta_{1}-2\theta_{2}-\phi)+\cos(2\theta_{1}-\theta_{2}+\phi)-17\cos(\theta_{2}+\phi)+7\cos(\theta_{1}+\phi)}{24(\cos(\theta_{1}-\theta_{2})-3)}\\\
\frac{1}{12}(\cos(\theta_{1}-\theta_{2})-3)\\\ 0\\\ 1\end{array}\right).$ (10)

In the sequel we let $D$ denote the distribution spanned by $F_{1}$ and
$F_{2}$. Note that we can obtain ${F}_{2}$ from ${F}_{1}$ by switching the
roles of $\theta_{1}$ and $\theta_{2}$. Moreover, system (8) is time-
reversible; indeed the transformation $t\mapsto 2\pi-t,u_{i}(t)\mapsto-
u_{i}(2\pi-t)$ sends $q(t)$ to $q(2\pi-t)$. This is a general consequence of
swimming at low Reynolds number. The system is also invariant under rigid body
transformations. It is obvious for translation as the $F_{i}$ do not depend on
$x$ or $y$, but it can also be verified that $\phi\mapsto\phi+\tau$ is
invariant under $x\mapsto x\cos\tau-y\sin\tau,y\mapsto x\sin\tau+y\cos\tau$.

In this paper we suppose that the copepod seeks to minimize the mechanical
energy expended when moving from one position to another. In [8], the authors
consider a two-legged copepod moving along an axis without rotation. They
describe a realistic but complicated mechanical energy functional, but show
that the resulting optimal trajectories are qualitatively similar to those
obtained when using the simplified energy

$E(u_{1},u_{2})=\int_{t_{0}}^{t_{f}}\left(u_{1}^{2}+u_{2}^{2}\right)\,dt.$
(11)

Here $t_{0}$ is a fixed initial time, while $t_{f}$ is associated to the
control $u$ in the following manner. Choose some terminal boundary manifold
$M_{1}\subseteq\mathbb{R}^{5}$, which is closed, and define the target set by
$M=[t_{0},\infty)\times M_{1}$. Then $t_{f}$ is the smallest time such that
$(t_{f},q(t_{f}))\in M$, where $q(t)$ is the state trajectory associated to
the control $u(t)$.

We therefore choose our cost function to be the energy $E$ given in (11)
corresponding to the orthonormal inner product for the two control vector
fields, yielding the following optimal control formulation. Provided certain
boundary conditions made explicit below, we seek solutions to the dynamical
system (8) which minimize the cost (11). Note that this is a sub-Riemannian
problem associated to the flat metric ([5]).

###  2.4. Maximum Principle

The Pontryagin maximum principle provides necessary conditions for a solution
to be optimal. The general statement can be found in the literature [5, 1, 17,
21], and we here only state it for our application.

Consider the optimal control problem stated in the previous section, defined
by the dynamics (8) and cost (11). Let $j,k\leq 5$ and consider the initial
and terminal boundary manifolds

$\displaystyle M_{0}$ $\displaystyle=\\{q\in\mathbb{R}^{5}\colon
g_{1}(q)=g_{2}(q)=\cdots=g_{j}(q)=0\\}$ $\displaystyle M_{1}$
$\displaystyle=\\{q\in\mathbb{R}^{5}\colon
h_{1}(q)=h_{2}(q)=\cdots=h_{k}(q)=0\\}.$

Define the Hamiltonian by

$H(q,u,p,p_{0})=\langle
p,u_{1}F_{1}+u_{2}F_{2}\rangle+p_{0}(u_{1}^{2}+u_{2}^{2}).$ (12)

######  Theorem 2.1 (Maximum Principle).

Let $u^{*}\colon[t_{0},t_{f}]\to\mathbb{R}^{2}$ be an optimal control and let
$q^{*}\colon[t_{0},t_{f}]\to\mathbb{R}^{5}$ be the corresponding optimal state
trajectory. Then there exists a function
$p^{*}\colon[t_{0},t_{f}]\to\mathbb{R}^{5}$ and a constant $p_{0}^{*}\leq 0$
such that $(p_{0}^{*},\,p^{*}(t))\neq(0,0)$ for all $t\in[t_{0},t_{f}]$ and
having the following properties:

  1. (1)

$q^{*}$ and $p^{*}$ satisfy Hamilton’s equations for the Hamiltonian (12) with
boundary conditions $q^{*}(t_{0})\in M_{0}$ and $q^{*}(t_{f})\in M_{1}$.

  2. (2)

$H(q^{*}(t),u^{*}(t),p^{*}(t),p^{*}_{0})\geq H(q^{*}(t),u,p^{*}(t),p^{*}_{0})$
for all $t\in[t_{0},t_{f}]$ and $u\in\mathbb{R}^{2}$.

  3. (3)

$H(q^{*}(t),u^{*}(t),p^{*}(t),p^{*}_{0})=0$ for all $t\in[t_{0},t_{f}]$.

  4. (4)

The vector $p^{*}(0)$ is orthogonal to the tangent space
$T_{q^{*}(t_{0})}M_{0}$ and the vector $p^{*}(t_{f})$ is orthogonal to the
tangent space $T_{q^{*}(t_{f})}M_{1}$.

Note that this version of the maximum principle closely follows Section 4.1.2
of [17], with the addition of the initial boundary conditions $M_{0}$ and the
associated transversality conditions. It important, however, to recognize that
Theorem 2.1 pertains to the free-time problem. In our simulations in Section
3.2 we fix the time $[t_{0},t_{f}]=[0,2\pi]$, which requires a minor
modification of the maximum principle. As described in Section 4.3.1 of [17],
one simply introduces an extra state variable to represent time, $q_{6}=t$,
and includes the fixed terminal time in the terminal manifold $M_{1}$. Then
the free-time Theorem 2.1 applies with only the following modification:
$H|_{*}=-p_{6}^{*}$, which is constant. In the sequel we will refer to a
trajectory satisfying the conclusions of the maximum principle as an extremal.

##  3\. Results

Most of the results in this section concern the two-legged copepod. Unless
explicitly stated otherwise, we assume $n=2$. The motivation for this choice
is contained in the next theorem, which says that one leg is insufficient for
producing rotation via periodic strokes, but two legs are sufficient. As in
[6], we define a _stroke_ to be a periodic deformation of the swimmer’s body.
That is, a stroke of period $T$ is any path in configuration space satisfying
$\theta_{i}(0)=\theta_{i}(T)$ for all $i=1,\dots n$. For simulations and
examples we also impose the following realistic constraint forcing each leg to
stay on one side of the copepod’s body:

$0\leq\theta_{1}\leq\pi\leq\theta_{2}\leq 2\pi.$ (13)

However, most of the mathematical analysis in this section is valid for the
configuration space $\mathbb{R}^{2}\times(S^{1})^{3}$. Note that [11] proves
that a two-legged copepod is incapable of producing net rotation via the
specific oscillating strokes considered in their work.

######  Theorem 3.1.

A one-legged copepod moving in strokes can neither produce a net rotation nor
net displacement. A two-legged copepod moving in strokes can produce net
rotation.

###### Proof.

For a copepod with one leg, we compute

$M=\left(\begin{array}[]{ccc}1+\sin^{2}\alpha&-\sin\alpha\cos\alpha&-\sin\alpha\\\
-\sin\alpha\cos\alpha&1+\cos^{2}\alpha&\cos\alpha\\\
-\sin\alpha&\cos\alpha&2\end{array}\right)\qquad\text{and}\qquad
K=\dot{\theta}\left(\begin{array}[]{ccc}\sin\alpha\\\ -\cos\alpha\\\
-\frac{2}{3}\end{array}\right).$

Then the equation $\dot{\hat{{q}}}=M^{-1}K$ simplifies to

$\left(\begin{array}[]{ccc}\dot{x}\\\ \dot{y}\\\
\dot{\phi}\end{array}\right)=\dot{\theta}\left(\begin{array}[]{ccc}\frac{4}{9}\sin\alpha\\\
-\frac{4}{9}\cos\alpha\\\ -\frac{1}{9}\end{array}\right).$

Suppose the copepod moves in strokes of period $T$, so $\theta(0)=\theta(T)$.
Without loss of generality assume $\phi(0)=0.$ Then
$\phi(t)=-\frac{1}{9}\theta(t)+\frac{1}{9}\theta(0)$. Because $\theta(t)$ is
periodic, $\phi(t)$ is also periodic and thus no net rotation is produced
after one period.

Now we also have that
$\alpha(t)=\theta(t)+\phi(t)=\frac{8}{9}\theta(t)+\frac{1}{9}\theta(0)$. To
calculate the total $x$ displacement over the period we compute the integral

$\int_{0}^{T}\dot{x}(t)\,dt=\int_{0}^{T}\frac{4}{9}\sin\left(\frac{8}{9}\theta(t)+\frac{1}{9}\theta(0)\right)\dot{\theta}(t)\,dt.$

Computing this integral gives us

$x(T)-x(0)=\left.-\frac{1}{2}\left(\cos\left(\frac{8}{9}\theta(t)+\frac{1}{9}\theta(0)\right)\right)\right|_{0}^{T},$

which equals zero again since $\theta$ is periodic. Thus the displacement in
the $x$ direction is 0; a similar argument gives the result for $y$.

The proof of the second statement in the theorem is contained in the following
example. ∎

The next example demonstrates that two legs moving in strokes can produce a
net change in displacement and orientation.

######  Example 3.2.

Consider the two-legged copepod with initial configuration

${q}(0)=(x(0),y(0),\phi(0),\theta_{1}(0),\theta_{2}(0))^{T}=(0,0,0,0,\pi)^{T}.$

We consider the following motion. Both legs rotate $\pi$ radians counter-
clockwise in $\pi$ time; then one at a time, each leg moves $\pi$ radians
clockwise in $\frac{\pi}{2}$ time. See Figure 4. Explicitly, we have

$\theta_{1}(t)=\begin{cases}t&0\leq t<\pi\\\ -2t+3\pi&\pi\leq
t<\frac{3\pi}{2}\\\ 0&\frac{3\pi}{2}\leq
t<2\pi\end{cases},\indent\theta_{2}(t)=\begin{cases}t+\pi&0\leq t<\pi\\\
2\pi&\pi\leq t<\frac{3\pi}{2}\\\ -2t+5\pi&\frac{3\pi}{2}\leq
t<2\pi.\end{cases}$

We can explicitly solve the equations of motion for the orientation over time,

$\phi(t)=\begin{cases}-\frac{2}{3}t&0\leq t<\pi\\\
\frac{1}{12}\sin\left(2t\right)+\frac{1}{2}t-\frac{7\pi}{6}&\pi\leq t\leq
2\pi,\\\ \end{cases}$

which implies a net rotation of $\phi(2\pi)-\phi(0)=-\frac{\pi}{6}\neq 0.$ The
orientation over time, along with the displacements over time obtained by
numerical integration, are shown in Figure 5. The net change in position,
equal to the final position, is given by

$\hat{{q}}\left(2\pi)=(x(2\pi),y(2\pi),\phi(2\pi)\right)^{T}=(0.0071,0.0019,-{\pi}/{6})^{T}.$

Figure 4. The motion described in Example 3.2 traces out the pink triangle in
the $\theta_{1}\theta_{2}$-plane, counter-clockwise from the bottom left
vertex. The black box represents the constraint (13).

Figure 5. Motion in Example 3.2; the three colors correspond to the three legs
of the triangle in Figure 4. (a) Orientation $\phi(t)$. (b) Displacements
$x(t)$ and $y(t)$.

Note that the energy (11) for the motion in Example 3.2 is equal to $6\pi$.
While this motion is dynamically valid, it is likely not minimizing the cost.
See Section 3.3 for further discussion.

###  3.1. Controllability

In this subsection we show that the two-legged copepod is indeed a
controllable system. That is, given any initial and final configuration,
controls exist which steer the copepod from the initial to the final
configuration. The main tool here is the Chow-Rashevskii theorem; a formal
statement, along with definitions of all the terminology in this subsection,
can be found in [18]. The proofs here are essentially just calculations, which
we performed using MATLAB and Mathematica.

Our two control vector fields $F_{1}$ and $F_{2}$ are given by (9) and (10).
Denote their iterated Lie brackets (which are too complicated to display) by

$F_{3}=[F_{1},F_{2}],\quad F_{4}=[F_{1},F_{3}],\quad F_{5}=[F_{2},F_{3}].$

Note that $F_{1}$ and $F_{2}$, and consequently their iterated brackets, only
depend on $\theta_{1},\theta_{2}$, and $\phi$. A computation shows that the
five vector fields $F_{1},F_{2},F_{3},F_{4},$ and $F_{5}$ are linearly
dependent if and only if

${\sin^{4}\left(\frac{\psi}{2}\right)\left(25\cos(2\psi)+120\cos(\psi)+79\right)}=0$
(14)

where $\psi=\theta_{1}-\theta_{2}$. Note that this is only a condition on two
of our five variables. Solving (14) yields two sets of solutions: $\psi=2\pi
n$ and $\psi=2\pi n\pm 2\arctan(2)$ for $n\in\mathbb{Z}$. We therefore have
two sets of configurations which are singular for our distribution $D$:

$\displaystyle S_{1}$ $\displaystyle=\\{{q}\,|\,\theta_{1}-\theta_{2}=2\pi
n\pm 2\arctan(2)\,\text{for some}\,n\in\mathbb{Z}\\}$ $\displaystyle S_{2}$
$\displaystyle=\\{{q}\,|\,\theta_{1}-\theta_{2}=2\pi n\,\text{for
some}\,n\in\mathbb{Z}\\}.$

These are illustrated in the $\theta_{1}\theta_{2}$-plane in Figure 6.

Figure 6. The singular sets for $D$. The green lines are given by $\psi=2\pi
n$, the red lines by $\psi=2\pi n-2\arctan(2)$, and the blue lines by $2\pi
n+2\arctan(2)$, where $n\in\mathbb{Z}$ and $\psi=\theta_{1}-\theta_{2}$. The
black box again shows the constraint (13).

######  Theorem 3.3.

If ${q}\in S_{1}$ then the small growth vector at ${q}$ is (2, 3, 4, 5). If
${q}\in S_{2}$ then the small growth vector at ${q}$ is (2, 2, 3, 4, 5). All
other points are regular with small growth vector (2, 3, 5).

###### Proof.

The fact that the small growth vector is (2, 3, 5) at generic points is
immediate from the fact that $F_{1},F_{2},F_{3},F_{4},$ and $F_{5}$ are
linearly independent there. These points are regular since the singular sets
are closed. Points in the singular sets are analyzed by computations. Let

$F_{6}=[F_{1},F_{4}],\quad F_{7}=[F_{1},F_{5}],\quad F_{8}=[F_{2},F_{4}],\quad
F_{9}=[F_{2},F_{5}],\quad F_{10}=[F_{1},F_{6}].$

For points in $S_{1}$, we find that $F_{4}=-F_{5}$, but
$F_{1},F_{2},F_{3},F_{4},$ and $F_{6}$ are linearly independent. For points in
$S_{2}$, we find that $F_{3}=F_{7}=F_{8}=0$ and $F_{4}=-F_{5}$ and
$F_{6}=-F_{9}$, but $F_{1},F_{2},F_{4},F_{6},$ and $F_{10}$ are linearly
independent. ∎

######  Corollary 3.4.

The two-legged copepod system is controllable at all points. At generic points
the degree of non-holonomy is 3. At points in $S_{1}$ the degree of non-
holonomy is 4. At points in $S_{2}$ the degree of non-holonomy is 5.

###### Proof.

Controllability follows from the Chow-Rashevskii theorem, as the vector fields
$F_{1}$ and $F_{2}$ Lie generate the tangent bundle at every point. The degree
of non-holonomy is simply the length of the small growth vector. ∎

###  3.2. Abnormal extremals.

Abnormal extremals are intrinsic to the dynamics; they do not depend on the
cost. It is well known that they play a very important role for the optimal
synthesis ([5]). They correspond to imposing $p_{0}=0$ in the Hamiltonian
(12). It follows that, for our application, the abnormal Hamiltonian is

$H_{a}(q,p,u)=u_{1}\langle p,F_{1}(q)\rangle+u_{2}\langle p,F_{2}(q)\rangle.$

According to the Pontryagin maximum principle, abnormal extremals are curves
$(q(t),p(t))$ satisfying Hamilton’s equations for $H_{a}$ as well as

$\displaystyle\langle p,F_{1}(q)\rangle$ $\displaystyle=0$ (15)
$\displaystyle\langle p,F_{2}(q)\rangle$ $\displaystyle=0.$ (16)

Differentiating these equations leads to the additional requirements

$\displaystyle\langle p,F_{3}(q)\rangle$ $\displaystyle=0$ (17) $\displaystyle
u_{1}\langle p,F_{4}(q)\rangle+u_{2}\langle p,F_{5}(q)\rangle$
$\displaystyle=0.$ (18)

The following results partially characterize the abnormal curves.

######  Proposition 3.5.

The horizontal lifts of the singular curves in Figure 6 are projections of
abnormal curves for the two-legged copepod. They are the integral curves for
the vector field $F_{1}+F_{2}$ restricted to the singular set $S=S_{1}\cup
S_{2}$ for the distribution $D$, and they project to uniform circular motion
in the $xy$-plane. More precisely, let $q_{0}=q^{i}(0)$, and take any
$n\in\mathbb{Z}$ and any $c_{1},c_{2}\in\mathbb{R}$ not both zero. Then the
curves are $(q^{1}(t),p^{1}(t))$ and $(q^{2}(t),p^{2}(t))$ are abnormal, where

$\displaystyle q^{1}(t)$
$\displaystyle=\begin{pmatrix}-\frac{\sqrt{5}}{6}\cos(\frac{2}{5}t\pm\arctan
2+\phi_{0})+c_{x}\\\ -\frac{\sqrt{5}}{6}\sin(\frac{2}{5}t\pm\arctan
2+\phi_{0})+c_{y}\\\ -\frac{3}{5}t+\phi_{0}\\\ t\\\ t\pm 2\arctan 2-2\pi
n\end{pmatrix},$ $\displaystyle p^{1}(t)$
$\displaystyle=\begin{pmatrix}c_{1}\\\ c_{2}\\\
\frac{1}{6}[(-2c_{1}+c_{2})\cos(\frac{2}{5}t+\phi_{0})-(c_{1}+2c_{2})\sin(\frac{2}{5}t+\phi_{0})]\\\
\frac{1}{36}[(-2c_{1}+11c_{2})\cos(\frac{2}{5}t+\phi_{0})-(11c_{1}+2c_{2})\sin(\frac{2}{5}t+\phi_{0})]\\\
\frac{1}{36}[(-10c_{1}-5c_{2})\cos(\frac{2}{5}t+\phi_{0})+(5c_{1}-102c_{2})\sin(\frac{2}{5}t+\phi_{0})]\end{pmatrix},$
$\displaystyle q^{2}(t)$
$\displaystyle=\begin{pmatrix}-\frac{1}{2}\cos(\frac{2}{3}t+\phi_{0})+x_{0}+\frac{1}{2}\cos(\phi_{0})\\\
-\frac{1}{2}\sin(\frac{2}{3}t+\phi_{0})+y_{0}+\frac{1}{2}\sin(\phi_{0})\\\
-\frac{t}{3}+\phi_{0}\\\ t\\\ t-2\pi
n\end{pmatrix},\,p^{2}(t)=\begin{pmatrix}0\\\ 0\\\ 6\\\ 1\\\
1\end{pmatrix},\,$

and

$c_{x}=x_{0}+\frac{\sqrt{5}}{6}\cos(\pm\arctan
2+\phi_{0})\qquad\text{and}\qquad
c_{y}=y_{0}+\frac{\sqrt{5}}{6}\sin(\pm\arctan 2+\phi_{0}).$

###### Proof.

Note that $q^{i}(t)$ is just the horizontal lift of the naive parametrization
of the lines which constitute the connected components of $S_{i}$ projected to
the $\theta_{1}\theta_{2}$-plane; these are the colored lines in Figure 6.
Thus for any time $t$ we have $q^{i}(t)\in S_{i}$. The curves $q^{i}$ are
integral curves for $F_{1}+F_{2}$ and therefore horizontal. In fact,
$\text{span}\\{F_{1}+F_{2}\\}$ is the intersection of $D$ and the set of
Cauchy characteristics for $D+[D,D]$. It is also the intersection of $D$ with
the tangent bundle $TS$.

It is straightforward to check that the pairs $(q^{1}(t),p^{1}(t))$ and
$(q^{2}(t),p^{2}(t))$ satisfy Hamilton’s equations. Observe that $p^{i}$ is
constant in the first two components since our control vector fields do not
depend on $x$ or $y$. We also have by construction that

$\dot{q}^{i}=F_{1}+F_{2}|_{q^{i}}=\frac{\partial H_{a}}{\partial
p}\bigg{|}_{(q^{i},p^{i})}.$

Interestingly, the abnormal Hamiltonian also satisfies

$\frac{\partial H_{a}}{\partial q}\bigg{|}_{(q^{2},p^{2})}=0.$

Note that $p^{1}$ is actually a 2-parameter family of curves, and $p^{2}$ is
constant and therefore an integral of motion.

It is similarly straightforward to check that $(q^{i}(t),p^{i}(t))$ satisfy
abnormal equations (15 – 18). Recall that $F_{3}|_{S_{2}}=0$ so equation (17)
is satisfied for any $p$ on $S_{2}$. Similarly, on any continuous curve within
$S$ we have $F_{4}=-F_{5}$ and $u_{1}=u_{2}$, so equation (18) is satisfied
for any $p$. ∎

######  Corollary 3.6.

Abnormal strokes contained in the singular set for $D$ can produce neither
rotation nor displacement.

###### Proof.

The only possible abnormal stroke lying in $S$ would require the legs tracing
a segment of one of the colored lines in Figure 6 first forward and then
backward. The symmetry of such a stroke prevents any net rotation or
displacement. ∎

######  Proposition 3.7.

Neither control can be zero along an abnormal extremal.

###### Proof.

Without loss of generality, assume $u_{2}=0$ along an abnormal extremal
$(q(t),p(t))$, so $\theta_{2}$ is constant. Assume $u_{1}\neq 0$; otherwise
the system is at a stationary point. Equation (18) then implies that $\langle
p,F_{4}(q)\rangle=0$. Differentiating this equation yields

$u_{1}\langle p,F_{6}(q)\rangle+u_{2}\langle p,F_{7}(q)\rangle=0,$ (19)

which reduces to $\langle p,F_{6}(q)\rangle$ since $u_{2}=0$. A
straightforward calculation shows that the vector fields
$F_{1},F_{3},F_{4},F_{6}$ are linearly independent along such curve. Moreover,
they all are identically zero in the fifth component, so the only $p$ mutually
orthogonal to these four vector fields is of the form $(0,0,0,0,p_{5})$. Then
using (10) and $\langle p,F_{2}\rangle=0$ we also have $p_{5}=0$. Thus $p=0$,
which contradicts the maximum principle. ∎

######  Corollary 3.8.

The motion in Example 3.2 is not abnormal. The boundary box (13) is not
abnormal.

######  Proposition 3.9.

The abnormal curves not lying in $S$ are solutions to the Hamiltonian system

$\dot{q}=\frac{\partial\tilde{H}_{a}}{\partial
p},\qquad\dot{p}=-\frac{\partial\tilde{H}_{a}}{\partial q}$

for the Hamiltonian

$\tilde{H}_{a}(q,p)=\langle p,F_{1}\rangle\langle p,F_{5}\rangle-\langle
p,F_{2}\rangle\langle p,F_{4}\rangle.$

###### Proof.

Note that Proposition 3.5 considers the abnormal curves within $S$. Assume
$(q(t),p(t))$ is an abnormal curve not lying within $S$ for any time interval.
By equations (15 – 17), we have that $p$ is orthogonal to
$F_{1}(q),F_{2}(q),F_{3}(q)$. But on the complement of $S$, we have that
$F_{1},F_{2},F_{3},F_{4},F_{5}$ are linearly independent, so $p$ cannot also
be orthogonal to both $F_{4}$ and $F_{5}$. Without loss of generality, assume
$p$ is not orthogonal to $F_{4}$. Then we can solve equation (18) for

$u_{1}=-u_{2}\frac{\langle p,F_{5}\rangle}{\langle p,F_{4}\rangle}$

obtaining

$H_{a}=-u_{2}\frac{\langle p,F_{5}\rangle}{\langle p,F_{4}\rangle}\langle
p,F_{1}\rangle+u_{2}\langle p,F_{2}\rangle.$

Thus we scale $H_{a}$ to obtain a new Hamiltonian in which the controls do not
appear at all:

$\tilde{H}_{a}(q,p)=\langle p,F_{1}\rangle\langle p,F_{5}\rangle-\langle
p,F_{2}\rangle\langle p,F_{4}\rangle.$

As the abnormal equations have been satisfied by construction, any solution to
Hamilton’s equations for this Hamiltonian will indeed be an abnormal curve. ∎

###  3.3. Normal Extremals

Taking $p_{0}=-\frac{1}{2}$, our normal Hamiltonian is

$H_{n}(q,p,u)=u_{1}\langle p,F_{1}(q)\rangle+u_{2}\langle
p,F_{2}(q)\rangle-\frac{1}{2}(u_{1}^{2}+u_{2}^{2}).$

We analyze the normal extremals indirectly, using the optimal control software
Bocop (see Section 2 and [27]). Our investigations have led to two interesting
observations regarding the behavior of normal extremals.

Figure 7. A simulated optimal normal trajectory showing elastica-like motion
in the plane. Here we set the following boundary conditions: $(x,y)(0)=(0,0)$,
$(x,y)(2\pi)=(0.1,0.5)$, and
$(\phi,\theta_{1},\theta_{2})(0)=(\phi,\theta_{1},\theta_{2})(2\pi)$. The
energy here is approximately 54.311. (a) The elastica-like path of the
microswimmer in the $xy$-plane. (b) Orientation $\phi(t)$. (c) The trajectory
in the $\theta_{1}\theta_{2}$-plane, within the constraint square (13). (d)
Angles of the legs $\theta_{1}(t)$ and $\theta_{2}(t)$.

First, certain trajectories seem to show the copepod moving along a curve in
the $xy$-plane which is a type of Euler elastica. See Figure 7 for one
example. In particular, we observe this behavior when fixing the start and end
positions in the plane, demanding that the net change in orientation is zero,
and demanding that the copepod completes a stroke, with no other imposed
boundary conditions. It can be observed that the legs follow a periodic
motion, and in turn the orientation of the copepod is periodic as well. The
motion in the angular phase plane $(\theta_{1},\theta_{2})$ is a perfect
ellipsoid within the constraint space, reflecting the symmetry of the motion
of the two legs.

A possible route to proving that this phenomena holds is suggested by [3]. Our
optimal control problem can be translated into a geodesic problem in sub-
Riemannian geometry. Our control vector fields $F_{1}$ and $F_{2}$ have dual
momenta $P_{1}=\langle p,F_{1}\rangle$ and $P_{2}=\langle p,F_{2}\rangle$, and
the sub-Riemannian Hamiltonian $H_{sR}=\frac{1}{2}(P_{1}^{2}+P_{2}^{2})$
generates normal geodesics corresponding to our normal optimal copepod
trajectories. These geodesics parametrized by arc length correspond to
solutions of Hamilton’s equations for $H_{sR}$ (geodesic equations) with
energy $H=1/2$. These equations could potentially allow us to show that the
curvature $\kappa$ of the projection $(x(t),y(t))$ satisfies one of the
defining differential equations for Euler elastica. The obstacles here are
that the computations are unwieldy, and it is not clear how to impose the
boundary conditions which seem to lead to elastica-like behavior in the
copepod.

Our second interesting observation concerns the triangle $\mathcal{T}$,
appearing in the lower right corner of the constraint square (13), consisting
of the boundary of the set $\\{(\theta_{1},\theta_{2}):\
0\leq\theta_{1}\leq\theta_{2}-\pi\leq\pi\\}$. Our simulations show that
following this triangle is optimal for a copepod desiring to rotate a
prescribed amount. More precisely, suppose we specify the net rotation
$\Delta\phi$ but impose no other boundary conditions: we do not specify the
start or end points in the plane, or require the motion be a stroke. Then the
optimal motion of the legs traces out the triangle $\mathcal{T}$ from the top
right corner counterclockwise; it may go around $\mathcal{T}$ more than once,
not necessarily an integer number of times. In fact, it will never go around
an integer number of times (which would constitute a stroke). Our observations
suggest the following characterization of the motion:

$\Delta\phi\in\begin{cases}[0,\frac{2\pi}{3}]\qquad\text{just follow
hypotenuse: 0 to .5 times around}\ \mathcal{T}\\\
(\frac{2\pi}{3},\frac{5\pi}{6}]\quad\,\text{once around $\mathcal{T}$, then
hypotenuse: 1 to 1.5 times around}\ \mathcal{T}\\\ (\frac{5\pi}{6},\pi]\quad\
\ \text{twice around $\mathcal{T}$, then hypotenuse: 2 to 2.5 times around}\
\mathcal{T}.\\\ \end{cases}$

In any of these cases, the hypotenuse need not be traced out completely. For
example, to rotate $\pi/3$ radians one would simply traverse half the
hypotenuse then stop. See Figure 8 for an example of the third case with
$\Delta\phi=\pi$. Note the symmetry of the legs in Figure 8(d), which is
implicit in the triangle $\mathcal{T}$ itself. Of course, to rotate more than
$\pi$ radians one simply reverses this process (starting at the bottom left
and following $\mathcal{T}$ clockwise).

Figure 8. A simulated normal trajectory showing that following the triangle
$\mathcal{T}$ is optimal for creating rotation. Our only boundary conditions
concern the orientation angle: $\phi(0)=0$ and $\phi(2\pi)=\pi$. The energy
here is approximately 107.735. (a) The path of the animal in the $xy$-plane.
(b) Orientation $\phi(t)$. (c) The trajectory in the
$\theta_{1}\theta_{2}$-plane, within the constraint square (13). The path
traces out the triangle $\mathcal{T}$ exactly 2.5 times counterclockwise
starting from the upper right vertex. (d) Angles of the legs $\theta_{1}(t)$
and $\theta_{2}(t)$.

Note that traveling along the hypotenuse induces no displacement and
traversing the complete triangle induces very small net displacement (see
Example 3.2). Thus these motions represent optimal swimming for a copepod
attempting to rotate any amount without much net displacement. Any rotation
amount less than or equal to $2\pi/3$ can be achieved optimally with no
displacement at all. Intuitively, this demonstrates the fact that traveling
along the hypotenuse gives the strongest possible power stroke for inducing
rotation – the legs of the triangle are simply necessary to move the copepod
legs back into position for another power stroke in a way that minimizes
backwards rotation. A motion which includes the legs of the triangle, as in
Figure 8, does require the copepod to move around in the plane, but it returns
to nearly its original position.

In Figure 9 we provide a catalog of the type of topological curves in the
$xy$-plane obtained when varying the boundary conditions. The boundary
conditions themselves appear in Table 1.

Figure 9. Gallery of simulated normal extremals showing copepod paths in the $xy$-plane. In all cases the copepod begins at the origin $(0,0)$ and $t\in[0,2\pi]$. Additionally, we specify boundary conditions for $x(2\pi),\,y(2\pi),\,\phi(0)$, and $\phi(2\pi)$; these are given in Table 1. No other boundary conditions are imposed.  Table 1. The boundary conditions corresponding to the paths in Figure 9. In all cases the initial position is $(0,0)$. In some simulations certain conditions were unspecified, labeled “Free” here. For some simulations we impose the stroke condition that $\theta_{i}(0)=\theta_{i}(2\pi)$. The total energy of each path is not imposed, but listed here for comparison.  Subfigure | Final position | Initial orientation | Final orientation | Stroke | Energy
---|---|---|---|---|---
(a) | $(0.1,0.1)$ | Free | Free | Yes | 3.868
(b) | $(0.1,0.1)$ | Free | $\phi(0)$ | Yes | 4.670
(c) | $(0.8,0.8)$ | Free | $\phi(0)$ | Yes | 215.560
(d) | $(0.8,0.8)$ | Free | $\phi(0)$ | No | 1.895
(e) | $(0.8,0.8)$ | 0 | $0$ | No | 19.968
(f) | $(0.5,0)$ | 0 | 0 | Yes | 122.723
(g) | $(0.5,0)$ | Free | Free | Yes | 27.745
(h) | $(0,0.5)$ | Free | $\phi(0)$ | Yes | 43.683
(i) | $(0,0.5)$ | Free | $\phi(0)$ | No | 0.319
(j) | $(0,0.5)$ | Free | Free | No | 0.319
(k) | $(0.9,0.3)$ | 0 | 0 | Yes | 197.796
(l) | $(0.9,0.1)$ | Free | $\phi(0)$ | Yes | 146.854
(m) | $(0.1,0.5)$ | 0 | Free | Yes | 30.321
(n) | $(0.7,0.1)$ | 0 | Free | Yes | 80.652
(o) | $(0.1,0.9)$ | Free | 0 | Yes | 91.949

##  4\. Discussion and Conclusions

Here we have provided a mathematical model of a swimming copepod nauplius with
two legs moving in a plane. This model allows for both rotation and two-
dimensional displacement by periodic deformation of the swimmer’s body. The
system was studied from the framework of optimal control theory, with a simple
cost function designed to approximate the mechanical energy expended by the
copepod. We have found that this model is sufficiently realistic to recreate
behavior similar to those of observed copepod nauplii, yet much of the
mathematical analysis is tractable. In particular, we have shown that the
system is controllable, but there exist singular configurations where the
degree of non-holonomy is non-generic. We have also partially characterized
the abnormal extremals and provided explicit examples of families of abnormal
curves. Finally, we have numerically simulated normal extremals and observed
some interesting and surprising phenomena.

This work suggests a plethora of interesting open problems and directions for
future research. First, there are a number of potential generalizations and
modifications to our model which may lead to even more realistic behaviors.
For example, one can study the model with four or six legs coupled with the
appropriate constraints. Real copepods have six legs. Further, this model has
the potential to design soft small-scale synthetic robots [8, 14, 25].
Alternatively, or perhaps additionally, one could work in a three-dimensional
environment, which is obviously more realistic. Moreover, there are other
reasonable cost functions to consider, including more complicated versions of
mechanical energy. Instead, it may be that copepods seek to minimize the time
needed to perform a given motion, or the total distance traveled, either to
evade predators or capture prey more effectively.

Without generalization, our current model already offers ideas for future
research. In particular, it would be quite interesting to find a mathematical,
physical, or biological explanation for the observed elastica-like paths, such
as the one shown in Figure 7. Even numerical verification that these paths are
indeed forms of elastica would be worth pursuing. A potential approach is
described in Section 3.3, which in turn leads to other questions of a
differential geometric flavor. Our optimal control problem can indeed be cast
as the geodesic problem for a particular sub-Riemannian geometry, which
appears geometrically interesting. The vast sub-Riemannian literature may
yield geometric or metric tools providing deeper insight into the copepod
system.

While the experimental approach in Section 3.3 led to some interesting
observations, the normal extremals are still largely not understood. It would
be particularly interesting to explore path-planning for the copepod system.
It is also important to recognize that Bocop, like any mathematical software,
has limitations, some of which we encountered. In particular, this is local
optimization software, and we are not working on a convex optimization problem
with one global extremum. Thus, despite the symmetries of the problem, the
numerical results were sensitive to transversality conditions (for example,
specifying that copepod start at the origin).

Finally, we consider how well our results approximate actual observations of
copepods in motion. In [9] the authors discuss copepod swimming and escape
behavior, based on observations of their swimming patterns and activity. In
particular, Figure 4 in that paper shows a helical pattern projecting onto the
$xy$-plane like an ellipse. This correlates with our motion presented in
Figure 8. Figure 8 in [9] depicts escape trajectories for nauplii and
copepodis which follow helical patterns that project on the $xy$-plane as
Euler elastica. In [10], the authors observe the positions of the appendages
during prey capture and prey handling; in their Figure 5 we see the leg
motions are oscillatory and mostly periodic, as in our simulations, during the
prey handling phase. The most striking comparison comes with the observed
behavior provided in [20]. Indeed, the projection of the 3D swimming motion of
the nauplii and early copepodid in their observations provide a similar
complexity to our simulated trajectories. Compare our Figure 9 to Figures 1
through 9 in their paper. It is quite remarkable that despite the simplified
assumptions made on the number of legs and the cost, our results still capture
the essence of swimming behavior for copepods.

## References

  * [1] A. Agrachev, D. Barilari, U. Boscain, A Comprehensive Introduction to sub-Riemannian Geometry, Cambridge University Press, 2019.
  * [2] C.M. Andersen-Borg , E. Bruno, T. Kiørboe, The kinematics of swimming and relocation jumps in copepod nauplii, PLoS ONE, 7 (2012), e47486.
  * [3] A. Ardentov, G. Bor, E. Le Donne, R. Montgomery, Y. Sachkov, Bicycle paths, elasticae and sub-Riemannian geometry, arXiv:2010.04201, (2021).
  * [4] P. Bettiol, B. Bonnard, J. Rouot, Optimal strokes at low Reynolds number: a geometric and numerical study of copepod and Purcell swimmers, SIAM Journal on Control and Optimization, 56 (2018), 1794–1822.
  * [5] B. Bonnard, M. Chyba, Singular Trajectories and Their Role in Control Theory, Springer, 2003.
  * [6] B. Bonnard, M. Chyba, J. Rouot, Geometric and Numerical Optimal Control: Application to Swimming at Low Reynolds Number and Magnetic Resonance Imaging, Springer, 2018.
  * [7] B. Bonnard, M. Chyba, J. Rouot, D. Takagi, A numerical approach to the optimal control and efficiency of the copepod swimmer, 2016 IEEE 55th Conference on Decision and Control, (2016), 4196–4201.
  * [8] B. Bonnard, M. Chyba, J. Rouot, D. Takagi, Sub-Riemannian geometry, Hamiltonian dynamics, micro-swimmers, copepod nauplii and copepod robot, Pac. J. Math. Ind., 10 (2018), 1–27.
  * [9] C. J. Bradley, J. R. Strickler, E. J. Buskey, P. H. Lenz, Swimming and escape behavior in two species of calanoid copepods from nauplius to adult, Journal of Plankton Research, 35 (2013), 49–65.
  * [10] E. Bruno, C. M. Andersen Borg, T. Kiørboe, Prey detection and prey capture in copepod nauplii, PLoS ONE, 7 (2012), e47906.
  * [11] D. Devine, Locomotion and rotation with three stiff legs at low Reynolds number, Master’s thesis, University of Hawai’i, Manoa, 2016.
  * [12] R. Dreyfus, J. Baudry, H.A. Stone, Purcell’s “rotator”: Mechanical rotation at low Reynolds number,  Eur. Phys. J. B, 47 (2005), 161–164.
  * [13] W. Gao, et al., Cargo-towing fuel-free magnetic nanoswimmers for targeted drug delivery, Small, 8 (2012), 460–467.
  * [14] W. Hu, G. Lum, M. Mastrangeli, M. Sitti, Small-scale soft-bodied robot with multimodal locomotion, Nature, 554 (2018), 81–85.
  * [15] M.A. Jalali, M. R. Alam, S.Mousavi, Versatile low-Reynolds-number swimmer with three-dimensional maneuverability. Phys. Rev. E , 90 (2014), 053006.
  * [16] P. Lenz, D. Takagi, D. Hartline, Choreographed swimming of copepod nauplii, J. R. Soc. Interface, 12 (2015), 20150776.
  * [17] D. Liberzon, Calculus of Variations and Optimal Control Theory: A Concise Introduction, Princeton University Press, 2012.
  * [18] R. Montgomery, A Tour of Subriemannian Geometries, Their Geodesics and Applications, American Mathematical Society, 2000.
  * [19] K.T.M. Niimoto, K.J. Kuball, L.N. Block, P.H. Lenz, D. Takagi, Rotational maneuvers of copepod nauplii at low Reynolds number, Fluids, 5 (2020).
  * [20] G.-A. Paffenhöfer, J.R. Strickler, K.D. Lewis, S. Richman, Motion behavior of nauplii and early copepodid stages of marine planktonic copepods, Journal of Plankton Research, 18 (1996), 1699–1715.
  * [21] L.S. Pontryagin, V.G. Boltyanskii, R.V. Gamkrelidze, E.F. Mishchenko, The Mathematical Theory of Optimal Processes, 1962.
  * [22] E.M. Purcell, Life at low Reynolds number, American Journal of Physics, 45 (1977), 3–11.
  * [23] M.S. Rizvi, A. Farutin, C. Misbah, Three-bead steering microswimmers, Phys. Rev. E, 97 (2018), 023102.
  * [24] H.E. Robinson, J.R. Strickler, M.J. Henderson, D.K. Hartline, P.H. Lenz, Predation strategies of larval clownfish capturing evasive copepod prey, Mar. Ecol. Prog. Ser., 614 (2019), 125–146.
  * [25] M. Sitti, Miniature soft robots — road to the clinic, Nature Reviews Materials, 3 (2018), 74–75.
  * [26] D. Takagi, Swimming with stiff legs at low Reynolds number, Physical Review, 92 (2015), 023020.
  * [27] Team Commands, Inria Saclay, BOCOP: an open source toolbox for optimal control, http://bocop.org, 2017.
  * [28] J. Turner, The importance of small planktonic copepods and their roles in pelagic marine food webs, Zoological Studies, 43 (2004), 255–266.


# Magnon spin transport through atomic ferrimagnetic domain walls

Zhaozhuo Zeng  Peng Yan<EMAIL_ADDRESS>School of Physics and State Key
Laboratory of Electronic Thin Films and Integrated Devices, University of
Electronic Science and Technology of China, Chengdu 611731, China

###### Abstract

It is a well-established notion that the spin of a magnon should be flipped
when it passes through a $180^{\circ}$ domain wall (DW) in both ferromagnets
and antiferromagnets, while the magnon spin transport through ferrimagnetic DW
is still elusive. In this work, we report that the magnon preserves its spin
after the transmission through an atomically sharp DW in ferrimagnets, due to
the intriguing interband magnon scattering at the domain interface. This
finding may provide significant insight to resolve the puzzling insensitivity
of magnon spin diffusion to the $180^{\circ}$ ferrimagnetic DWs observed by
recent experiments. Our results reveal the unique role of ferrimagnetic DWs in
manipulating the magnon spin and may facilitate the design of novel magnonic
devices based on ferrimagnets.

Magnons (or spin waves) are quasiparticles associated with wave-like
disturbances in ordered magnets, and are able to carry both linear and angular
momentum. The rise of the emerging magnonics is largely due to the low energy
dissipation and long coherence length of magnons serga2010yig ;
chumak2015magnon ; yu2021magnetic ; barman20212021 ; yuan2022quantum . One
fundamental issue is to understand the interaction between magnons and
magnetic textures yan2011all ; wang2012domain ; yu2016magnetic ;
hamalainen2018control ; oh2019bidirectional ; han2019mutual ;
liang2022nonreciprocal ; lan2017antiferromagnetic ; faridi2022atomic ;
iwasaki2014theory ; schutte2014magnon ; kim2019tunable ; wang2021magnonic ;
li2022interaction . Magnetic domain wall (DW), the transition region
separating two magnetic domains, is a prominent example of spin textures that
are promising for information industry schryer1974motion ; wang2018theory . It
has been known that a magnon will transfer its angular momentum when it passes
through the DW, leading to the so-called magnonic spin transfer torque
yan2011all . When the DW is wide, a continuum model is convenient to describe
the magnon transport since it allows analytical solutions. For a narrow DW,
however, the continuum model does not suffice to characterize the spin texture
and an atomistic approach is demanded to interpret, for instance, the atomic
DW pinning novoselov2003subatomic , the strong magnon reflection
yan2012magnonic , and the magnonic Cherenkov emission Yang2019 . Narrow DWs
have been observed in antiferromagnets as well, like CuMnAs
krizek2022atomically  and $\mathrm{FePS_{3}}$ lee2023giant . As their
ferromagnetic counterpart, the spin of a magnon undergoes a sign flipping as
it traverses antiferromagnetic DWs tveten2014antiferromagnetic ;
kim2014propulsion .

Ferrimagnets (FiMs), which exhibit antiferromagnetic coupling but non-zero net
magnetization, can be readily manipulated for ultrafast devices
kim2022ferrimagnetic ; zhang2023ferrimagnets . Due to different angular
momentum densities of two sublattices in FiMs, there are two non-degenerate
circularly polarized states with contrast magnon dispersions even without the
external magnetic field, providing a new degree of spin freedom oh2017coherent
; kim2021current ; kim2020distinct ; okamoto2020flipping . Very recently, a
puzzling insensitivity of magnon spin diffusion in multidomain ferrimagnet has
been reported li2023puzzling . Following the conventional wisdom for both
ferromagnets and antiferromagnets that magnons flip their spin after passing
through the DW, one may expect that the angular momentum carried by magnons in
multi-domain state should quickly decay to zero because of the cancellation
effect in opposite domains. However, non-local measurements show almost the
same signal strength as that in a single domain state li2023puzzling . This
thus bring about a critical question: Do magnon reverse or reserve its spin
after passing through the DW in FiMs?

In this work, we theoretically investigate the angular momentum transport of
magnons through ferrimagnetic DWs in atomic scales. For wide DWs, magnons
indeed switch their spin and meanwhile reserve the linear momentum after
passing through the DW, which echoes the case for both ferromagnet and
antiferromagnet. Interestingly, we find that, for atomically sharp DWs, the
magnon spin remains after the transmission due to the intriguing interband
magnon scattering, which is accompanied by a linear momentum jump. Our results
advance the understanding of the interaction between magnons and magnetic
textures, and may offer stimulating insight to clarify the aforementioned
puzzling insensitivity of magnon spin diffusion to ferrimagnetic DWs.

Figure 1: Spin wave transport through a single domain, wide DW, and sharp DW
in ferrimagnets, respectively. Green circles are the top view of the
magnetization precession. $A$ and $B$ label two sublattices. $\pm\hbar$
represents the magnon spin. $k$ and $k^{\prime}$ indicate the magnon wave
vector.

Let’s consider the following one-dimensional atomistic spin Hamiltonian

$\mathcal{H}=J{\sum\limits_{n}}{\bf{S}}_{n}\cdot{\bf{S}}_{n+1}-K{\sum\limits_{n}}({S}_{n}^{z})^{2},$
(1)

where $J>0$ and $K>0$ represent the antiferromagnetic exchange couple and
easy-axis magnetic anisotropy, respectively. Due to the competition between
these two energies and the boundary condition, there are three kinds of
magnetization profiles, i.e., a single domain, a wide DW, and an atomically
sharp DW, as shown in Fig. 1. It is noted that the width of DW is determined
by the radio $K/J$ (see below).

The equation of motion for each spin with index $n$ can be expressed as
$s_{l}\partial_{t}{\bf{S}}_{n}=-{\bf{S}}_{n}\times{\bf{H}}_{\mathrm{eff},n}$,
where $s_{l}=\mu_{l}/\gamma_{l}$ $(l=\mathrm{A,B})$ is the angular momentum
density, $\mu_{l}$ is the magnetic moment, $\gamma_{l}$ is the gyromagnetic
ratio, and
${\bf{H}}_{\mathrm{eff},n}=-\delta\mathcal{H}_{n}/\delta{\bf{S}}_{n}$ is the
effective field. When the index $n$ is odd (even), $l$ represents A (B)
sublattice. It is worth mentioning that $s_{\mathrm{A}}\neq s_{\mathrm{B}}$ in
general for FiMs, because of which the degeneracy of $\pm\hbar$ magnons is
broken even without the external magnetic field.

We first linearize the equation of motion to obtain the dispersion relation of
magnons. To this end, we express the spin $(S_{n}^{x},S_{n}^{y},S_{n}^{z})$ in
local frames as

$\left(\begin{array}[]{c}S_{n}^{x}\\\ S_{n}^{y}\\\
S_{n}^{z}\end{array}\right)=\left(\begin{array}[]{cccc}1&0&0\\\
0&\cos\theta_{n}&\sin\theta_{n}\\\
0&-\sin\theta_{n}&\cos\theta_{n}\end{array}\right)\left(\begin{array}[]{c}S_{n}^{X}\\\
S_{n}^{Y}\\\ S_{n}^{Z}\end{array}\right),$ (2)

where $S_{n}^{Z}\approx 1,|S_{n}^{X}|,|S_{n}^{Y}|\ll 1$, and $\theta_{n}$ is
the equilibrium deflection angle of the magnetization ${\bf{S}}_{n}$ from the
$z$-axis, satisfying

$J\sin(\theta_{n-1}-\theta_{n})+J\sin(\theta_{n+1}-\theta_{n})+K\sin(2\theta_{n})=0.$
(3)

We then obtain

$\displaystyle s_{l}\partial_{t}S_{n}^{X}=$
$\displaystyle\left[Jc_{n-1}+Jc_{n}-2K(\cos^{2}\theta_{n}-\sin^{2}\theta_{n})\right]S_{n}^{Y}$
(4) $\displaystyle-Jc_{n-1}S_{n-1}^{Y}-Jc_{n}S_{n+1}^{Y},$ $\displaystyle
s_{l}\partial_{t}S_{n}^{Y}=$
$\displaystyle(-Jc_{n-1}-Jc_{n}-2K\cos^{2}\theta_{n})S_{n}^{X}$
$\displaystyle+JS_{n-1}^{X}+JS_{n+1}^{X},$

where $c_{n}=\cos(\theta_{n}-\theta_{n+1})$. To facilitate the analysis, we
define wave functions of right-handed (RH) and left-handed (LH) precession of
local magnetic moments as $\psi_{n}^{\pm}=S_{n}^{X}\pm
iS_{n}^{Y}=\Xi_{s_{l}}^{\pm}\text{exp}[i(\omega t-kna)]$, where
$\Xi_{s_{l}}^{+}$ and $\Xi_{s_{l}}^{-}$ are the amplitudes of RH and LH
precessions, respectively, $k$ is the wave vector, $a$ is the lattice
constant, and $\omega/2\pi$ is the frequency.

For a single domain shown in the top panel of Fig. 1, where
$\theta_{2n-1}(\theta_{2n})=0(\pi)$ with the integer $n$ from $-N$ to $N$, we
derive the dispersion and amplitude ratio for the RH (LH) magnon SM

$\omega_{\pm}=\frac{\mp\zeta+\sqrt{\zeta^{2}+\eta}}{s_{\mathrm{A}}s_{\mathrm{B}}},$
(5)
$\rho^{\mathrm{\pm}}=\frac{\Xi_{s_{\mathrm{A}}}^{\pm}}{\Xi_{s_{\mathrm{B}}}^{\mp}}=-\frac{2J\cos
ka}{2J+2K\mp s_{\mathrm{A}}\omega_{\mathrm{\pm}}},$ (6)

where the RH (LH) magnon takes the negative (positive) sign,
$\zeta=(s_{\mathrm{A}}-s_{\mathrm{B}})(J+K)$ and
$\eta=4s_{\mathrm{A}}s_{\mathrm{B}}[(J+K)^{2}-J^{2}\cos^{2}(ka)]$.

The wide DW.—When $K\ll J$, we have a wide FiM DW. The atomistic model then
can be reduced to a continuum form. We then obtain the equation of motion of
the Néel vector
${\bf{n}}=\left({\bf{S}}_{\mathrm{A}}-{\bf{S}}_{\mathrm{B}}\right)/2$
kim2022ferrimagnetic

$\partial_{t}(\sigma_{s}{\bf{n}}\times\dot{{\bf{n}}}-\delta_{s}{\bf{n}})=\nabla\cdot{\bf{J}}+K_{u}n_{z}{\bf{n}}\times\hat{z},$
(7)

where ${\bf{J}}=A{\bf{n}}\times\nabla{\bf{n}}$ is the magnon spin current
tveten2014antiferromagnetic ; kim2014propulsion , $\sigma_{s}=s^{2}a^{3}/4J$,
$\delta_{s}=s_{\mathrm{A}}-s_{\mathrm{B}}$, $s=s_{\mathrm{A}}+s_{\mathrm{B}}$,
$A=J/a$, and $K_{u}=2K/a^{3}$.

Figure 2:  The distributions of the variation of the $x$\- and $y$-component
of Néel vector $\mathbf{n}$ under $K/J=0.01$ (a) and $K/J=0.3$ (b) for an
incoming RH magnon with $k=0.5$ $\rm{nm^{-1}}$. $T_{\mathrm{M}}$ as a function
of $k$ under different $K/J$ for an incoming RH (c) and LH (d) magnon.

In the local frame, we assume ${\bf{n}}=(\delta n^{X},\delta n^{Y},1)$ in Eq.
(7) with $|\delta n^{X}|,|\delta n^{Y}|\ll 1$, which is justified for small-
angle precession of local magnetic moments. Substituting the Walker ansatz
$\theta(x)=2\arctan[\mathrm{exp}(x/\Delta)]$ into Eq. (7), we obtain

$q^{2}\varphi(\xi)=\Big{(}-\frac{d^{2}}{d\xi^{2}}-2\mathrm{sech}^{2}\xi\Big{)}\varphi(\xi),$
(8)

where $\Delta=\sqrt{A/K_{u}}$ is the DW width,
$q^{2}=(\sigma_{s}\omega^{2}-\delta_{s}\omega)/{K_{u}}-1$, $\varphi=\delta
n^{X}\pm i\delta n^{Y}$, and $\xi=(x-x_{0})/{\Delta}$. The hyperbolic
potential $V=-2\mathrm{sech}^{2}\xi$ yan2011all ; kim2014propulsion  allows a
total transmission of all magnon modes, as shown in Fig. 2(a). We also observe
that $\delta n_{x}$ is ahead (behind) of $\delta n_{y}$ in the left (right)
domain, which clearly indicates an angular momentum switching of magnon, as
schematically plotted in the middle panel of Fig. (1). Figure 2(b) shows that
the amplitude of spin wave decreases when it passes through the domain wall
for $K/J=0.3$, indicating that the practical scalar potential is no longer
reflectionless SM . Meanwhile, when an incident spin wave with the same wave
vector passes through DW of different widths, there are different
transmittances. As the width of DW decreases, the transmittance of the magnon
decreases.

Figure 3:  Spatial profiles of $S^{X}$ with an incoming (a) RH and (b) LH spin
wave at $\omega/2\pi=4.45$ THz with $K/J=1$. Inset: green curves shows the
magnon dispersion SM  where the black (red) dot indicates the spin wave in the
left (right) side of sharp DW. The orange arrow shows the interband scattering
direction. The dependence of $T_{\mathrm{M}}$ on $\omega/2\pi$ with an
incoming (c) RH and (d) LH spin wave. Curves and circles represent the
analytical results (13) and numerical simulations, respectively.

To quantitatively describe the scattering between spin waves and DWs, we
define the transmittance based on $z$-component of magnon spin current
$T_{\mathrm{M}}=J_{z}^{t}/J_{z}$, where $J_{z}^{t}$ is the $z$-component of
transmitted magnon spin current. In atomic lattices, the $z$-component of
magnon spin current $J_{z}=A\left({\bf{n}}\times\nabla{\bf{n}}\right)_{z}$
reads

$J_{z}=\pm\frac{J}{2a^{2}}[(\mathcal{R}_{s_{\mathrm{A}}})^{2}+(\mathcal{R}_{s_{\mathrm{B}}})^{2}-2\mathcal{R}_{s_{\mathrm{A}}}\mathcal{R}_{s_{\mathrm{B}}}\cos(ka)]\sin(2ka),$
(9)

where $\mathcal{R}_{s_{\mathrm{A}}}$ and $\mathcal{R}_{s_{\mathrm{B}}}$ are
the amplitudes of magnon in each sublattice, and the RH (LH) magnon takes on
the negative (positive) sign. Using the micromagnetic simulations, we
determine the transmittance $T_{\mathrm{M}}$ under different $K/J$, as shown
in Figs. 2(c) and 2(d). The negative values of $T_{\mathrm{M}}$ represent
spin-wave chiral flipping after passing through the DW. Meanwhile, for
$K/J\leq 0.5$, $T_{\mathrm{M}}$ is similar for an incoming RH and LH spin
waves, indicating the insensitivity of polarization. However, in the case of
$K/J=0.6$, there exists noticeable difference between $T_{\mathrm{M}}$ of the
RH and LH magnon, which results from the polarization-dependent interband
magnon scattering and preserves the chirality of magnon (see analysis below).

The sharp DW.—It is noted that, when $K/J\geq 2/3$, an abruptly sharp DW
($\theta_{n}=0$ or $\pi$) forms with the spins pointing towards the $\pm z$
axis [see Eq. (15) below], which is plotted in the bottom panel of Fig. (1).
Figures 3(a) and 3(b) show the profiles of $S^{X}$ for an incoming RH and LH
spin wave, respectively, with $\omega/2\pi=4.45$ THz under $K/J=1$. Strikingly
different from the case of wide DWs, we observe that the transmission
coefficient sensitively depends on the polarization of incoming spin waves and
the wave vector $k$ varies when it passes through the DW. As shown in the
inset of Fig. 3(a) and 3(b), the interband magnon scattering induces the
variation of magnon wave vector. Specifically, it is found that the amplitude
of spin wave in sublattice A decreases, while it increases in sublattice B for
an incoming RH spin wave. For the case of an LH spin wave, the amplitudes in
both sublattices decrease after the transmission. To reach a quantitative
understanding, we derive both the magnon reflection and transmission
coefficient below.

First of all, the four spins in the dashed green box of the bottom panel of
Fig. (1) [or the top-right inset in Fig. (3)(b)] are used as the boundary
condition. For an incoming RH spin wave, we obtain

$\begin{array}[]{c}(2K+s_{\mathrm{B}}\omega_{\mathrm{+}})\psi_{0}^{-}+J\psi_{-1}^{+}+J\psi_{1}^{-}=0\\\
(2K+s_{\mathrm{A}}\omega_{\mathrm{+}})\psi_{1}^{-}+J\psi_{0}^{-}+J\psi_{2}^{+}=0\end{array}.$
(10)

From the $\psi_{0}^{-}$ and $\psi_{1}^{-}$ in Eq. (10), it is evident that the
two spins at the domain interface automatically exhibit identical LH
precession xing2022comparison . Hence, the angular momentum of magnon remains
when it passes through the sharp DW. However, due to the different dispersions
for the same polarized spin wave in two sides of DW, there is an interband
magnon scattering, associating with a linear momentum transfer. Meanwhile, it
is worth noting that, because of the huge atomic pinning, the sharp DW cannot
be driven by the linear force yang2019atomic . For $n\leq 0$,
$\psi_{n}^{+}=\Xi_{s_{\mathrm{A}}}^{+}(e^{ikna}+r_{\mathrm{R}}e^{-ikna})$ and
$\psi_{n}^{-}=\Xi_{s_{\mathrm{B}}}^{-}(e^{ikna}+r_{\mathrm{R}}e^{-ikna})$. As
for $n>0$,
$\psi_{n}^{-}=t_{\mathrm{R}}\Xi_{s_{\mathrm{A}}}^{-}e^{ik^{\prime}na}$ and
$\psi_{n}^{+}=t_{\mathrm{R}}\Xi_{s_{\mathrm{B}}}^{+}e^{ik^{\prime}na}$, where
$k$ and $k^{\prime}$ are the wave vector of RH spin wave at the left and right
sides of the sharp DW, respectively. We derive the reflection and transmission
coefficients for RH magnons

$r_{\mathrm{R}}=\frac{\mathcal{A}\mathcal{B}^{\ast}-J^{2}\rho^{\mathrm{-}}}{J^{2}\rho^{\mathrm{-}}-\mathcal{A}\mathcal{B}},$
(11)

and

$t_{\mathrm{R}}=\frac{J^{2}\rho^{\mathrm{+}}(e^{ika}-e^{-ika})}{\beta(J^{2}\rho^{\mathrm{-}}-\mathcal{A}\mathcal{B})e^{ik^{\prime}a}},$
(12)

where $\beta=\Xi_{s_{\mathrm{B}}}^{-}/\Xi_{s_{\mathrm{B}}}^{+}$,
$\mathcal{A}=2K\rho^{\mathrm{-}}+s_{\mathrm{A}}\omega_{\mathrm{+}}\rho^{\mathrm{-}}+Je^{ik^{\prime}a}$
and
$\mathcal{B}=2K+s_{\mathrm{B}}\omega_{\mathrm{+}}+J\rho^{\mathrm{+}}e^{ik^{\prime}a}$.
For a LH spin wave, one can do the following substitution:
$\omega_{\mathrm{+}}\rightarrow-\omega_{\mathrm{-}}$,
$\rho^{\mathrm{+}}\leftrightarrow\rho^{\mathrm{-}}$ in Eq. (12).

In the case of an incident RH spin wave, the transmission of the $z$-component
of magnon spin current $T_{\mathrm{M}}$ can be expressed as

$T_{\mathrm{M}}=\frac{\left|\beta
t_{\mathrm{R}}\right|^{2}[1+(\rho^{-})^{2}-2\rho^{-}\cos(k^{\prime}a)]\sin(2k^{\prime}a)}{[1+(\rho^{+})^{2}-2\rho^{+}\cos(ka)]\sin(2ka)}.$
(13)

Figure 4:  Dependence of $T_{\mathrm{M}}$ (a) and $P_{\mathrm{T}}$ (b) on
$K/J$ for an incident RH magnon with $k=2.0$ $\rm{nm^{-1}}$. Blue and orange
curves represent the spin-wave chiral flipping and reserving process,
respectively. The time evolution of the $x$\- and $y$-components of Néel
vector $\mathbf{n}$ under $K/J=0.1$ (c), $0.487$ (d), and $1.0$ (e) at $x=200$
nm, respectively. The color bar shows the accumulation of the time.

It can be seen that the $T_{\mathrm{M}}$ first increases and then decreases
with the frequency, as illustrated in Figs. 3(c) and 3(d). The two zero points
result from the zero group velocity of the spin wave. Meanwhile, the RH magnon
passes more easily than the LH one, which can be attributed to the distinct
bound states of spin waves associated with the DW faridi2022atomic . In
addition, when the angular momentum densities $s_{\mathrm{A}}$ and
$s_{\mathrm{B}}$ switch, there exist only slight differences of
$T_{\mathrm{M}}$ SM .

To further understand the scattering between the magnon and DW, we fix the
wavelength of the incident RH spin wave and adjust the DW width, as shown in
Fig. 4. It can be seen that the $T_{\mathrm{M}}$ for spin wave chirality
flipping process is negative and saturates to zero as the DW width decreases.
However, for the chirality conserving process, $T_{\mathrm{M}}$ reaches the
maximum at $K/J=2/3$ and then decreases. To interpret it, we consider
fluctuations of a single domain barbara1994magnetization  as

$\mathcal{W}(\theta_{2n}-\pi)=\theta_{2n-1}+\theta_{2n+1},$ (14)

where $\mathcal{W}=2(K/J+1)$. The following solution is found:
$\theta_{2n}=\theta_{\mathrm{s}}e^{-2n\phi}+\pi$ and
$\theta_{2n-1}=\theta_{\mathrm{s}}e^{-(2n-1)\phi}$, where
$\phi=\cosh^{-1}(\mathcal{W}/2)$ and $\theta_{\mathrm{s}}$ is the tilted angle
for the upward spin. Equation (3) in sharp DW then can be recast as

$\frac{4}{3}\Big{[}\frac{K}{J}-1-\frac{(\lambda-1)^{3}}{8}\Big{]}\theta_{\mathrm{s}}^{2}=2\frac{K}{J}-1-\lambda,$
(15)

where $\lambda=e^{-\phi}$. The phase transition emerges from
$\theta_{\mathrm{s}}\neq 0$ to $\theta_{\mathrm{s}}=0$ with $K/J=2/3$. To
describe the ratio of two transmitted spin waves with different chiralities,
we define the proportion
$P_{\mathrm{T}}=|T_{\mathrm{M}}^{s}|/\sum_{s}|T_{\mathrm{M}}^{s}|$, where
$T_{\mathrm{M}}^{s}$ represents the transmittance of the spin-wave chirality
flipping or conserving process, as shown in Fig. 4(b). For $K/J\leq 0.2$, the
chirality flipping process is dominant, which transforms the incoming RH
magnon to a LH one, as depicted in Fig. 4(c). For $K/J\geq 2/3$, the spin-wave
chirality conserving process is dominant, which reserves the chirality of
magnon, as presented in Fig. 4(e). In the region between the wide and sharp
DW, the magnon follows a chiral asymptote. Noteworthily, at $K/J=0.487$, we
observe $P_{\mathrm{T}}=1/2$ that represents a zero magnon spin current, and
corresponds to a linearly polarized magnon state, as shown in Fig. 4(d).
Hence, we can conveniently manipulate the chirality of magnon by tuning the DW
width.

Figure 5:  $T_{\mathrm{M}}$ as a function of the wave vector $k$ under
different net angular momentum densities for an incident RH (a) and LH (b)
magnons with $K/J=1$. The orange dots label the maximum of each curve.
$T_{\mathrm{M}}$ at various net angular momentum densities $\delta_{s}$ for an
incident RH (c) and LH (d) magnons with $k=2.0$ nm-1. Red curves and purple
circles show analytical results (13) and numerical simulations, respectively.

In previous calculations, we have fixed the net spin angular momentum density
$\delta_{s}$ of FiMs. To further modulate the magnon transmittance, we vary
$\delta_{s}$ under a constant $K/J=1$, as shown in Figs. 5(a) and 5(b). As
$\delta_{s}$ increases, the maximum value of transmittance for an incoming RH
and LH magnon first increases and then decreases, as indicated by the orange
dots. Then, we fix the wave vector $k=2.0$ nm-1 of incident spin waves, which
is indicated by the red dotted line in Figs. 5(a) and 5(b). It can be seen
that for an incoming RH magnon the transmittance increases as $\delta_{s}$
increases, see Fig. 5(c). As for the incident LH magnon, the transmittance
$T_{\mathrm{M}}$ first increases and then decreases as $\delta_{s}$ increases,
see Fig. 5(d). Simulation results agree well with our analytical formula.

It is crucial to point out that the fundamental reason for the interband
magnon scattering is merely $\theta_{s}=0$. Consequently, an artificially
sharp DW, stabilized by the RKKY coupling, is sufficient induce a significant
interband magnon scattering SM . The roles of dipolar interaction and
magnetoelastic coupling are interesting issues for future study. Generalizing
the present formalism to higher dimensions is an open question.

In conclusion, we investigated the spin transport of magnons through atomic
ferrimagnetic DWs. In the case of a wide DW, it was found that the magnon
flips its spin, echoing its ferromagnetic and antiferromagnetic counterparts.
However, for sharp DWs, the magnon reserves the spin due to the emerging
interband magnon scattering. In the intermediate region, the transmission
magnons with two spins mix and can even carry zero magnon spin current under
some circumstances. Moreover, we identified that the net angular momentum
density of FiMs can effectively modulate the magnon transmittance. These
findings may greatly advance our understanding of the scattering between
magnons and magnetic textures in FiMs.

###### Acknowledgements.

We thank Z. Jin and Y. Cao for useful discussions. This work was funded by the
National Key R$\&$D Program under Contract No. 2022YFA1402802 and the National
Natural Science Foundation of China (NSFC) (Grants No. 12374103 and No.
12074057).

## References

  * (1) A. Serga, A. Chumak, and B. Hillebrands, YIG magnonics, J. Phys. D 43, 264002 (2010).
  * (2) A. V. Chumak, V. I. Vasyuchka, A. A. Serga, and B. Hillebrands, Magnon spintronics, Nat. Phys. 11, 453 (2015).
  * (3) H. Yu, J. Xiao, and H. Schultheiss, Magnetic texture based magnonics, Phys. Rep. 905, 1 (2021).
  * (4) A. Barman, G. Gubbiotti, S. Ladak, A. O. Adeyeye, M. Krawczyk, J. Gräfe, C. Adelmann, S. Cotofana, A. Naeemi, V. I. Vasyuchka, _et al_., The 2021 magnonics roadmap, J. Phys.: Condens. Matter 33, 413001 (2021).
  * (5) H. Yuan, Y. Cao, A. Kamra, R. A. Duine, and P. Yan, Quantum magnonics: When magnon spintronics meets quantum information science, Phys. Rep. 965, 1 (2022).
  * (6) P. Yan, X. S. Wang, and X. R. Wang, All-magnonic spin-transfer torque and domain wall propagation, Phys. Rev. Lett. 107, 177207 (2011).
  * (7) X. S. Wang, P. Yan, Y. Shen, G. E. Bauer, and X. R. Wang, Domain wall propagation through spin wave emission, Phys. Rev. Lett. 109, 167209 (2012).
  * (8) W. Yu, J. Lan, R. Wu, and J. Xiao, Magnetic Snell’s law and spin-wave fiber with Dzyaloshinskii-Moriya interaction, Phys. Rev. B 94, 140410(R) (2016).
  * (9) S. J. Hämäläinen, M. Madami, H. Qin, G. Gubbiotti, and S. van Dijken, Control of spin-wave transmission by a programmable domain wall, Nat. Commun. 9, 4853 (2018).
  * (10) S.-H. Oh, S. K. Kim, J. Xiao, and K.-J. Lee, Bidirectional spin-wave-driven domain wall motion in ferrimagnets, Phys. Rev. B 100, 174403 (2019).
  * (11) J. Han, P. Zhang, J. T. Hou, S. A. Siddiqui, and L. Liu, Mutual control of coherent spin waves and magnetic domain walls in a magnonic device, Science 366, 1121 (2019).
  * (12) X. Liang, Z. Wang, P. Yan, and Y. Zhou, Nonreciprocal spin waves in ferrimagnetic domain-wall channels, Phys. Rev. B 106, 224413 (2022).
  * (13) J. Lan, W. Yu, and J. Xiao, Antiferromagnetic domain wall as spin wave polarizer and retarder, Nat. Commun. 8, 178 (2017).
  * (14) E. Faridi, S. K. Kim, and G. Vignale, Atomic-scale spin-wave polarizer based on a sharp antiferromagnetic domain wall, Phys. Rev. B 106, 094411 (2022).
  * (15) J. Iwasaki, A. J. Beekman, and N. Nagaosa, Theory of magnon-skyrmion scattering in chiral magnets, Phys. Rev. B 89, 064412 (2014).
  * (16) C. Schütte and M. Garst, Magnon-skyrmion scattering in chiral magnets, Phys. Rev. B 90, 094423 (2014).
  * (17) S. K. Kim, K. Nakata, D. Loss, and Y. Tserkovnyak, Tunable magnonic thermal hall effect in skyrmion crystal phases of ferrimagnets, Phys. Rev. Lett. 122, 057204 (2019).
  * (18) Z. Wang, H. Yuan, Y. Cao, Z.-X. Li, R. A. Duine, and P. Yan, Magnonic frequency comb through nonlinear magnon-skyrmion scattering, Phys. Rev. Lett. 127, 037202 (2021).
  * (19) Z. Li, M. Ma, Z. Chen, K. Xie, and F. Ma, Interaction between magnon and skyrmion: Toward quantum magnonics, J. Appl. Phys. 132, 210702 (2022).
  * (20) N. L. Schryer and L. R. Walker, The motion of 180∘ domain walls in uniform dc magnetic fields, J. Appl. Phys. 45, 5406 (1974).
  * (21) X. S. Wang, H. Y. Yuan, and X. R. Wang, A theory on skyrmion size, Commun. Phys. 1, 31 (2018).
  * (22) K. S. Novoselov, A. K. Geim, S. V. Dubonos, E. W. Hill, and I. V. Grigorieva, Subatomic movements of a domain wall in the Peierls potential, Nature 426, 812 (2003).
  * (23) P. Yan and G. E. Bauer, Magnonic domain wall heat conductance in ferromagnetic wires, Phys. Rev. Lett. 109, 087202 (2012).
  * (24) H. Yang, H. Y. Yuan, M. Yan, H. W. Zhang, and P. Yan, Atomic antiferromagnetic domain wall propagation beyond the relativistic limit, Phys. Rev. B 100, 024407 (2019).
  * (25) F. Krizek, S. Reimers, Z. Kašpar, A. Marmodoro, J. Michalička, O. Man, A. Edström, O. J. Amin, K. W. Edmonds, R. P. Campion, _et al_., Atomically sharp domain walls in an antiferromagnet, Sci. Adv. 8, eabn3535 (2022).
  * (26) Y. Lee, S. Son, C. Kim, S. Kang, J. Shen, M. Kenzelmann, B. Delley, T. Savchenko, S. Parchenko, W. Na, _et al_., Giant magnetic anisotropy in the atomically thin van der Waals antiferromagnet FePS3, Adv. Electron. Mater. 9, 2200650 (2023).
  * (27) E. G. Tveten, A. Qaiumzadeh, and A. Brataas, Antiferromagnetic domain wall motion induced by spin waves, Phys. Rev. Lett. 112, 147204 (2014).
  * (28) S. K. Kim, Y. Tserkovnyak, and O. Tchernyshyov, Propulsion of a domain wall in an antiferromagnet by magnons, Phys. Rev. B 90, 104406 (2014).
  * (29) S. K. Kim, G. S. Beach, K.-J. Lee, T. Ono, T. Rasing, and H. Yang, Ferrimagnetic spintronics, Nat. Mater. 21, 24 (2022).
  * (30) Y. Zhang, X. Feng, Z. Zheng, Z. Zhang, K. Lin, X. Sun, G. Wang, J. Wang, J. Wei, P. Vallobra, _et al_., Ferrimagnets for spintronic devices: From materials to applications, Appl. Phys. Rev. 10, 011301 (2023).
  * (31) S.-H. Oh, S. K. Kim, D.-K. Lee, G. Go, K.-J. Kim, T. Ono, Y. Tserkovnyak, and K.-J. Lee, Coherent terahertz spin-wave emission associated with ferrimagnetic domain wall dynamics, Phys. Rev. B 96, 100407(R) (2017).
  * (32) D.-H. Kim, S.-H. Oh, D.-K. Lee, S. K. Kim, and K.-J. Lee, Current-induced spin-wave Doppler shift and attenuation in compensated ferrimagnets, Phys. Rev. B 103, 014433 (2021).
  * (33) C. Kim, S. Lee, H.-G. Kim, J.-H. Park, K.-W. Moon, J. Y. Park, J. M. Yuk, K.-J. Lee, B.-G. Park, S. K. Kim, _et al_., Distinct handedness of spin wave across the compensation temperatures of ferrimagnets, Nat. Mater. 19, 980 (2020).
  * (34) S. Okamoto, Flipping handedness in ferrimagnets, Nat. Mater. 19, 929 (2020).
  * (35) R. Li, L. J. Riddiford, Y. Chai, M. Dai, H. Zhong, B. Li, P. Li, D. Yi, Y. Zhang, D. A. Broadway, _et al_., A puzzling insensitivity of magnon spin diffusion to the presence of 180-degree domain walls, Nat. Commun. 14, 2393 (2023).
  * (36) See Supplemental Material at http://link.aps.org/supplemental/ for materials parameters adopted in calculations, the derivation of the dispersion of spin wave in ferrimagnets, the magnon transport through sharp FiM domain wall by switching $s_{\mathrm{A}}$ and $s_{\mathrm{B}}$, and the influence of interfacial exchange interaction $J_{12}$.
  * (37) Y. W. Xing, Z. R. Yan, and X. F. Han, Comparison of spin-wave transmission in parallel and antiparallel magnetic configurations, Phys. Rev. B 105, 064427 (2022).
  * (38) H. Yang, H. Yuan, M. Yan, H. Zhang, and P. Yan, Atomic antiferromagnetic domain wall propagation beyond the relativistic limit, Phys. Rev. B 100, 024407 (2019).
  * (39) B. Barbara, Magnetization processes in high anisotropy systems, J. Magn. Magn. Mater. 129, 79 (1994).


# Nonequlibrium Third Law of Thermodynamics

Faezeh Khodabandehlou<EMAIL_ADDRESS>Department of Physics
and Astronomy, KU Leuven, Belgium

###### Abstract

The extended Third Law for nonequilibrium jump processes is studied in [1, 2],
where the nonequilibrium heat capacity and the excess heat with corresponding
quasipotential are introduced. The extended Third Law states that under two
conditions, the nonequlibrium heat capacity vanishes as the temperature
approaches zero. In this paper, we present a concise overview of key papers
addressing the nonequilibrium Nernst postulate (the nonequilibrium Third Law
of thermodaynamics) and nonequilibrium calorimetry. The new result presented
in this work include an interpretation of quasipotential in terms of the mean
time taken by the Markov jump process, starting from a given state, to reach
other states for the first time.

##  1 Introduction

The Third Law of Thermodynamics, initially formulated by Nernst and Planck,
examines fundamental questions about entropy and heat capacity at low
temperatures [3, 4]. Nernst’s version suggests that entropy becomes a constant
at zero temperature, independent of external parameters, and that heat
capacity vanishes as temperature approaches zero [5, 6]. However, the Third
Law or the Nernst postulate is not universally applicable to all substances
and requires nontrivial conditions [7, 8].
Over time, the emphasis has been on developing calorimetry and establishing a
general structure for studying thermodynamic processes in nonequilibrium
systems [9, 10, 11]. The objective is not solely to extend the formal
definition of entropy. Instead, it aims to obtain quantitative measures of
heat changes during quasistatic relaxations in small open systems under new
nonequilibrium conditions. If the change is induced by modifying the
temperature of the environment, the excess heat per change in temperature is
appropriately termed nonequilibrium heat capacity [12, 13, 14]. This need not
be restricted to the context of linear irreversible thermodynamics or close-
to-equilibrium processes. The present paper specifically delves into the low-
temperature behavior of excess heat resulting from changes in system
parameters during such quasistatic relaxations [9, 10, 11, 15, 16].

Our main result indicates that, for a large and physically motivated class of
nonequilibrium systems modeled as Markov jump processes, the excess heat
vanishes at absolute zero. Similar to its equilibrium version, this result
requires two conditions for extending the Third Law: the uniqueness of the
ground state and the sufficient accessibility of states to each other. For the
second condition we will show that in low temperature, it is crucial that the
maximum mean time for the random walker starting from any state to reach other
states remains finite.

##  2 Setup

Nonequilibrium calorimetry, a relatively new subject in nonequilibrium
statistical physics, has been explored in other works, see [14, 9]. However,
in this paper, our focus is on discussing the Nernst Heat Theorem. The
appropriate and reasonable model for describing low-temperature thermal
properties is the Markov jump process. We denote the process at time $t$ as
$X_{t}$.
Consider a finite set of physical states represented by variables like
$x,y,z,\ldots$, each corresponding to energy levels $E(x)$; for instance,
those linked to molecules or lattice gases. In the Markov jump process,
transitions occur from state $x$ to $y$ with a transition rate denoted as
$k(x,y)$. Each state can be assigned a vertex, and an edge exists between $x$
and $y$ only if there is a possible jump between these two states. As the
result, we have a simple finite graph. Let $L$ backward generator of the
Markov jump process, where

$L\,h(x)=\sum_{y}k(x,y)(h(y)-h(x));\quad e^{tL}h\,(x)=\langle
h(X_{t})|X_{0}=x\rangle$ (1)

for any function of $h$. The transition rates contain kinetic and
thermodynamic details, summarizing the weak coupling of the system to a
thermal bath. Typically, these rates may depend on a multi-dimensional vector
in parameter space, see [1]. Nevertheless, in this context, our emphasis is on
the thermal response and excess heat release to the bath when perturbing the
temperature. The rates are taken to satisfy local detailed balance,

$\log\frac{k_{\beta}(x,y)}{k_{\beta}(y,x)}=\beta\,q(x,y)$ (2)

where the antisymmetric $q(x,y)=-q(y,x)$ is the heat dissipated in the thermal
bath at inverse temperature $\beta$; see [17] and references therein.
The rates appear in the master equation, which, in the stationary case, is
expressed as:

$\sum_{y}k_{\beta}(y,x)\rho_{\beta}(y)-k_{\beta}(x,y)\rho_{\beta}(x)=0$ (3)

where $\rho_{\beta}(x)$ is the final solution and representing the unique
stationary distribution of state $x$. The expectation of dissipated power or
heat flux into the thermal bath is

$\mathcal{P}_{\beta}(x)=\sum_{y}k_{\beta}(x,y)\,q(x,y)$ (4)

when in state $x$. The stationary average is denoted by
$\langle\mathcal{P}_{\beta}\rangle_{\beta}=\sum_{x}\mathcal{P}_{\beta}(x)\rho_{\beta}(x)\geq
0$, which we assume to be strictly positive to signify the nonequilibrium
condition.

##  3 Excess heat and heat capacity

Suppose the initial temperature of the thermal bath at time $t=0$ is $\beta$.
If we perturb $\beta$ slightly, $\beta^{\prime}=\beta+\textrm{d}\beta$, the
system undergoes a relaxation from state $x$ to the new steady state,
resulting in the exchange of new excess heat with the thermal bath. We define
the quasi-potential as:

$V_{\beta^{\prime}}(x):=\int_{0}^{\infty}\textrm{d}t\,[\langle\mathcal{P}_{\beta^{\prime}}(X_{t}\,|\,X_{0}=x)\rangle_{\beta^{\prime}}-\langle\mathcal{P}_{\beta^{\prime}}\rangle_{\beta^{\prime}}]$
(5)

where expectations refer to relaxation, resp. stationarity under the
nonequilibrium process for control $\beta^{\prime}$. The time-integral
converges by well-known results on finite-state Markov processes. The excess
heat towards the system is then
$Q^{\text{exc}}_{\beta^{\prime}}(x)=-V_{\beta^{\prime}}(x)$, and we must draw
$x$ from the initial stationary distribution at the original
$\beta=\beta^{\prime}-\textrm{d}\beta$ . Taking that zero-time expectation
over $x$ we find

$\langle Q^{\text{exc}}_{\beta^{\prime}}\rangle_{\beta}:=\langle
V_{\beta}\rangle_{\beta}-\langle V_{\beta^{\prime}}\rangle_{\beta}$

where $\langle V_{\beta}\rangle_{\beta}=0$. Thus, it follows that the change
in excess heat in response to a temperature change corresponds to the same
change in the quasi-potential, or in other words

$\delta Q^{\text{exc}}=-\langle\textrm{d}V\rangle$ (6)

where we left out the subscripts. It is important to see that under detailed
balance, $q_{\beta}(x,y)=E_{\beta}(x)-E_{\beta}(y)$ and $\delta
Q^{\text{eq}}=\textrm{d}\langle
E\rangle^{\text{eq}}-\langle\textrm{d}E\rangle^{\text{eq}}$ as in the First
Law. Since $\langle V\rangle=0$, we have
$-\langle\,\textrm{d}V\,\rangle=\sum_{y}\textrm{d}\rho(y)\,V(y)$ where
$\textrm{d}\rho(y)$ is the change of the (unique) stationary distribution
$\rho$ under $\beta\rightarrow\beta+\textrm{d}\beta$. In other words,

$\delta Q^{\text{exc}}=\sum_{y}\textrm{d}\rho(y)\,V(y)=\langle
V\,;\,\textrm{d}\log\rho\rangle$ (7)

where $\langle f\,;\,g\rangle=\langle fg\rangle-\langle f\rangle\langle
g\rangle$ denotes the stationary covariance, and we used that
$\langle\textrm{d}\log\rho\rangle=0$. Note that in equilibrium, $V=E$ and
$\log\rho=-\beta E-\log Z$ for the partition function $Z$, and the excess heat
(the total heat) gets expressed in terms of the energy-variance, which assures
its positivity. In nonequilibrium regimes, we get two ‘potentials’ (formally,
$V$ and $\log\rho$) and their structure and mutual correlation shape the
excess heat following the response relation (7).

The concept of excess heat naturally leads to the definition of non-
equilibrium heat capacity. The foundational theory and initial examples of
nonequilibrium heat capacities can be found in [12, 13]. Nonequlibrium heat
capacity is

$C(T)=-\Big{\langle}\frac{\textrm{d}}{\textrm{d}T}V\Big{\rangle}_{T}\quad\text{or}\quad
C(\beta)=\beta^{2}\Big{\langle}\frac{\textrm{d}}{\textrm{d}\beta}V\Big{\rangle}_{\beta}$
(8)

which measures the excess heat towards the system per temperature ($=\beta$ in
that case) change. The starting point is simply (6), but now the change in
$\beta$ corresponds to a temperature shift, denoted as $T\rightarrow
T^{\prime}=T+\textrm{d}T$, while keeping other controls constant, such as
energy levels (similar to fixed volume specific heats).

The quasipotential (5) and the heat capacity (8) are calculated in papers [18,
19, 20], where different aspects of non-equilibrium calorimetry are explored.
The quasipotential $V$ can also be expressed as the solution of the Poisson
equation $LV=\mathcal{P}-\langle\mathcal{P}\rangle$. The graphical
representation of $V$ is given in the paper [21], which is utilized for
calculating heat capacity in various physical models, as discussed in [22,
23].

##  4 Nonequlibrium Third Law of thermodynamics

The Third Law asserts that the excess heat, as defined for quasistatic
processes between nonequilibrium conditions, vanishes at absolute zero:

$\delta Q^{\text{exc}}\rightarrow 0$ (9)

As a consequence, the heat capacity needs to vanish at zero temperature as
well, $C(T\downarrow 0)\rightarrow 0$. The right-hand side in (7) depends on
the temperature $T$ via $V$ and via $\rho$. In order to show (9), two things
suffice: 1) $\textrm{d}\rho\rightarrow 0$ for $T\downarrow 0$ and 2) bounded
$V$.
Suppose there exists a unique state $x^{*}$ such that the stationary
distribution $\rho(x)$ converges to $\delta_{x,x^{*}}$ as $T$ approaches zero.
Additionally, assume the smoothness of $\rho$, allowing the interchangeability
of the limit $T$ with the derivative with respect to $\beta$ (or $T$).
Considering both assumption $\textrm{d}\rho$ goes to zero in the limit of
$T\downarrow 0$. If the dominant state is unique, the convergence of
$\textrm{d}\rho$ at low temperatures is exponentially fast. However, if there
are multiple dominant states, $\textrm{d}\rho\rightarrow 0$ may exhibit
polynomial behavior with respect to temperature. The way the stationary
condition $\rho$ influences the excess heat at low temperatures is a
noteworthy example of how heat, even in nonequilibrium statistical mechanics,
reveals insights into low-temperature degeneracy.

Secondly, from (5), it is evident that at low temperatures, the boundedness of
$V$ is linked to relaxation properties and the accessibility of states. To
show the boundedness of $V(x)$ it is enough we show that for every edge
$e=(x,y),$ the diffrences of quasipotentials
$V(e)=V_{\beta^{\prime}}(x)-V_{\beta^{\prime}}(y)$ is bounded [1]. The
subsequent section will demonstrate that the differences in quasipotentials
over an edge are bounded only when the mean time to reach any other states for
the first time from the edge is also bounded. Notably, under global detailed
balanced condition, if in (4) $q(x,y)=E(x)-E(y),$ then $V(x)=E(x)-\langle
E\rangle,$ and the boundedness is immediately guaranteed. In [1], the
boundedness of the quasipotential is investigated through graphical
representation, where the second condition is precisely studied. In the limit
of close to equilibrium the quasi potential is always bounded and the heat
capacity goes to zero as temperature goes to zero, see [24].

##  5 Bound for quasipotential in terms of mean first-passage times

Let $\tau(x,y)$ represent the mean time it takes for a random walker starting
from $x$ to reach $y$ for the first time. In [21], the graphical
representations of the quasipotential and mean first passage time within the
framework of Markov jump processes are provided, where $\tau$ is satisfying
the poisson equation $L\tau=-1$. The graphical representations, in the same
paper, establishe the relationship between quasipotential and mean first
passage, which is

$\displaystyle V(x)-V(y)=\sum_{u}\,\rho(u)\,f(u)[\tau(y,u)-\tau(x,u)]$ (10)

where $f=\mathcal{P}-\langle\mathcal{P}\rangle^{s}$. As a result, the
differences of quasipotential over an edge is bounded if the mean time to
reach any state from the edge, for the first time, is bounded. In simpler
terms, the absence of local traps in the graph ensures this bound. This
interpretation again leads to the second condition of the Third Law, which
pertains to the accessibility of all states in the graph.
In the same paper an upper bound is given to differences of quasi potential,
put $||f||:=\max_{z\in K}|f(z)|$, and the bound

$|V(x)-V(y)|\leq||f||\min\\{\tau(x,y),\tau(y,x)\\}$ (11)

The bound (11) is not always optimal. It still can be improved as is given in
[21].

##  6 Conclusions

Excess heat in nonequilibrium systems, representing absorbed heat during
quasistatic relaxation, defines nonequilibrium heat capacity. Considering
Markov jump processes on a connected graph, we express excess heat in the
quasistatic relaxation to a new nonequilibrium conditionas. The excess heat is
the product of $\textrm{d}\rho$ in stationary probabilities and quasipotential
$V$. To establish our main theorem, the boundedness of the quasipotential in
low temperature is a prerequisite. In the last section, we provide an
expression of quasipotential in terms of mean first passage time, offering an
interpretation of the accessibility of states. Consequently, in the low
temperature, $T\downarrow 0$, the quasipotential remains bounded if the states
are enough accessible together. Our main result asserts that the
nonequilibrium heat capacity tends to zero at low temperatures, contingent on
two conditions: the unique ground state and a bounded quasipotential.

## References

## References

  * [1] Khodabandehlou F, Maes C, Maes I and Netočný K 2023 Annales Henri Poincaré
  * [2] Khodabandehlou F, Maes C and Netočný K 2023 The Journal of Chemical Physics 158 ISSN 1089-7690
  * [3] Nernst W 1926 The New Heat Theorem (Methuen and Company, Ltd.)
  * [4] Planck M 1911 Thermodynamik 3rd ed (De Gruyter)
  * [5] Callen H B 1991 Thermodynamics and an Introduction to Thermostatistics 2nd ed (Wiley)
  * [6] Masanes L and Oppenheim J 2017 Nat. Commun. 8 14538
  * [7] Lieb E 1967 Phys. Rev. 162
  * [8] Aizenman M and Lieb E 1981 J. Stat. Phys. 24 279–297
  * [9] Glansdorff P, Nicolis G and Prigogine I 1974 Proc. Nat. Acad. Sci. USA 71 197–199
  * [10] Oono Y and Paniconi M 1998 Prog. Theor. Phys. Suppl. 130
  * [11] Komatsu T, Nakagawa N, Sasa S and Tasaki H 2008 Phys. Rev. Lett. 100
  * [12] Boksenbojm E, Maes C, Netočný K and Pešek J 2011 Europhys. Lett. 96
  * [13] Pešek J, Boksenbojm E and Netočný K 2012 Cent. Eur. J. Phys. 10 692–701
  * [14] Maes C and Netočný K 2019 J. Stat. Mech. 114004
  * [15] Saito K and Tasaki H 2011 J. Stat. Phys. 145 1275–1290
  * [16] Maes C and Netočný K 2014 J. Stat. Phys. 154 188–203
  * [17] Maes C 2021 SciPost Phys. Lect. Notes 32 17pp
  * [18] Beyen A and Maes I 2023 Phase diagram and specific heat of a nonequilibrium curie-weiss model (Preprint 2307.01795)
  * [19] Dolai P, Maes C and Netočný K 2023 SciPost Physics 14
  * [20] Dolai P and Maes C 2024 Annals of Physics 460 169546 ISSN 0003-4916
  * [21] Khodabandehlou F, Maes C and Netočný K 2023 On the poisson equation for nonreversible markov jump processes (Preprint 2310.19219)
  * [22] Khodabandehlou F and Maes I 2023 Stochastic Processes and their Applications 164 337–356 ISSN 0304-4149
  * [23] Khodabandehlou F, Krekels S and Maes I 2022 Journal of Statistical Mechanics: Theory and Experiment 2022 123208 ISSN 1742-5468
  * [24] Khodabandehlou F and Maes C 2023 Close-to-equilibrium heat capacity (Preprint 2307.11551)