# Odd elastohydrodynamics: non-reciprocal living material in a viscous fluid

Kenta Ishimoto<EMAIL_ADDRESS>Research Institute for
Mathematical Sciences, Kyoto University, Kyoto 606-8502, Japan  Clément Moreau
<EMAIL_ADDRESS>Research Institute for Mathematical Sciences,
Kyoto University, Kyoto 606-8502, Japan  Kento Yasuda
<EMAIL_ADDRESS>Research Institute for Mathematical Sciences,
Kyoto University, Kyoto 606-8502, Japan

(August 28, 2024)

###### Abstract

Motility is a fundamental feature of living matter, encompassing single cells
and collective behavior. Such living systems are characterized by non-
conservativity of energy and a large diversity of spatio-temporal patterns.
Thus, fundamental physical principles to formulate their behavior are not yet
fully understood. This study explores a violation of Newton’s third law in
motile active agents, by considering non-reciprocal mechanical interactions
known as odd elasticity. By extending the description of odd elasticity to a
nonlinear regime, we present a general framework for the swimming dynamics of
active elastic materials in low-Reynolds-number fluids, such as wave-like
patterns observed in eukaryotic cilia and flagella. We investigate the non-
local interactions within a swimmer using generalized material elasticity and
apply these concepts to biological flagellar motion. Through simple solvable
models and the analysis of Chlamydomonas flagella waveforms and experimental
data for human sperm, we demonstrate the wide applicability of a non-local and
non-reciprocal description of internal interactions within living materials in
viscous fluids, offering a unified framework for active and living matter
physics.

##  I Introduction

Motility is one of the main features of living matter, from a single cell to a
swarm of birds or a human crowd [1, 2, 3]. In the last few decades, the
dynamics of motile active agents, both individual and collective behavior,
have been intensively studied, giving rise to a rapidly expanding research
field in physics bridging non-equilibrium statistical physics, biophysics, and
continuum mechanics, now known as active matter and living matter physics. A
crucial feature of these systems is that inner activity units convert energy
into mechanical forces. In turn, Newton’s third law may be violated when we
regard it as an open system, with its mechanical energy being injected from
microscopic active units. Therefore, the mechanical interactions between the
units can be non-reciprocal [4, 5].

The concept of reciprocity is also widely used in continuum mechanics.
Recently, violation of the Maxwell-Betti reciprocity in elasticity has been
discovered in an active system, and termed odd elasticity [6, 7, 8]. The
elastic matrix in the constitutive stress-strain relation is then allowed to
contain non-symmetric components, and it generates a self-sustained
propagating wave. Odd elasticity reflects the non-conservative forces
generated by microscopic active units and provides an effective material
constitutive relation for active and living matter. This formulation was shown
to effectively describe active locomotion as an autonomous system without
controlled, tuned actuation [9].

Motile agents at the cellular scale are usually immersed in viscous fluids and
are self-propelled by their deformation, as seen in swimming microorganisms
[10, 11]. The motility of microswimmers, a term used for active agents in a
low-Reynolds-number fluid, is, however, only possible when their deformations
are non-reciprocal, which is known as the scallop theorem [12, 13, 11].

Recent theoretical studies on the swimming dynamics of odd-elastic materials
[14, 15] revealed the relations between the violation of Maxwell-Betti
reciprocity and the non-reciprocal deformation for microswimming around an
equilibrium configuration, demonstrating that the swimming velocity is
proportional to the magnitude of odd elasticity.

A traveling wave is a typical example of non-reciprocal deformation
ubiquitously observed in biological microswimmers. Indeed, many eukaryotic
cells use a flexible slender appendage, called a flagellum or cilium, for
propulsion by generating a wave. Examples include tail motions of sperm cells
and breaststrokes of Chlamydomonas green algae [16]. This evolutionarily
conserved filament is actuated by inner molecular motors in coordination,
resulting in a periodic traveling wave with a self-organized nature. The
flagellar whip-like motion is therefore regarded as a limit cycle oscillator,
and the generic form of flagellar swimming is provided by Hopf bifurcation
[17]. Recent theoretical and numerical studies using elaborate
elastohydrodynamic models also found the emergence of the various flagellar
waveform patterns via Hopf bifurcation [18, 19, 20, 21]. Moreover, refinements
of videomicroscopy of biological flagella have enabled the detailed analyses
of waveforms, and found that the flagellar shape dynamics are well described
by a noisy limit cycle that reflects internal activity [22, 23, 24, 25, 26,
27].

The self-sustained wave for an odd-elastic material, however, is insufficient
to describe the flagellar waveform, because the odd-elastic waves are
dissipated rather than sustained by the fluid viscosity, similar to the
classical (passive) elastic response in a viscous medium [28]. Hence,
nonlinearity is required for an odd-elastic system to exhibit a stable limit
cycle [15], calling for a more general, nonlinear odd constitutive relation to
deal with biological flagellar swimming. In fact, the importance of nonlinear
odd elasticity has been reported as a topical challenge within the field of
active matter studies.[8]

The primary aim of this study is therefore to extend the odd-elastic
description of microswimmers to a nonlinear regime to deal with stable
periodic deformations, as seen in biological flagellar motion. This theory,
which we call odd elastohydrodynamics, therefore provides a unified framework
for the study of non-local, non-reciprocal interactions of an elastic material
in a viscous fluid.

Using this generic formulation, we can access the interactions inside an
active elastic material, while these are masked by fluid dynamic coupling when
observing flagellar motion under a microscope. To distinguish the non-
reciprocal activity from the passive elastic response, we introduce a new
concept, the odd-elastic modulus, as a spatial Fourier transform in an
extended space. The real and imaginary parts of this complex function possess
proper symmetry and characterize the reciprocal and non-reciprocal
interactions, respectively.

The secondary aim of this study is then to apply our theory to biological
flagellar swimmers. By examining the odd-elastic modulus based on simple
mathematical models and biological experimental data, we show the wide
applicability of a non-local and non-reciprocal description of internal
interactions within living materials.

The contents of this paper are summarized as follows. In Section II, we
provide a setup for the theoretical formulation of odd elastohydrodynamics to
describe an active elastic material in a viscous fluid. We also discuss the
connection between Hopf bifurcation and nonlinear odd elasticity and express
the dynamics of a microswimmer undergoing periodic deformation. In Section
III, we introduce the concept of the odd-elastic modulus.

In Sections IV and V, we apply our theory to understand the inner mechanical
interactions that biological flagellar motion exhibit. To gain physical
intuition regarding non-local, non-reciprocal interactions encoded by
nonlinear odd elasticity, we start with simple and solvable models in Section
IV. We also discuss how the odd-elastic modulus captures the inner
interactions of these example models. In Section V, we numerically investigate
the extended bending modulus in model flagellar waveforms for Chlamydomonas
and sperm cells, together with experimental data. With these, we propose a new
continuum description of living soft matter in a viscous fluid by means of
nonlinear odd elasticity. The discussion and conclusions are provided in
Section VI.

One of the advantages of the odd-elastic description of activity is the
application of the autonomous equations of motion. These allow us to analyze
some general features of microswimming with periodic deformation, including
theoretical formulae for the average swimming velocity. In Appendix A, to
complete our general theory of odd elastohydrodynamics, we further extend our
framework to encompass fluctuations in shape gaits by internal actuation,
following biological observations of a noisy limit cycle in shape space.
Exploiting the autonomous structure of the odd-elastic formulation and the
gauge-field formulation for microswimming, we investigate the effects of
internal active noise on swimming velocity. The role of odd elasticity is
further discussed in terms of non-equilibrium thermodynamics.

##  II Odd elastohydrodynamics of microswimmers

###  II.1 Shape and deformation of a swimmer

To describe the motion of a deforming microswimmer in a fluid, we need to
specify the position and orientation together with the instantaneous shape of
the swimmer. In Fig. 1, we present a schematic of a general elastic
microswimmer. The rigid body motion is defined by the translation and rotation
between the laboratory frame $\\{\bm{e}_{1},\cdots,\bm{e}_{d}\\}$ and swimmer-
fixed frame $\\{\bm{e}^{(s)}_{1},\cdots,\bm{e}^{(s)}_{d}\\}$. We assume that
the swimmer moves in a $d$-dimensional space, where $d=1,2$, or $3$. $d=1$
indicates linear motion, $d=2$ indicates planar motion, and $d=3$ corresponds
to general three-dimensional motion in space. The origin of the swimmer-fixed
frame is set to be the swimmer’s position and is denoted by
$\bm{x}=(x_{1},\cdots,x_{d})^{\textrm{T}}$. The number of angular degrees of
freedom to specify the orientation in $d$-dimensional space is
$d^{\prime}=d(d-1)/2$. We let $n$ be the number of degrees of freedom for
rigid motion, that is, $n=d+d^{\prime}$, and introduce a $n$-dimensional
vector to represent the position and orientation as
$\bm{z}_{0}=(x_{1},\cdots,x_{d},\theta_{1},\cdots,\theta_{d^{\prime}})^{\textrm{T}}\in\mathbb{R}^{n}$.

\begin{overpic}[width=241.84842pt]{fig_model.eps} \end{overpic} FIG. 1:
Schematic of general odd-elastic microswimmer. This example swimmer moves in a
three-dimensional space ($d=3$). The position and rotation of the swimmer are
represented by the relative motions between the laboratory frame
$\\{\bm{e}_{x},\bm{e}_{y},\bm{e}_{z}\\}$ and the swimmer-fixed frame
$\\{\bm{e}^{(s)}_{x},\bm{e}^{(s)}_{y},\bm{e}^{(s)}_{z}\\}$. The shape of the
swimmer is parameterized by $N$ shape variables
$\bm{\sigma}=(\sigma_{1},\cdots,\sigma_{N})$ In this schematic, we use
displacements of material units from equilibrium positions as the shape
variables. In typical linear elastic theory, a recovery force is applied that
is proportional to the displacement using a spring constant, as indicated by
$k$ in this schematic. In the current study, however, we generalize this
elastic force to include internal actuation, which is represented by odd
elasticity.

We assume that the shape of the swimmer is parameterized by $N$ shape
coordinates as
$\bm{\sigma}=(\sigma_{1},\dots,\sigma_{N})^{\textrm{T}}\in\mathbb{R}^{N}$. For
the shape coordinates, we employ, for example, displacements from the
equilibria of the material units or relative angles between neighboring
material units. We will later introduce generalized elastic forces and torques
associated with the shape coordinates (see also Fig. 1).

Let us denote representations in the swimmer-fixed coordinates by superscript
${(s)}$ and introduce the extended coordinates vector
$\bm{z}=(\bm{z}^{(s)}_{0};\bm{\sigma})\in\mathbb{R}^{n+N}$ and its associated
velocity vector,
$\dot{\bm{z}}=(\dot{\bm{z}}^{(s)}_{0};\dot{\bm{\sigma}})\in\mathbb{R}^{n+N}$,
where the semicolon indicates vertical concatenation and the dot symbol
indicates a time derivative. The velocity in the swimmer-fixed coordinates
$\dot{\bm{z}}^{(s)}_{0}$ is a physical quantity that is obtained from the
force and torque balance equations as explained below. The vector
$\bm{z}^{(s)}_{0}$, computed by integrating $\dot{\bm{z}}^{(s)}_{0}$, is
introduced for later use but does not represent a physical position or
orientation when $d=3$ due to the non-commutative nature of the dynamics. We
set the origin of the shape coordinates, $\bm{\sigma}=\bm{0}$, to be the
equilibrium configuration without any internal or external forces.

This description includes, in particular, several minimal mathematical models
of swimmers. For example, Najafi-Golestanian’s three-sphere model [29, 30, 31]
is a swimmer consisting of three spheres connected in a straight line by two
rods and moving in one direction by changing the lengths of the rods; thus,
the degrees of freedom are $(n,N)=(1,2)$. Purcell’s three-link swimmer [12,
32, 33] is another minimal model, which consists of three rods connected by
two hinges to form a snake-like robot, and can swim in a plane by changing the
angles of the hinges. The degrees of freedom are therefore $(n,N)=(3,2)$ for
this model. The shape parameters are the lengths of two arms for the three-
sphere model and the two relative angles for the three-link swimmer.

###  II.2 Odd-elastohydrodynamic equations

The dynamics of a three-dimensional self-deforming elastic object in a viscous
fluid are well represented by the Stokes equation:

$\eta\nabla^{2}\bm{u}=\nabla p,$ (1)

where the velocity field $\bm{u}$ satisfies the incompressibility condition
$\nabla\cdot\bm{u}=0$. Here, $p$ is the pressure field, and the viscosity
$\eta$ is assumed to be constant. Due to the linearity of the Stokes
equations, the hydrodynamic forces and torques conjugate to the extended
coordinates, denoted symbolically by $\bm{f}^{\textrm{hyd}}$, and are
proportional to the time derivative of the extended coordinates. This linear
relation is represented by a positive-definite matrix, called a generalized
grand resistance matrix ${\bf M}$ [34]; hence, $\bm{f}^{\textrm{hyd}}=-{\bf
M}\dot{\bm{z}}$. Due to the negligible inertia, these forces and torques are
balanced by internal or external forces and torques, which we denote by
$\bm{f}$ and introduce below.

We now define an “elasticity” matrix (or equivalently an elastic matrix)
through a general stress-strain constitutive relation as a function of the
shape parameters, ${\bf K}(\bm{\sigma})\in\mathbb{R}^{N\times N}$, to
represent all the internal forces and torques, including the internal activity
force as well as the ordinary passive elastic response. To be more precise,
this generalized elastic matrix is defined by mapping from shape coordinates
to internal forces and torques, given by $\bm{f}=-{\bf
K}(\bm{\sigma})\bm{\sigma}$. This generalized elasticity is reduced to that of
an elastic spring when we take the displacement of the material point for the
shape coordinates and to that of a torque spring (torsion spring) when we
employ the relative angle along a filament as the shape coordinates. At the
equilibrium configuration ($\bm{\sigma}=\bm{0})$, the generalized elastic
force vanishes ($\bm{f}=\bm{0}$.) The non-symmetric part may have non-zero
values, that is, ${\bf K}\neq{\bf K}^{\textrm{T}}$; this corresponds to odd
elasticity and effectively represents the non-conservative, internal activity
of the self-deforming material. If it is linearly odd-elastic, the elastic
matrix ${\bf K}$ is a constant matrix, although it is, in general, determined
by the instantaneous shape of the object. The balance equations for the forces
and torque, $\bm{f}^{\textrm{hyd}}+\bm{f}=\bm{0}$, are therefore summarized in
the following form [15]:

$-{\bf M}(\bm{\sigma})\dot{\bm{z}}={\bf L}(\bm{\sigma})\bm{z}.$ (2)

Note that the matrix ${\bf M}$ only depends on the spontaneous shape of the
swimmer. The right-hand side of Eq. (2) represents a general elastic force,
including both the passive elastic response and internal actuation, and ${\bf
L}$ is a $(n+N)\times(n+N)$ matrix containing the elastic matrix ${\bf K}$ as
$L_{n+\alpha,n+\beta}=K_{\alpha\beta}$ with the other components being zero,
namely, $L_{ij}=L_{i\alpha}=L_{\alpha j}=0$. Throughout this paper, we use
Roman indices such as $i,j=\\{1,\dots,n\\}$ for the translation and
orientation of the object, Greek indices such as
$\alpha,\beta=\\{1,\dots,N\\}$ for the shape coordinates, and the Einstein
summation convention for repeated indices.

By inverting the resistance matrix, we can decompose the shape dynamics from
the rigid body motion in the form

$\dot{\bm{z}_{0}}=-{\bf P}{\bf K}\bm{\sigma}\leavevmode\nobreak\
,\leavevmode\nobreak\ \leavevmode\nobreak\ \dot{\bm{\sigma}}=-{\bf Q}{\bf
K}\bm{\sigma}.$ (3)

The matrices ${\bf P}$ and ${\bf Q}$ are respectively given by
$P_{i\alpha}=N_{i,n+\alpha}$ and $Q_{\alpha\beta}=N_{n+\alpha,n+\beta}$, with
${\bf N}={\bf M}^{-1}$. Note that the second equation (3) provides an
autonomous dynamical system in shape space, and the non-symmetric part of the
elastic matrix plays the role of an internal actuation to drive the
deformation. The first equation determines the translation and rotation of the
swimmer and coincides with the equation of the kinematic swimming problem, in
which the shape gait is a given function.

\begin{overpic}[width=128.0374pt]{fig_hvm_pca1.eps} \put(-4.0,70.0){(a)}
\end{overpic}\begin{overpic}[width=99.58464pt]{fig_hvm_pca2.eps}
\put(-4.0,90.0){(b)} \end{overpic}

\begin{overpic}[width=221.93158pt]{fig_hvm_shape.eps} \put(-9.0,24.0){(c)}
\end{overpic}

FIG. 2: Human sperm swimming as example of microswimmer with noisy limit
cycle. (a) The lowest two PCA modes were obtained from experimental data. The
horizontal axis indicates the normalized arc length along the flagellum from
the head-tail junction. (b) Projections of the shape onto the two-dimensional
PCA shape space. Each circle indicates the shape at a different time. (c)
Superposed snapshots of a swimming human sperm with a time-periodic beat,
obtained by a direct numerical simulation of the Stokes equations. The
waveform was extracted from the experimental data for swimming human sperm as
a limit cycle in the two-dimensional PCA shape space. Figures reproduced from
[26] with permission under the creative commons license
(http://creativecommons.org/licenses/by/4.0/).

###  II.3 Periodic swimming by nonlinear odd elasticity

We now consider a microswimmer undergoing a periodic deformation with a
particular focus on flagellar-like filament dynamics. For eukaryotic flagella,
internal molecular motors synchronously actuate the elastic filament to
generate a periodic waveform. While the emergent waveform is obtained by
elastohydrodynamic mechanical coupling, the onset of wave generation from a
straight equilibrium configuration is well formulated by a Hopf bifurcation
[17, 22, 19, 20].

To illustrate the limit cycle behavior in the shape space, we reproduce in
Fig. 2 the figures on human sperm swimming from Ishimoto et al. [26], in which
principal component analysis (PCA) was performed to reduce the dimensionality
of the flagellar waveform obtained from experimental observations. The authors
found that the flagellar waveforms are well represented by noisy limit cycle
orbits in the two-dimensional shape space spanned by the lowest PCA modes
[Fig. 2(a, b)]. The embedded limit cycle orbit was then extracted and used to
analyze the time-periodic swimming dynamics of human sperm [Fig. 2(c)].

To derive a generic description for such time-periodic swimming, we employ the
normal form of the Hopf bifurcation, which may be written as [35]

$\frac{dz}{dt}=cz+b|z|^{2}z,$ (4)

where $z\in\mathbb{C}$. The parameters $c$ and $b$ are both complex numbers:
$c=\lambda+i\omega\in\mathbb{C}$ and $b=\mu+i\xi\in\mathbb{C}$ with real-
valued parameters $\lambda,\omega,\mu$, and $\xi$. Let us introduce the
apparent shape space in which the shape dynamics are described by the normal
form (4) and denote the shape coordinates in this shape space by $\bm{q}$. The
apparent shape coordinates do not always coincide with the shape coordinates
$\bm{\sigma}$ used in the stress-strain relation. To distinguish $\bm{\sigma}$
from $\bm{q}$, we now refer to $\bm{\sigma}$ as intrinsic shape coordinates.
We then assume the existence of a transformation from the apparent shape
coordinates to the intrinsic shape coordinates given by a full-rank matrix
${\bf W}\in\mathbb{R}^{N\times N}$, that is, $\bm{\sigma}={\bf W}\bm{q}$. The
matrix ${\bf W}$ may be obtained by PCA and we will later examine detailed
construction of the matrix with some examples (Sections IV and V).

From the normal form of Eq. (4), the dynamics in the apparent shape space are
separated into the limit cycle in the $q_{1}-q_{2}$ space and the damping
dynamics in the remaining $(N-2)$ dimensions. Let us introduce the apparent
elastic matrix $\hat{{\bf K}}\in\mathbb{R}^{N\times N}$ to distinguish the
apparent elasticity from the intrinsic elasticity ${\bf K}$ in Eq. (3), and
write the dynamics in the form

$\dot{\bm{q}}=-\hat{\bf K}\bm{q}\leavevmode\nobreak\ \leavevmode\nobreak\
\textrm{with}\leavevmode\nobreak\ \leavevmode\nobreak\
{\bf\hat{K}}(\bm{q})=\begin{pmatrix}{\bf\hat{K}}^{\textrm{LC}}&{\bf O}\\\ {\bf
O}&{\bf\hat{K}}^{\textrm{d}}\end{pmatrix},$ (5)

where the two-dimensional nonlinear elastic matrix
${\bf\hat{K}}^{\textrm{LC}}\in\mathbb{R}^{2\times 2}$ represents the limit
cycle in Eq. (4), and the $(N-2)$-dimensional matrix
${\bf\hat{K}}_{d}\in\mathbb{R}^{(N-2)\times(N-2)}$ in the right-bottom block
expresses the stable modes around the Hopf bifurcation. All the eigenvalues of
${\bf\hat{K}}^{\textrm{d}}$ therefore have non-negative real parts.

After relabeling the parameters $\lambda,\omega,\mu,$ and $\xi$ in Eq. (5) as
$k^{\textrm{e}},k^{\textrm{o}},k^{\textrm{ne}},$ and $k^{\textrm{no}}$ with an
additional minus sign, we may write the components of
${\bf\hat{K}^{\textrm{LC}}}$ as

$\hat{K}^{\textrm{LC}}_{\alpha\beta}=(k^{\textrm{e}}+k^{\textrm{ne}}r^{2})\delta_{\alpha\beta}+(k^{\textrm{o}}+k^{\textrm{no}}r^{2})\epsilon_{\alpha\beta}$
(6)

for $\alpha,\beta\in\\{1,2\\}$. Here, $\delta_{\alpha\beta}$ is the Kronecker
delta, $\epsilon_{\alpha\beta}$ is the two-dimensional Levi-Civita permutation
symbol, and $r=(q_{1}^{2}+q_{2}^{2})^{1/2}$. With these terms, the dynamics in
the apparent shape coordinates are translated from the normal form of the Hopf
bifurcation into dynamics described by odd-elastic interactions. The four
parameters $k^{\textrm{e}},k^{\textrm{o}},k^{\textrm{ne}}$, and
$k^{\textrm{no}}$ are then interpreted as even linear elasticity, odd linear
elasticity, even nonlinear elasticity, and odd nonlinear elasticity,
respectively.

In this Section II.3, we introduced the normal form for the limit cycle in
shape space, which then couples with hydrodynamics to generate the net
displacement, i.e., locomotion or swimming. In this theoretical framework, the
swimming dynamics are fully described by an autonomous system. Hence, by
integrating this system over the cycle of shape deformation, we may obtain a
general formula for the average swimming velocity for a small-amplitude
swimmer [36, 15].

The position and the orientation in $d$ dimensions are represented by an
element of the $d$-dimensional Euclidean group $\textrm{SE}(d)$. We represent
these by $\mathcal{R}\in\mathbb{R}^{n\times n}$ and the time evolution is
provided by its generator, $\mathcal{A}$, via
$\dot{\mathcal{R}}=\mathcal{R}\mathcal{A}$. With the linearity of the Stokes
equation, we may rewrite this generator as
$\mathcal{A}=\mathcal{A}_{\alpha}q_{\alpha}$, and the third-rank tensor
$[\mathcal{A}_{\alpha}]_{ij}=A_{ij\alpha}$ is the connection of the gauge
group $\textrm{SE}(d)$. If the swimmer exhibits a periodic motion with period
$T_{c}$, the displacement and rotation after one beat cycle are obtained by a
loop integral in shape space [36, 37, 38] as

$\mathcal{R}(T)=\mathcal{R}_{0}\bar{\textrm{P}}\exp\left[\int_{0}^{T_{c}}\mathcal{A}(t)\,dt\right]=\mathcal{R}_{0}\overline{\textrm{P}}\exp\left[\oint\mathcal{A}_{\alpha}\,dq_{\alpha}\right],$
(7)

where we write $\mathcal{R}(t=0)=\mathcal{R}_{0}$ and introduce a path-
ordering operator $\overline{\textrm{P}}$. The integral in the last term is
performed over a closed loop in shape space. After expanding for a small
$\bm{q}$ up to its quadratic term, the swimming velocity over one beat cycle
is obtained for a small-amplitude deformation by using a fourth-rank tensor
$\mathcal{F}$, called the curvature of the gauge field, as

$\overline{A_{ij}}=\frac{1}{2}F_{ij\alpha\beta}\overline{q_{\alpha}\dot{q}_{\beta}},$
(8)

where the overline indicates the average over one deformation cycle. Hence,
the average swimming velocity is proportional to the areal velocity enclosed
by the limit cycle in shape space. Using our odd-elastic representation, Eqs.
(5)-(6), the limit cycle exists only when $k^{\textrm{e}}<0$. Then, the
swimming formula can be computed as

$\overline{A_{ij}}=F_{ij12}\left[\frac{k^{\textrm{o}}}{2}\frac{|k^{\textrm{e}}|}{k^{\textrm{ne}}}+\frac{k^{\textrm{no}}}{2}\left(\frac{|k^{\textrm{e}}|}{k^{\textrm{ne}}}\right)^{2}\right].$
(9)

The terms on the right-hand side are proportional to the odd-elastic
coefficients. This equation generalizes the swimming formula for an odd-
elastic swimmer from the linear to the non-linear regime [15].

In addition, as observed in many biological systems, the shape gait, or the
coefficients for the odd elasticity in our framework, temporally fluctuates.
Further, these active fluctuations provides an important link to non-
equilibrium statistical physics. Therefore, internal noise should be taken
into consideration to precisely evaluate locomotion. Thus, to complete our
theory, we thoroughly investigated the impact of noise on the swimming
velocity, extending the swimming formula (9) for noisy limit cycles. The
detailed calculations can be found in Appendix A.

##  III Non-local, non-reciprocal interactions and odd-elastic modulus

In this section, we focus on the interactions between the units of active
material undergoing periodic deformation in a viscous fluid. To characterize
its non-local, non-reciprocal interactions, we introduce the concept of the
odd-elastic modulus.

By changing the variables from $\bm{q}$ to $\bm{\sigma}$, the intrinsic
elastic matrix can be derived from the apparent elastic matrix as

${\bf K}={\bf Q}^{-1}{\bf W}{\bf\hat{K}}{\bf W}^{-1}.$ (10)

As already described, the elastic force or torque including the passive and
active elastic response is symbolically given by

$f_{\alpha}=-K_{\alpha\beta}\sigma_{\beta}.$ (11)

When considering the intrinsic shape coordinates, we usually chose the
displacement from the equilibrium for a material unit and the relative
distance or angle between neighboring material units.

In the example of an odd three-sphere swimmer [14], the intrinsic elastic
matrix is given by

$K_{\alpha\beta}=k^{\textrm{e}}\delta_{\alpha\beta}+k^{\textrm{o}}\epsilon_{\alpha\beta},$
(12)

for $\alpha,\beta\in\\{1,2\\}$. For a Purcell swimmer with odd-elastic hinges
[15], the same intrinsic elastic matrix is considered, where $f_{i}$ indicates
the torque at the $i$-th hinge and $\sigma_{i}$ is the $i$-th relative angle
between neighboring rods.

\begin{overpic}[width=433.62pt]{fig_mat_schmat.eps} \put(-3.0,80.0){{ (a)}}
\put(-3.0,32.0){{ (b)}} \end{overpic} FIG. 3: Schematics of elastic matrix
$K_{\alpha\beta}$ and its continuum representation
$\bar{\kappa}(s,s^{\prime})$ with decomposition into symmetric even elasticity
(reciprocal interaction) and anti-symmetric odd elasticity (non-reciprocal
interaction), (a) ${\bf K}={\bf K}_{\textrm{e}}+{\bf K}_{\textrm{o}}$ and (b)
$\bar{\kappa}=\bar{\kappa}_{\textrm{e}}+\bar{\kappa}_{\textrm{o}}$. The two-
dimensional Fourier transform is associated with a two-dimensional wave vector
$(\nu_{s},\nu_{s^{\prime}})$. To characterize the non-reciprocal nature of the
interactions encoded in the intrinsic elasticity, we consider the Fourier
modes along the diagonal components and in the perpendicular direction,
respectively indicated by $\nu$ and $\hat{\nu}$ in the schematics. The odd-
elastic modulus, Eq. (18), is defined by the Fourier modes in the $\hat{\nu}$
direction.

We further extend this form into a $N\times N$ matrix representation. A
schematic is shown in Fig. 3(a), where the size of the matrix is set to
$N=10$, with the colors indicating the values of the matrix components, which
are chosen arbitrarily for illustration purposes. The off-diagonal components
express non-local interactions. The matrix is then decomposed into its
symmetric and anti-symmetric parts, i.e., ${\bf K}={\bf K}_{\textrm{e}}+{\bf
K}_{\textrm{o}}$. These correspond to the even and odd elastic matrices,
respectively, and the anti-symmetric matrix represents the non-reciprocal
interactions between the units of material.

Let us introduce the Lagrangian coordinate of the material for a point with
shape index $\alpha$, which is denoted by position
$\bm{s}_{\alpha}\in\mathbb{R}^{d}$. For simplicity and later use, here we
focus on a one-dimensional elastic object such as a filament or rod, and
assume its length at rest to be $\ell$. For the Lagrangian coordinates, we
take an equally-spaced arclength along the material at rest and represent it
by $s_{\alpha}\in[0,\ell]$ with its separation being $\Delta\ell=\ell/N$. We
then rewrite Eq. (11), representing the force or torque acting on a material
point with a Lagrange label $s_{\alpha}$, by a non-local interaction
represented by a kernel, $\kappa(s_{\alpha},s_{\beta})$, as

$f(s_{\alpha})=-\sum_{\beta}\kappa(s_{\alpha},s_{\beta})\sigma_{j}(s_{\beta}).$
(13)

For large $N$, it is also useful to consider a continuum representation. By
dividing both terms in (13) by the spatial discretization, we obtain

$\frac{f(s_{\alpha})}{\Delta\ell}=-\sum_{\beta}\frac{\kappa(s_{\alpha},s_{\beta})}{\Delta\ell}\frac{\sigma_{j}(s_{\beta})}{\Delta\ell}\,\Delta\ell.$
(14)

Each term in Eq.(14) represents the density (quantity per unit length):
$\bar{f}=f/\Delta\ell$ and $\bar{\sigma}=\sigma/\Delta\ell$ are the force or
torque per unit length and the displacement per unit length, respectively,
whereas $\bar{\kappa}=\kappa/\Delta\ell$ is interpreted as an extension of the
elastic modulus.

In the continuum representation, as is usually considered, we assume that
these densities are well-behaved functions. Hence, by replacing the summation
with an integral, the continuous form in the large-$N$ limit is obtained as

$\bar{f}(s)=-\int_{0}^{\ell}\bar{\kappa}(s,s^{\prime})\bar{\sigma}(s^{\prime})\,ds^{\prime},$
(15)

where the integral is performed over the material.

A non-energy-conserving active material may possess non-symmetric components
in the elastic matrix ($K_{\alpha\beta}\neq K_{\beta\alpha}$), and this
corresponds to non-reciprocity of non-local elastic interactions:

${\color[rgb]{0,0,0}\bar{\kappa}(s,s^{\prime})\neq\bar{\kappa}(s^{\prime},s)}.$
(16)

As illustrated in Fig. 3(b), we then decompose the kernel function into
reciprocal (even) and non-reciprocal (odd) components as
$\bar{\kappa}(s,s^{\prime})=\bar{\kappa}_{\textrm{e}}(s,s^{\prime})+\bar{\kappa}_{\textrm{o}}(s,s^{\prime})$.
Note that the non-reciprocity expressed by Eq. (16) refers to the breakdown of
the symmetry under the change of two material points,
$(s,s^{\prime})\mapsto(s^{\prime},s)$, and differs from the reciprocity with
respect to the physical coordinates, $(i,j)\mapsto(j,i)$. The macroscopic odd
elasticity often refers to non-reciprocity in the physical coordinates, while
the odd-elastic modulus here refers to microscopic odd-elastic interactions
[6, 8].

To extract the non-reciprocal interactions encoded in the kernel in Eqs. (13)
and (15), we introduce a new physical quantity. Let us first consider the two-
dimensional Fourier transform of the non-local elastic modulus in a space
spanned by $s$ and $s^{\prime}$. The transformation of
${\color[rgb]{0,0,0}\bar{\kappa}}(s,s^{\prime})$ with a wave vector
$(\nu_{s},\nu_{s^{\prime}})$ is then given by

$\tilde{\kappa}{\color[rgb]{0,0,0}(\nu_{s},\nu_{s^{\prime}})}=\int\int{\color[rgb]{0,0,0}\bar{\kappa}}(s,s^{\prime})e^{-i(\nu_{s}s+\nu_{s^{\prime}}s^{\prime})}\,dsds^{\prime},$
(17)

where the integral is performed over the two Lagrangian coordinates, i.e.,
$(s,s^{\prime})\in[0,\ell]\times[0,\ell]$.

To characterize the interactions between the units of material, it is useful
to decompose the wave vector into diagonal and perpendicular components (See
Fig. 3), rather than horizontal and vertical components. The diagonal part of
the elastic matrix, along the $\nu$-direction in Fig. 3, represents spatial
variations of the ordinary elastic response, while the off-diagonal parts
represent non-local interactions. To characterize this non-local behavior, we
consider the wave vector along the $\hat{\nu}$-direction in Fig. 3, which is
perpendicular to the diagonal direction, given by
$\nu_{s}+\nu_{s^{\prime}}=0$.

By plugging this relation into Eq. (17) and introducing
$\hat{\nu}=(\nu_{s}-\nu_{s^{\prime}})/2$, we can obtain a complex function,

$\tilde{\kappa}(\hat{\nu})=\int\int{\color[rgb]{0,0,0}\bar{\kappa}}(s,s^{\prime})e^{-i\hat{\nu}(s-s^{\prime})}\,dsds^{\prime},$
(18)

which we hereafter call the odd-elastic modulus. The physical meaning of this
quantity becomes clearer when the kernel function is decomposed into
reciprocal and non-reciprocal components, which are symmetric and anti-
symmetric, respectively, with respect to the exchange of two material points
$(s,s^{\prime})\mapsto(s^{\prime},s)$ as

$\bar{\kappa}_{\textrm{e}}(s,s^{\prime})=\bar{\kappa}_{\textrm{e}}(s^{\prime},s)\leavevmode\nobreak\
\textrm{and}\leavevmode\nobreak\
\bar{\kappa}_{\textrm{o}}(s,s^{\prime})=-\bar{\kappa}_{\textrm{o}}(s^{\prime},s).$
(19)

Substituting this decomposition into Eq.(18), we have
$\tilde{\kappa}(\hat{\nu})=\tilde{\kappa}_{\textrm{e}}(\hat{\nu})+\tilde{\kappa}_{\textrm{o}}(\hat{\nu})$
with

$\displaystyle\tilde{\kappa}_{\textrm{e}}$ $\displaystyle=$
$\displaystyle\int\int\bar{\kappa}_{\textrm{e}}(s,s^{\prime})e^{-i\hat{\nu}(s-s^{\prime})}\,dsds^{\prime},$
(20) $\displaystyle\tilde{\kappa}_{\textrm{o}}$ $\displaystyle=$
$\displaystyle\int\int\bar{\kappa}_{\textrm{o}}(s,s^{\prime})e^{-i\hat{\nu}(s-s^{\prime})}\,dsds^{\prime}.$
(21)

From the relations in Eq. (19), we find that $\tilde{\kappa}_{\textrm{e}}$ is
a real function, whereas $\tilde{\kappa}_{\textrm{o}}$ is a pure imaginary
function. Hence, the newly introduced odd-elastic modulus allows us to
characterize reciprocal and non-reciprocal interactions using their real and
imaginary components.

The odd-elastic modulus $\tilde{\kappa}$ is equivalent to the Fourier spectrum
if the elastic interactions only depend on $s-s^{\prime}$. By definition, it
is readily found that the real part of the odd-elastic modulus is an even
function of $\hat{\nu}$, and this represents the even elasticity,
characterizing the non-local, reciprocal elastic interactions. The imaginary
part, in contrast, encodes the odd elasticity and non-local, non-reciprocal
elastic interactions, and is an odd function of $\hat{\nu}$.

In later sections, we will examine an active elastic filament, where
$f_{\alpha}$ in Eq. (11) indicates the torque and $\sigma_{\beta}$ encodes the
relative angle as in the case of the Purcell swimmer. Then, the density
$\bar{\kappa}(s,s^{\prime})$ represents an extension of the bending modulus,
while $\bar{\sigma}(s^{\prime})$ corresponds to the local curvature. In this
case, we will call the odd-elastic modulus $\tilde{\kappa}$ an odd-bending
modulus, because it generalizes the linear relation between the torque and
local curvature of the filament.

##  IV Odd-elastic modulus for active filaments

In the previous sections, we examined the properties of a general odd-elastic
material around the Hopf bifurcation. In the following sections, focusing on
one-dimensional active filaments such as the flagella of Chlamydomonas and
sperm cells, we will exploit this framework to analyze the intrinsic elastic
interactions that result in a limit cycle oscillation in the apparent shape
space. We therefore assume that the swimmer shape gait is a known function
obtained, for instance, from experimental observations.

###  IV.1 One-dimensional elastic sphere-spring system

To gain insights into the non-local, non-reciprocal interactions inside an
elastic material, we start with a solvable one-dimensional model. Our first
example is a one-dimensional sphere-spring system, where $N$ spheres of radius
$a$ are connected by elastic springs and form a one-dimensional array as
schematically shown in Fig. 4(a). The system is immersed in a viscous medium
with viscosity $\eta$, and each sphere experiences viscous drag with a drag
coefficient $\gamma=6\pi a\eta$, whereas we neglect the viscous drag on the
elastic springs. We use the Lagrangian label for the sphere, $s_{\alpha}$,
which corresponds to the position of a sphere at rest. We assume a periodic
boundary condition, $s_{N+1}=s_{1}$, and the material points are equally
spaced with discretization, $\Delta\ell=\ell/N$, where $\ell$ is the total
length of the spatial period of the problem.

Let ${\color[rgb]{0,0,0}\zeta}_{\alpha}$ be the displacement of the sphere
labeled by $\alpha$ as the shape variable $\sigma_{\alpha}$. We first consider
a local elastic interaction using elastic springs with a spring constant $k$
between neighboring spheres. The equation of motion is then given by

$m\ddot{{\color[rgb]{0,0,0}\zeta}}_{\alpha}=-\gamma\dot{{\color[rgb]{0,0,0}\zeta}}_{\alpha}+k\left({\color[rgb]{0,0,0}\zeta}_{\alpha+1}-2{\color[rgb]{0,0,0}\zeta}_{\alpha}+{\color[rgb]{0,0,0}\zeta}_{\alpha-1}\right).$
(22)

In the large-$N$ limit, its continuous representation can be obtained, where
we take the $\Delta\ell\rightarrow 0$ limit with the mass density
$\bar{m}=m/\Delta\ell$, drag per unit length $\bar{\gamma}=\gamma/\Delta\ell$
and $\bar{k}=k\Delta\ell$ kept constant. This argument leads to the well-known
continuum equation for the displacement field $\zeta(s,t)$ via

$\bar{m}\frac{\partial^{2}\zeta}{\partial
t^{2}}=-\bar{\gamma}\frac{\partial\zeta}{\partial
t}+\bar{k}\frac{\partial^{2}\zeta}{\partial s^{2}}.$ (23)

What we examine by odd-elastic modulus takes a different approach; we specify
the interactions for a given wave pattern, rather than deriving differential
equations from local microscopic interactions. Hence, instead of Eq. (22), we
consider a general non-local elastic force given by Eq. (11), and our non-
local sphere-spring system as

$m\ddot{{\color[rgb]{0,0,0}\zeta}}_{\alpha}=-\gamma\dot{{\color[rgb]{0,0,0}\zeta}}_{\alpha}{\color[rgb]{0,0,0}-}K_{\alpha\beta}{\color[rgb]{0,0,0}\zeta}_{\beta}.$
(24)

We then solve $K_{\alpha\beta}$ for a given wave pattern
$\zeta_{\alpha}(s_{\alpha},t)$.

\begin{overpic}[width=250.38434pt]{fig_solvable.eps}
\put(0.0,50.0){{\large(a)}} \put(0.0,25.0){{\large(b)}} \end{overpic} FIG. 4:
Schematics of example models for intrinsic elasticity. (a) One-dimensional
sphere-spring system. Spheres of radius $a$ are connected by linear springs.
(b) Elastic filament immersed in a viscous fluid. Assuming a small-amplitude
deformation from a straight line, we parameterize the filament shape by the
height $h$.

We first consider the frictionless case, where $\gamma=0$. Our goal here is to
calculate the kernel function $K_{\alpha\beta}$ that sustains a given
traveling wave. To do so, we consider a traveling wave pattern with wavenumber
$\nu$ and angular frequency $\omega$ as

${\color[rgb]{0,0,0}\zeta}_{\alpha}=A\sin(\nu s_{\alpha}-\omega t)$ (25)

with an arbitrary amplitude $A$. Then, we set an orthogonal basis ${\bf
W}=(\bm{w}^{(1)},\bm{w}^{(2)},\cdots)$, such that $w^{(1)}_{\alpha}=\sin(\nu
s_{\alpha})$ and $w^{(2)}_{\alpha}=\cos(\nu s_{\alpha})$. Through direct
calculations of ${\color[rgb]{0,0,0}\zeta}_{\alpha}$ and
$\ddot{{\color[rgb]{0,0,0}\zeta}}_{\alpha}$ and comparison of both terms in
Eq. (24), we readily obtain

${\bf K}={\color[rgb]{0,0,0}m\omega^{2}}{\bf W}\begin{pmatrix}1&&\\\ &1&\\\
&&{\bf O}\\\ \end{pmatrix}{\bf W}^{\textrm{T}},$ (26)

where the symbol in the bottom right block ${\bf O}$ indicates that non-
designated components are all zero. This leads to a matrix representation of
the non-local elasticity,

$K_{\alpha\beta}={\color[rgb]{0,0,0}m\omega^{2}}\cos[\nu(s_{\alpha}-s_{\beta})].$
(27)

The continuum representation is obtained by taking the large-$N$ limit, with
the mass density $\bar{m}$ (mass per unit length) kept constant, as

${\color[rgb]{0,0,0}\bar{\kappa}}(s,s^{\prime})={\color[rgb]{0,0,0}\bar{m}\omega^{2}}\cos[\nu(s-s^{\prime})].$
(28)

Here, the scaling $\bar{\kappa}=\kappa/\Delta\ell$ is different from that used
for Eq. (23), where $\bar{k}=k\Delta\ell$ is assumed to be constant.

The wavenumber of the propagating wave is selective in the non-local sphere-
spring system, which is qualitatively different from the wave equation (24),
where the interactions are local and waves with an arbitrary frequency may
propagate. These differences become clearer when we solve the non-local
sphere-spring system using the kernel representation of Eq. (27).

To do so, we first consider a time-frequency Fourier transform,
$\zeta_{\alpha}(t)=\int\xi_{\alpha}(p)e^{-ipt}dp$, where $\xi_{\alpha}(p)$ is
a Fourier component with frequency $p$. We then have

$\omega^{2}{\bf W}\begin{pmatrix}1&&\\\ &1&\\\ &&{\bf O}\\\ \end{pmatrix}{\bf
W}^{\textrm{T}}\bm{\xi}(p)=p^{2}\bm{\xi}(p),$ (29)

which forms an eigenvalue problem. This may be exactly solved and the inverse
Fourier transform leads to a general solution, given by

$\displaystyle\zeta_{\alpha}$ $\displaystyle=$ $\displaystyle c_{1}\sin(\nu
s_{\alpha}-\omega t+\varphi_{1})$ (30) $\displaystyle+c_{2}\cos(\nu
s_{\alpha}-\omega
t+\varphi_{2})+\sum_{\alpha=3}^{N}c_{\alpha}w_{\alpha}^{(k)},$

with constants $\varphi_{1}$, $\varphi_{2}$ and $c_{\alpha}$
($\alpha=1,\dots,N)$ determined by initial conditions. As the final term in
Eq. (30) is a time-independent constant, only a sinusoidal wave with
wavenumber $\nu$ can propagate, and thus the wave pattern is robust against
disturbance.

We then consider the same oscillatory behavior for the overdamped limit of the
dynamics in Eq. (24), that is, $m=0$. Due to viscous drag, to sustain the
traveling wave, we need to inject some energy into the system. We therefore
expect that the elastic representation of $K_{\alpha\beta}$ should contain
non-reciprocal, odd components.

Using an analysis similar to that for the frictionless dynamics above, we then
have

${\bf K}={\color[rgb]{0,0,0}\gamma\omega}{\bf W}\begin{pmatrix}&1&\\\ -1&&\\\
&&{\bf O}\\\ \end{pmatrix}{\bf W}^{\textrm{T}},$ (31)

yielding the matrix representation,

$K_{\alpha\beta}={\color[rgb]{0,0,0}\gamma\omega}\sin[\nu(s_{\alpha}-s_{\beta})].$
(32)

The continuum representation is also obtained using similar arguments as

${\color[rgb]{0,0,0}\bar{\kappa}}(s,s^{\prime})={\color[rgb]{0,0,0}\bar{\gamma}\omega}\sin[\nu(s-s^{\prime})],$
(33)

where the drag per unit length $\bar{\gamma}=\gamma/\Delta l$ is assumed to be
constant in the large-$N$ limit.

For overdamped dynamics ($m=0$), the elastic interactions are non-local, as in
the previous case, but no longer reciprocal. The matrix form, Eq. (32), not
only selects a specific wavenumber but also sustains the associated traveling
wave by the non-reciprocal components. Moreover, the elasticity is purely odd
in the sense that the matrix is anti-symmetric, that is, ${\bf K}=-{\bf
K}^{\textrm{T}}$ and
${\color[rgb]{0,0,0}\bar{\kappa}}(s,s^{\prime})=-{\color[rgb]{0,0,0}\bar{\kappa}}(s^{\prime},s)$.

In the underdamped case with non-zero $m$ and $\gamma$, the results are the
sum of the even and odd terms. The necessity for energy injection is
consistent with the existence of the odd elasticity in the interactions, and
the injected energy is dissipated by the viscosity of the medium.

The odd-elastic moduli for the frictionless system and the overdamped system
are simply calculated from the Fourier transform of Eqs. (28) and (33) as

$\tilde{\kappa}(\hat{\nu})=\frac{{\color[rgb]{0,0,0}\bar{m}}\omega^{2}}{2}\left[\delta(\hat{\nu}-\nu)+\delta(\hat{\nu}+\nu)\right]$
(34)

and

$\tilde{\kappa}(\hat{\nu})=\frac{i{\color[rgb]{0,0,0}\bar{\gamma}}\omega}{2}\left[\delta(\hat{\nu}-\nu)-\delta(\hat{\nu}+\nu)\right],$
(35)

respectively, where $\nu$ is the wavenumber for a given wave pattern in Eq.
(25) and $\delta(x)$ is again the Dirac delta function. $\tilde{\kappa}$ is a
real, even function in the frictionless system, while it is a purely
imaginary, odd function in the overdamped system. The wavenumber of the
traveling wave is clearly captured by the singular peaks.

In summary, the sphere-spring example shows that the odd-elastic modulus
captures the internal mechanisms that robustly and selectively sustain a given
wave pattern. Moreover, by considering its real and imaginary parts, we can
distinguish the conservative, reciprocal elastic force from the non-reciprocal
force, the latter of which is associated with the energy input required in a
dissipative environment.

###  IV.2 Small-amplitude elastohydrodynamic filament

As a second example of the use of the intrinsic non-local elastic matrix, we
examine the small-amplitude dynamics of an active filament at low Reynolds
number. This is one of the simplest cases in elastohydrodynamics and has been
studied with regard to various aspects of filament dynamics in viscous fluid
for more than half a century [39, 40, 41, 21]. In particular, the model in
this section is known to be generic near the critical point of the instability
[17].

As shown schematically in Fig. 4(b), we consider an infinitely long filament
and let $h(x,t)$ be the height of the elastic filament from the $x$-axis, with
its projection to the $x$-axis. We impose a periodic boundary condition with a
spatial period of $\ell$ and discretize it by $N$ equally spaced points as in
the previous example. The Lagrange label for the filament is taken by a
projection onto the $x$-axis so that we have $s_{\alpha}=x_{\alpha}$. The
shape variables are therefore the local height,
$\sigma_{\alpha}=h(x_{\alpha})$. Under the assumption of a small-amplitude
oscillation, the continuous limit of the elastohydrodynamics of the filament
follows a partial differential equation given by [17, 42]

$\xi_{\perp}\frac{\partial h}{\partial t}=-\kappa\frac{\partial^{4}h}{\partial
s^{4}}+\frac{\partial f}{\partial s},$ (36)

where $\xi_{\perp}$ is the perpendicular drag coefficient per unit length for
a slender filament, $\kappa$ is the bending modulus, and $f(s,t)$ is the
pairwise force acting on the filament per unit length .

In elastohydrodynamics studies, the driving force $f(s,t)$ is often a given
function; otherwise, finer models for molecular activity are required to
describe the dynamics of $f(s,t)$. In contrast, here, we remove the driving
force term from Eq. (36) and consider instead a non-local, non-reciprocal
bending modulus, via

$\xi_{\perp}\frac{\partial h}{\partial
t}=-\int_{-\infty}^{\infty}{\color[rgb]{0,0,0}\bar{\kappa}}(s,s^{\prime})\frac{\partial^{4}h}{\partial
s^{4}}(s^{\prime})\,ds^{\prime},$ (37)

that will effectively play the role of driving actuation. We then again
examine the kernel function ${\color[rgb]{0,0,0}\bar{\kappa}}(s,s^{\prime})$
that produces a sinusoidal traveling wave $h(s,t)=A\sin(\nu s-\omega t)$. To
sustain this wave, through calculations similar to those given above, we
obtain the non-local elastic kernel as

${\color[rgb]{0,0,0}\bar{\kappa}}(s,s^{\prime})=\frac{\xi_{\perp}\omega}{\nu^{4}}\sin[\nu(s-s^{\prime})].$
(38)

As in the overdamped case of the sphere-spring system, this is non-reciprocal
and purely odd in terms of the exchange of the positions, while no even
components emerge.

Although we use a continuum equation in this example, we can start from its
discrete version, where the relative angles between neighboring angles are
used as shape parameters. This formulation is used in the next section to deal
with more general elastohydrodynamic interactions, such as finite-amplitude
flagellar waveforms, where we perform numerical estimations.

##  V Odd-elastic modulus for biological swimmers

\begin{overpic}[width=241.84842pt]{fig_config.eps} \end{overpic} FIG. 5:
Schematic of coarse-grained representation of elastic filament and typical
dynamics in $q_{1}-q_{2}$ shape space. We represent the flagellum by $N+1$
rods of length $\Delta\ell$, which are connected at each end by elastic
hinges. The shape configuration is then specified by the relative angles
between neighboring rods, denoted by $\sigma_{\alpha}$
($\alpha=1,2,\cdots,N$). The rotation of the filament is described by the
angle between the horizontal axis and the first rod denoted by $\theta$.
[Inset] The shape gait is given by the autonomous system and we show its
typical trajectory in the $q_{1}-q_{2}$ shape space. The flagellar waveform
approaches a periodic pattern described by the stable limit cycle, after
starting from the initial point near the origin.

\begin{overpic}[width=122.34692pt]{fig_sin_shape.eps} \put(-5.0,75.0){(a)}
\end{overpic}\begin{overpic}[width=122.34692pt]{fig_sin1.eps}
\put(-3.0,75.0){(b)} \put(47.0,80.0){$r=0$} \end{overpic}

\begin{overpic}[width=122.34692pt]{fig_sin2.eps} \put(-3.0,75.0){(c)}
\put(47.0,80.0){$r=1$}
\end{overpic}\begin{overpic}[width=122.34692pt]{fig_sin3.eps}
\put(-3.0,75.0){(d)} \put(46.0,80.0){$r=1.5$} \end{overpic}

FIG. 6: Flagellar waveforms and odd-bending modulus along circle with radius
$r=\sqrt{q_{1}^{2}+q_{2}^{2}}$ in two-dimensional shape space for sinusoidal
flagellum with wavenumber $\nu/2\pi=1.5$: (a) superposed waveforms for a
swimming flagellum during one beat cycle with its left end initially located
at $(x,y)=(0,0)$, (b) $r=0$ at the origin, (c) $r=1$ on the limit cycle orbit,
and (d) $r=1.5$ outside the limit cycle. In each panel, the real and imaginary
parts are shown by the blue broken curve and the red dotted curve,
respectively. The stable limit cycle corresponds to a circle of radius $r=1$.

In this section, to further clarify intrinsic non-local interactions in
flagellar swimming, we investigate finite-amplitude biological models of
flagellar waveforms in a representative sperm cell and Chlamydomonas beat, in
addition to experimental data for human sperm flagella provided by Ishimoto et
al. [26].

The elastic properties of a flagellum are often modeled as an Euler-Bernoulli
beam, in which the local elastic moments are proportional to the local
curvature of the filament. To generate the flagellum wave, as seen in the
small-amplitude example, we require internal actuation in the flagellar model.
Our theory of odd elastohydrodynamics developed in Section II, however,
incorporates the internal actuation as non-local generalized elasticity,
together with the passive elastic response. In this section, we numerically
examine the intrinsic elasticity of a swimming filament by analyzing the given
limit cycle dynamics in shape space (see the inset of Fig. 5).

###  V.1 Numerical methods

To implement the flagellar elastohydrodynamics, we model the swimming
flagellum as a slender filament moving in a plane and represent its dynamics
by the linkages model [43], where $N+1$ rods of length $\Delta\ell$ are
connected at each end by elastic hinges to form a single filament. The shape
configuration is then specified by the relative angles between neighboring
rods, denoted by $\sigma_{\alpha}$ ($\alpha=1,2,\cdots,N$; Fig. 5). These
relative angles are discretized representations of the local curvature, and
for a passive elastic filament, a linear elastic torque is applied at each
hinge. Here, however, we generalize this torque to non-local, nonlinear
interactions, analogous to the dynamics of Eq. (5).

The shape gait of the filament is described by a stable limit cycle in the
$q_{1}-q_{2}$ apparent shape space (inset of Fig. 5). To transform the two
representations between $\bm{\sigma}$ and $\bm{q}$, let us introduce
$\bm{w}^{(\alpha)}$ as the lowest $N$ PCA modes obtained either from data
generated through a mathematical model or experimental data. The number of
intrinsic shape coordinates should be smaller than the PCA modes
$N_{\textrm{flag}}$ in the original data. We then expand the intrinsic shape
coordinates in the PCA modes $\bm{w}^{(\alpha)}$ as
$\bm{\sigma}=q_{\alpha}\bm{w}^{(\alpha)}$ or $\bm{\sigma}={\bf W}\bm{q}$ in a
matrix form, where ${\bf
W}=(\bm{w}^{(1)},\bm{w}^{(2)},\cdots,\bm{w}^{(N)})\in\textrm{O}(N)$ form an
orthonormal basis of the shape space. To apply our theory, we consider the
flagellar waveforms described by Eqs. (5)–(6), and assume the form of the
damping modes appearing in the right-bottom block in Eq. (5), given by
${\bf\hat{K}}^{\textrm{d}}=k^{\textrm{d}}{\bf I}_{N-2}$, where
$k^{\textrm{d}}$ is a non-negative constant and ${\bf I}_{N-2}$ is the
$(N-2)$-dimensional identity matrix.

When the flagellar waveform possesses non-zero time-averaged curvature,
$\bm{\sigma}_{0}$, this part cannot be captured by the PCA modes. The odd-
elastic representation of Eq. (3) is therefore extended to the form

$\dot{\bm{\sigma}}=-{\bf Q}{\bf W}{\bf\hat{K}}{\bf
W}^{\textrm{T}}(\bm{\sigma}-\bm{\sigma}_{0}).$ (39)

Also, this change affects the factor $r^{2}$ in Eq. (6) as

$r^{2}=\left(\bm{w}^{(1)}\cdot(\bm{\sigma}-\bm{\sigma}_{0})\right)^{2}+\left(\bm{w}^{(2)}\cdot(\bm{\sigma}-\bm{\sigma}_{0})\right)^{2}.$
(40)

In the following examples, we consider a freely swimming sperm flagellum in
Secs. V.2 and V.4, and a clamped Chlamydomonas flagellum in Sec. V.3. We
neglect the sperm head and Chlamydomonas cell body in these examples in order
to showcase the odd-bending modulus for different wave patterns. The
elastohydrodynamic coupling, represented by the matrix ${\bf Q}$, is
numerically computed by the coarse-grained method based on resistive force
theory [43], for which we used $N=80$ links.

###  V.2 Sinusoidal flagellum model

We start with a representation used as a simple but canonical model of a sperm
flagellum [44], where the local curvature, or relative angle in the
discretized model, at the arc length $s_{\alpha}\in[0,\ell]$ is given by a
sinusoidal function in the form

$\sigma_{\alpha}=C_{1}\sin(\nu s_{\alpha}-\omega t).$ (41)

Here, the constants $C_{1},\nu,$ and $\omega$ are the curvature amplitude,
wavenumber, and beat angular frequency, respectively. This simple sinusoidal
function is not only theoretically useful, but is also representative of many
sperm flagella of marine species [45, 27]. Here, as a reasonable choice to
match the flagellar waveforms, we set $\nu=3\pi$. The corresponding swimming
dynamics are shown in Fig. 6(a).

We then non-dimensionalize the system. We employ the flagellar length for the
unit of the length scale ($\ell=1$) The linear odd elasticity $k^{\textrm{o}}$
is identical to the phase velocity and hence characterizes the timescale of
the beat cycle, and we set $k^{\textrm{o}}=1$. After we fix the length scale
and timescale, the only remaining physical unit is the force scale. We use the
elastic force as the unit for the force scale by setting $k^{\textrm{e}}=1$.
There is, therefore, one dimensionless parameter remaining in the system that
characterizes the ratio between the timescales for the elastic and viscous
responses, and it is known as the sperm number,
$\mathrm{Sp}=\ell(\xi_{\perp}k^{\textrm{o}}/|k^{\textrm{e}}|)^{1/4}$ [43],
once we set $k^{\textrm{ne}}=1$ to fix the value $C_{1}$ and the radius of the
limit cycle. We also find from Eq. (41) that
$k^{\textrm{no}}=k^{\textrm{d}}=0$.

As a typical parameter for swimming sperm flagella, we set $\textrm{Sp}=3$ and
compute the odd-bending modulus obtained from the intrinsic elastic matrix.
Due to the nonlinear nature of the elasticity, the odd-bending modulus
generally depends on the instantaneous configuration. Nonetheless, due to the
rotational symmetry of the dynamics in shape space, the odd-bending modulus is
almost constant along a circle with a constant radius
$r=\sqrt{q_{1}^{2}+q_{2}^{2}}$.

In Fig. 6(b-d), we plot the odd-bending modulus at $r=0$, and its average over
circles of radius $r=1$ and $1.5$, where $r=1$ corresponds to the stable limit
cycle orbit. The odd-bending modulus for the straight configuration ($r=0$)
possesses a negative bending modulus in its reciprocal part around
$\tilde{\nu}/2\pi=\pm 1.5$ [Fig. 6(b)], indicating that the wave pattern
emerges as an instability for the straight rod. The non-reciprocal part also
has a peak around $\tilde{\nu}/2\pi=\pm 1.5$, and this corresponds to the wave
traveling along the flagellum.

On the limit cycle orbit [Fig. 6(c)], remarkably, the real part of the odd-
bending modulus almost vanishes and the non-reciprocal interactions are
dominant, as observed in the overdamped sphere-spring system described in
Section IV. When the orbit moves outside the limit cycle, as in Fig. 6(d), due
to the nonlinearity, the even elasticity is enhanced for the wavenumber
$\tilde{\nu}/2\pi\approx\pm 1.5$, while the non-reciprocal interactions that
cause the flagellar wave to propagate are also strengthened.

###  V.3 Chlamydomonas flagellum model

We now proceed to another type of simple model, reproducing a Chlamydomonas
flagellum, which is one of the most studied flagella beat patterns and
characterized by its asymmetric ciliary beat pattern. According to Geyer et
al. [46], the local curvature of the C. reinhardtii flagellum is well
represented by a simple function,

$\sigma_{\alpha}=C_{0}+C_{1}\sin(\nu s_{\alpha}-\omega t),$ (42)

as in the previous model, but with a non-zero constant $C_{0}$ for the mean
curvature. We set $C_{0}=-1/40$, $\nu=2\pi$, $k^{\textrm{ne}}=25/4$, and the
sperm number $\textrm{Sp}=1$ for a biologically reasonable waveform, as
plotted in Fig. 7(a). Other parameters are the same as in Section V.2. Here,
we clamped the proximal end of the Chlamydomonas flagellum at $(x,y)=(0,0)$
and neglected the cell body and the other flagellum for simplicity. In
computing the elastohydrodynamics, we further imposed the clamped boundary
condition by removing the rows for the rigid body motion in Eq. (2).

\begin{overpic}[width=113.81102pt]{fig_chlamy_shape.eps} \put(-2.0,76.0){(a)}
\end{overpic}\begin{overpic}[width=119.50148pt]{fig_chlamy2.eps}
\put(-4.0,72.0){(b)} \end{overpic}

FIG. 7: Flagellar waveforms and odd-bending modulus along stable limit cycle
for Chlamydomonas model with proximal end clamped at $(x,y)=(0,0)$. (a)
Superposed waveforms during one beat cycle. (b) Real and imaginary parts of
the odd-bending modulus.

The resulting odd-bending modulus averaged over the limit cycle is plotted in
Fig. 7(b). In this asymmetric beat pattern, the even part of the non-local
elastic interactions remains non-zero and takes both positive and negative
signs with peaks at $\hat{\nu}/2\pi\approx 0.5$ and $\hat{\nu}/2\pi\approx
1.5$, respectively. In contrast, the peak for the non-reciprocal interactions
coincides with the wavenumber of the flagellar waveform, indicating that the
odd elasticity drives the wave as non-conservational internal actuation.

###  V.4 Human sperm flagellum model

We now proceed to analyze human sperm data for our active elastic filament.
Here, we approximate the limit cycle orbit as a unit circle in the two-
dimensional flagellar PCA space using Eq. (5) by expanding the shape variable
using the flagellar PCA modes. The parameters used in this section are the
same as those used in Section V.2, except that we use a non-zero value of
$k^{\textrm{d}}=0.1$ to ensure the existence of a stable limit cycle in the
$N$-dimensional shape space. The free swimming behavior is shown in Fig. 8(a)
as superposed snapshots of the flagellum.

\begin{overpic}[width=119.50148pt]{fig_sperm_shape.eps} \put(-2.0,73.0){(a)}
\end{overpic}\begin{overpic}[width=119.50148pt]{fig_spec2.eps}
\put(-4.0,72.0){(b)} \end{overpic}

FIG. 8: Flagellar waveforms and odd-bending modulus along stable limit cycle
in data-driven human sperm model. (a) Superposed waveforms for a simulated
sperm flagellum during one beat cycle with its left end initially located at
$(x,y)=(0,0)$. (b) Real and imaginary parts of the odd-bending modulus.

As in the previous examples, we plot the odd-bending modulus averaged over the
limit cycle in Fig. 8(b). The reciprocal interactions represented as the real
part of the odd-bending modulus can take both positive and negative values
depending on the wavenumber $\hat{\nu}$. The peaks around $\hat{\nu}\approx 0$
indicate the local passive elastic response of the flagellum. The negative
reciprocal elasticity has a peak around a wavenumber $\hat{\nu}/2\pi\approx
2$, where the peaks for the non-reciprocal elastic interactions are also
located. This are similar to the generation of the wave as a Hopf instability,
as shown in Fig. 6(b). These observations, therefore, imply the following
mechanical balance. The internal activity, shown by the peak in the odd
elasticity and negative even elasticity, generates a flagellar wave, which is
relaxed by the passive elasticity characterized by a local even elastic
response.

##  VI Discussion and conclusions

In this study, to formulate the dynamics of a living material in a viscous
fluid, we investigated a general description of swimming under a periodic
limit cycle oscillation by extending the concept of odd elasticity to a
nonlinear regime. By means of a change of basis from intrinsic to apparent
shape coordinates, we reduced the shape dynamics to the normal form of a Hopf
bifurcation, which is in turn mapped to nonlinear odd elasticity. This
formulation, which we refer to as odd elastohydrodynamics, then enables us to
access the internal non-local, non-reciprocal interactions in the intrinsic
shape space. Further, to characterize the internal activity as well as the
passive elastic response, we introduced a new concept, the odd-elastic
modulus, defined by a spatial Fourier transform in an extended space. Of note,
this odd-elastic modulus is distinct from the widely used complex modulus or
dynamic modulus defined in frequency Fourier space [34].

With the help of the autonomous odd-elastic dynamics of an active system, we
were able to examine the general aspects of microswimmer dynamics.
Furthermore, in the Appendices, we  examined in detail the effects of noise
from the swimming gait on the swimming performance by extending the well-known
swimming formula that provides the average swimming velocity to a general
noisy limit cycle. By calculating the probabilistic areal velocity in shape
space, we found that the effect of noise on the odd elasticity is negligibly
small in a small-deformation regime. Further analysis of the noisy limit cycle
in the shape space allowed us to bridge the entropy production and work done
by the odd elasticity in the nonlinear regime, and making it consistent with
the physical interpretation that the internal actuation inside the elastic
material is described by odd elasticity.

Then, we applied our theory to the analysis of the internal interactions of
living organisms, focusing on flagellar swimming. From solvable simple models
to biological flagellar waveforms for Chlamydomonas and sperm cells, we
studied the odd-bending modulus to decipher the non-local, non-reciprocal
inner interactions within the material. In particular, we found that the
swimmers can possess negative reciprocal even elasticity at some spatial
frequencies, indicating mechanical instability by internal actuation. The
imaginary part of the odd-bending modulus is the material non-reciprocal
response and corresponds to the odd elasticity, which represents the speed of
the generated flagellar wave.

For the limit cycle, we found that the even elasticity ceases for some simple
models, suggesting a simple nonlinear description of the material. To
illustrate its usefulness in a biological context, we further analyzed the
intrinsic odd-elastic response by using Chlamydomonas and human sperm flagella
models, deciphering non-local elastic interactions in biological flagella.

It is useful to point out that the current description of active elastic
material includes elastohydrodynamic coupling with the outer viscous
environment. We have in turn expanded the notion of odd elasticity as a
stress-strain linear relation to an effective material constitutive relation
that deals with the activity, elasticity, and fluid dynamics. In this paper,
we have assumed a circular trajectory in the apparent shape space as a limit
cycle. However, by definition, the odd-elastic modulus can be calculated for
any closed loop, indicating potential applicability to a wide range of
biological data.

Not limited to a one-dimensional filament, our methodology is applicable to
higher-dimensional materials such as active elastic membranes and bulk
dynamics. The odd response of materials has also been examined in terms of odd
viscosity and odd viscoelasticity [47, 48, 8], and natural extensions of the
current methods to these odd materials may be expressed as a viscoelastic
force representation, $\bm{f}=-{\bf K}\bm{\sigma}-{\bf J}\bm{\dot{\sigma}}$,
where odd viscosity is encoded by a non-symmetric matrix ${\bf J}$ coupled to
the rate of deformation. Extension to an active elastic matrix in a
viscoelastic medium is also an interesting future direction, where it is
necessary to numerically calculate the hydrodynamic force ${\bf Q}$ from the
viscoelastic fluid equation.

The current methodology is also applicable to wet active matter systems, as
these extended descriptions of active materials could be useful for
simplifying the modeling of elastohydrodynamic interactions between cells.
These modeling methods will therefore contribute to a better understanding of
the underlying principles of collective behavior, in particular when
elastohydrodynamics play an essential role, as reported for sperm population
dynamics [49, 50]. Furthermore, microswimmers change their swimming pattern in
response to the external environment, such as mammalian spermatozoa before and
after capacitation [51, 52], marine sperm cells in chemoattractant gradients
[45, 53], phototactic algae such as Chlamydomonas and Volvox under a light
source[54, 55], and ciliates in response to mechanical stimuli [56, 57, 58].
Our description of odd elasticity will therefore enable unified comparisons
for such diverse waveform morphologies of microswimmers to characterize the
differences among species and individual cells.

## Acknowledgments

K.I. acknowledges the Japan Society for the Promotion of Science (JSPS)
KAKENHI for Transformative Research Areas A (Grant No. 21H05309), and the
Japan Science and Technology Agency (JST), FOREST (Grant No. JPMJFR212N). C.M.
is a JSPS International Research Fellow (PE22023) and acknowledges funding
support by JSPS (Grant No. 22KF0197). K.Y. acknowledges support by a JSPS
Grant-in-Aid for JSPS Fellows (Grant No. 22KJ1640). K.I., C.M., and K.Y. were
supported in part by the Research Institute for Mathematical Sciences, an
International Joint Usage/Research Center located at Kyoto University.

##  Appendix A Swimming with noisy limit cycle

In this Appendix, to complete the general theory of odd elastohydrodynamics in
the presence of noise caused by internal actuation, we first extend the
swimming formula for the average swimming velocity to a temporally fluctuating
swimming gait, following biological observations of a noisy limit cycle in
shape space. Using the gauge-field formulation for microswimming, we
investigate the effects of internal active noise on swimming velocity. The
role of odd elasticity is then further discussed in terms of non-equilibrium
thermodynamics.

###  A.1 Swimming with probability current

We now consider the swimming formula in Eq. (8) in a statistical sense [15].
With a bracket indicating an ensemble average, the average swimming formula
becomes

$\langle\mathcal{A}\rangle=\frac{1}{2}F_{ij\alpha\beta}\langle
q_{\alpha}\dot{q}_{\beta}\rangle=\textrm{Tr}(\mathcal{F}{\bf J}),$ (43)

where the trace is taken over the shape components. The anti-symmetric matrix
${\bf J}$ is the probabilistic areal velocity matrix, given by

$J_{\alpha\beta}=\bigg{\langle}\oint q_{\alpha}\,dq_{\beta}\bigg{\rangle}.$
(44)

In the two-dimensional shape space, this statistical swimming formula (43)
simply reads as

$\langle A_{ij}\rangle=2F_{ij12}J$ (45)

if we write $J_{\alpha\beta}=J\epsilon_{\alpha\beta}$. The swimming speed in
the form of a gauge potential is proportional to the probabilistic areal
velocity $J$.

To examine the effects of a noisy shape gait on the swimming velocity, we
therefore need to evaluate the value of ${\bf J}$ by introducing stochastic
dynamics in shape space.

We consider an $N$-dimensional autonomous system in the apparent shape space
(5) with Gaussian white noise, given by stochastic differential equations
(SDE) in the sense of Stratonovich in the form

$\frac{dq_{\alpha}}{dt}=f_{\alpha}(\bm{q})+G_{\alpha\beta}(\bm{q})\zeta_{\beta}(t),$
(46)

where the noise has a zero-mean normal Gaussian form, that is,
$\langle\zeta_{i}(t)\rangle=0$ and
$\langle\zeta_{\alpha}(t)\zeta_{\beta}(0)\rangle=\delta_{\alpha\beta}\delta(t)$
with $\delta(t)$ denoting the Dirac delta function. The function $f_{\alpha}$
corresponds to the generalized force and torque in the apparent shape
coordinates and is provided by Eq. (5) as
$f_{\alpha}=-\hat{K}_{\alpha\beta}q_{\beta}$. The diffusion tensor is
introduced as ${\bf D}=(1/2){\bf G}{\bf G}^{\textrm{T}}$ and is symmetric and
positive-definite by definition. The corresponding Fokker-Planck equation for
the probability distribution function $P(\bm{q},t)$ and the probability
current $\bm{j}(\bm{q},t)$ is

$\frac{\partial P}{\partial t}+\nabla\cdot\bm{j}=0,$ (47)

which is explicitly given in the sense of Stratonovich by

$\frac{\partial P}{\partial t}=-\frac{\partial}{\partial
q_{\alpha}}\left[f_{\alpha}P-\frac{1}{2}\left(G_{\alpha\gamma}\frac{\partial}{\partial
q_{\beta}}\left(G_{\beta\gamma}P\right)\right)\right],$ (48)

where $\bm{j}$ is provided by the terms in the bracket on the right-hand side
of Eq. (48).

For brevity, we here only consider the dynamics in the two-dimensional
$q_{1}-q_{2}$ shape space on which the stable limit cycle is located. This
simplification is equivalent to the assumption of ${\bf K}^{\textrm{d}}={\bf
0}$ in Eq. (5). When we relax this assumption to include $N-2$ stable modes,
the swimming formula in Eq. (45) is calculated as the sum of contributions
from other dimensions as discussed in detail in Appendix C of Ishimoto et al.
[15]. We may then write the noise tensor ${\bf G}$ as

${\bf
G}=g_{r}(r)\sqrt{2D_{r}}\,\bm{e}_{r}\otimes\bm{e}_{r}+rg_{\theta}(r)\sqrt{2D_{\theta}}\,\bm{e}_{\theta}\otimes\bm{e}_{\theta},$
(49)

with the unit bases denoted by $\bm{e}_{r}$ and $\bm{e}_{\theta}$ for the
radial and angle directions, respectively. The dynamics of Eq. (46),
represented in Cartesian apparent shape coordinates, is then reduced to a set
of equations in polar coordinates, in the sense of Stratonovich, as

$\frac{dr}{dt}=f_{r}(r)+g_{r}(r)\sqrt{2D_{r}}\zeta_{r},\leavevmode\nobreak\
\frac{d\theta}{dt}=f_{\theta}(r)+g_{\theta}(r)\sqrt{2D_{\theta}}\,\zeta_{\theta}.$
(50)

Here, the suffixes $r$ and $\theta$ indicate the radial and angle coordinates,
respectively, and the zero-mean noise satisfies the relation
$\langle\zeta_{a}(t)\zeta_{b}(0)\rangle=\delta_{ab}\delta(t)$ for
$a,b\in\\{r,\theta\\}$. We then rewrite the probability current $\bm{j}$ in
polar coordinates as

$\displaystyle\bm{j}$ $\displaystyle=$
$\displaystyle\left[\left(f_{r}(r)+D_{r}g_{r}(r)g_{r}^{\prime}(r)\right)P-D_{r}\frac{\partial}{\partial
r}\left([g_{r}(r)]^{2}P\right)\right]\bm{e}_{r}$ (51)
$\displaystyle+\left[rf_{\theta}(r)P-D_{\theta}\frac{\partial}{\partial\theta}\left(rg_{\theta}(r)P\right)\right]\bm{e}_{\theta},$

yielding the Fokker-Planck equation (48) in the form

$\displaystyle\frac{\partial P}{\partial t}$ $\displaystyle=$
$\displaystyle-\frac{\partial}{\partial
r}\left[\left(f_{r}+D_{r}g_{r}g_{r}^{\prime}\right)P-D_{r}\frac{\partial}{\partial
r}\left(g_{r}^{2}P\right)\right]$ (52) $\displaystyle-f_{\theta}\frac{\partial
P}{\partial\theta}+D_{\theta}g_{\theta}\frac{\partial^{2}P}{\partial\theta^{2}},$

where the prime symbol denotes the derivative with respect to $r$. From the
rotational symmetry of the system (50), the steady distribution
$P_{\textrm{st}}$ is independent of the angle coordinate. By using the
boundary condition for the steady distribution satisfying a zero distribution
at infinity, we obtain the following relation for the steady distribution,
after once integrating over $r$:

$\left(f_{r}(r)+D_{r}g_{r}(r)g_{r}^{\prime}(r)\right)P_{\textrm{st}}-D_{r}\frac{d}{dr}\left(\\{g_{r}(r)\\}^{2}P_{\textrm{st}}\right)=0.$
(53)

The steady distribution is then formally solved as

$P_{\textrm{st}}(r)=\frac{N_{\textrm{st}}}{g_{r}(r)}\exp\left[\int^{r}\frac{f_{r}(x)}{D_{r}\\{g_{r}(x)\\}^{2}}dx\right],$
(54)

where $N_{\textrm{st}}$ is the normalization factor. Note that the steady
distribution is known to be unique [59] and here it is  independent of both
the angle coordinates and odd elasticity because of the form
$f_{r}=-k^{\textrm{e}}r-k^{\textrm{ne}}r^{3}$. For a steady distribution, the
probability current $\bm{j}$ (51) possesses only angular components due to the
condition (53), leading to the form

$\bm{j}_{\textrm{st}}=rf_{\theta}(r)P_{\textrm{st}}(r)\bm{e}_{\theta},$ (55)

with $f_{\theta}(r)=k^{\textrm{o}}+k^{\textrm{no}}r^{2}$ in our dynamics. The
non-vanishing probability current characterizes the violation of the detailed
balance or the non-reciprocity of the non-equilibrium system, and has been
studied from several perspectives, such as irreversible circulation [60], curl
flux [61, 62], and non-reciprocity [63, 64].

To calculate the value of $J$, we use the internal noise model proposed by Ma
et al. [22] based on experimental observation of bull sperm. Their study found
that the fluctuating sperm flagellar waveforms were well described by the SDE
model with multiplicative white noise given by

$g_{r}=r\leavevmode\nobreak\ \leavevmode\nobreak\
\textrm{and}\leavevmode\nobreak\ \leavevmode\nobreak\ g_{\theta}=1,$ (56)

where $D_{r}$ and $D_{\theta}$ correspond to the diffusion constant in the
amplitude and angle coordinates, respectively. This noise term in Eq. (56) is
obtained from additive noise in the (linear) even and odd elastic constants
$k^{\textrm{e}}$ and $k^{\textrm{ne}}$ by the transformations
$k^{\textrm{e}}\mapsto k^{\textrm{o}}+\sqrt{2D_{r}}\,\zeta_{r}(t)$ and
$k^{\textrm{o}}\mapsto k^{\textrm{e}}+\sqrt{2D_{\theta}}\,\zeta_{\theta}(t)$
in the deterministic dynamics of Eq. (5) [65]. Because $k^{\textrm{e}}$ and
$k^{\textrm{o}}$ represent the amplitude and phase speed of the limit cycle,
respectively, the noise strengths $D_{r}$ and $D_{\theta}$ therefore
correspond to the amplitude and angle diffusion.

\begin{overpic}[width=136.5733pt]{fig_dist1.eps} \put(5.0,68.0){(a)}
\end{overpic}\begin{overpic}[width=136.5733pt]{fig_dist2.eps}
\put(5.0,68.0){(b)} \end{overpic}

FIG. 9: Steady probability distribution and probabilistic current for system
with noisy limit cycle with biologically relevant multiplicative noise.
Parameters are set as
$k^{\textrm{e}}=-1,k^{\textrm{ne}}=1,k^{\textrm{o}}=1,k^{\textrm{no}}=1$, and
$D_{r}=0.3$. (a) Steady distribution function with a volcano shape. The
distribution function has a degenerate peak and reaches zero at the origin and
infinity. The function is normalized so that the integrated probability
becomes unity. (b) Distribution function projected onto the $q_{1}-q_{2}$
plane and superposed with probabilistic current vectors. The circular current
has a maximum strength around the peak of the distribution.

The detailed form of the steady solution depends on the sign of
$k^{\textrm{e}}$ [65]. Indeed, when $k^{\textrm{e}}\geq 0$, the origin is
stable and the possible steady solution is simply the Dirac delta distribution
at the origin,

$P_{\textrm{st}}(r)=\frac{1}{\pi}\delta(r).$ (57)

In the case with $k^{\textrm{e}}<0$, however, the origin is unstable and a
steady limit cycle emerges, leading to a steady distribution in the form

$P_{\textrm{st}}=N_{\textrm{st}}r^{\frac{|k^{\textrm{e}}|}{D_{r}}-1}\exp\left[-\frac{k^{\textrm{ne}}}{2D_{r}}r^{2}\right]$
(58)

with a normalization prefactor $N_{\textrm{st}}$. This probability function
vanishes at the origin when the noise is sufficiently small,
$D_{r}<|k^{\textrm{e}}|$, whereas the distribution becomes singular at the
origin once the magnitude of the noise is increased to
$D_{r}>|k^{\textrm{e}}|$. From the steady distribution with
$|k^{\textrm{e}}|>D_{r}$, as plotted in Fig. 9(a), it can be seen that for the
volcano-like function, the maximum values are located on a circle of radius
$r^{\ast}=\sqrt{(|k^{\textrm{e}}|-D_{r})/k^{\textrm{ne}}}$, which
monotonically decreases as the diffusion $D_{r}$ increases.

Note that the distribution in Eq. (57) is always possible because the origin
is a stationary solution of the stochastic system with our choice of noise
(56). The resulting steady distribution is therefore obtained by the summation
of (57) and (58), with their weights depending on the initial distribution.

Let us assume $k^{\textrm{e}}<0$ and that the initial distribution does not
contain the $r=0$ state. Then, the steady distribution can be further
expressed with an explicit form of the normalization factor $N_{\textrm{st}}$,
given by

$N_{\textrm{st}}=\frac{\sqrt{2}}{\pi}\left(\sqrt{\frac{k^{\textrm{ne}}}{2D_{r}}}\right)^{\frac{|k^{\textrm{e}}|}{D_{r}}-1}\bigg{/}\leavevmode\nobreak\
\Gamma\left(\frac{|k^{\textrm{e}}|}{2D_{r}}\right),$ (59)

where $\Gamma(x)$ in the denominator indicates the gamma function. Then, the
steady-state probability current $\bm{j}_{\textrm{st}}$ in Eq. 55 possesses a
non-zero value. The plots in Fig. 9(b) show the rotational probability current
superposed on the two-dimensional projection of the probability distribution.
The size of $\bm{j}$ has a maximum value on a circle, which roughly overlaps
with the ridge of the distribution function. Also, we can calculate $J$ as

$\displaystyle J$ $\displaystyle=$ $\displaystyle\bigg{\langle}\oint
q_{1}\,dq_{2}\bigg{\rangle}=\frac{k^{\textrm{o}}}{2}\langle
r^{2}\rangle+\frac{k^{\textrm{no}}}{2}\langle r^{4}\rangle$ (60)
$\displaystyle=$
$\displaystyle\frac{k^{\textrm{o}}}{2}\frac{|k^{\textrm{e}}|}{k^{\textrm{ne}}}+\frac{k^{\textrm{no}}}{2}\frac{|k^{\textrm{e}}|}{k^{\textrm{ne}}}\left(\frac{|k^{\textrm{e}}|}{k^{\textrm{ne}}}+\frac{2D_{r}}{k^{\textrm{ne}}}\right).$

The first term is proportional to the linear odd elasticity $k^{\textrm{o}}$
and is equivalent to the deterministic case. The noise effects only appear in
the nonlinear odd-elastic term and increase (decrease) by amplifying the noise
magnitude for $k^{\textrm{no}}>0$ ($k^{\textrm{no}}<0$), although the
nonlinear effects are sub-dominant compared with the linear odd-elastic term.
Note that this dependence is a characteristic feature of the noise form in Eq.
(56), and the stochastic model with simple additive noise provides
qualitatively similar results for a small noise case but different results in
general (see Appendix B for details).

The non-zero probability current due to the noisy limit cycle is generated
only by the odd parts of the elasticity matrix, and from the swimming formula
(45), the average swimming velocity is also proportional to the size of the
odd elasticity. When $k^{\textrm{e}}>0$, however, the steady probability
distribution shrinks to the origin and the steady-state probability current
vanishes, yielding an average swimming velocity of zero. Moreover, the noise
effect does not appear in the dominant linear odd-elastic term, indicating
that the fluctuating shape gait does not impact swimming under the noisy limit
cycle. However, the effects of noise depend on the stochastic model, and some
qualitative differences are discussed in Appendix B.

In summary, the autonomous odd-elastic dynamics of an active system and the
gauge-field formulation of microswimming enable us to examine the effects of
noise due to the swimming gait on the swimming performance. Notably, we have
found that the effect of noise on the odd elasticity is negligibly small in a
small-deformation regime.

###  A.2 Non-reciprocity, irreversibility, and entropy production

To conclude the general theory of odd elastohydrodynamics, we now examine the
irreversible stochastic dynamics driven by the odd elasticity from the point
of view of thermodynamics. The non-zero probability current is due to the
violation of the detailed balance, and this can be characterized by a non-zero
(positive) entropy production rate in the non-equilibrium statistical physics
[66]. The averaged entropy production rate, $\dot{e}_{p}$, is defined by using
non-equilibrium entropy
$S=\int_{-\infty}^{\infty}\int_{-\infty}^{\infty}P_{\textrm{st}}\log
P_{\textrm{st}}\,dq_{1}\,dq_{2}$, the average rate of heat $\dot{Q}$, and the
system temperature $T$ as

$\dot{e}_{p}\mathrel{\mathop{\mathchar 58\relax}}=\dot{S}-\frac{\dot{Q}}{T},$
(61)

which is, in general, non-negative by the thermodynamic law. For the Langevin
system given by Eq. (46), the entropy production rate is provided by [67]

$\dot{e}_{p}=\langle v_{\alpha}D^{-1}_{\alpha\beta}v_{\beta}\rangle,$ (62)

where $\bm{v}$ is the probabilistic velocity for the steady state defined as
$\bm{v}=\bm{j}_{\textrm{st}}/P_{\textrm{st}}$; hence, in our problem, simply
$\bm{v}=rf_{\theta}(r)\bm{e}_{\theta}$ from Eq. (55). By directly calculating
the diffusion tensor via Eq.(49), we obtain

$\dot{e}_{p}=\frac{1}{D_{\theta}}\langle f^{2}_{\theta}\rangle,$ (63)

which can be calculated, with
$f_{\theta}(r)=k^{\textrm{o}}+k^{\textrm{no}}r^{2}$, as

$\dot{e}_{p}=\frac{1}{D_{\theta}}\left[(k^{\textrm{o}})^{2}+2k^{\textrm{o}}k^{\textrm{no}}\frac{|k^{\textrm{e}}|}{k^{\textrm{ne}}}+(k^{\textrm{no}})^{2}\left(\frac{|k^{\textrm{e}}|}{k^{\textrm{ne}}}+\frac{2D_{r}}{k^{\textrm{ne}}}\right)\right]$
(64)

for $k^{\textrm{e}}<0$ by using the steady-state distribution and Eqs. (58)
and (59). When $k^{\textrm{e}}\geq 0$, we obtain
$\dot{e}_{p}=(k^{\textrm{o}})^{2}/D_{\theta}$.

These results confirm that the average entropy production rate is generated
only by the odd elasticity, which agrees with the physical interpretation of
the odd elasticity as an internal non-conservative activity. The average work
done by the elastic force is also calculated using

$\dot{W}=-\bigg{\langle}\oint\hat{K}_{\alpha\beta}q_{\beta}\,dq_{\alpha}\bigg{\rangle}=\langle
f^{2}_{\theta}\rangle.$ (65)

We note that the even elasticity is represented by a potential force and does
not contribute to the average work. Comparing Eqs. (63) and (65) shows that
the second law of thermodynamics reads as
$\dot{e}_{p}=\dot{W}/k_{\textrm{B}}T$, with the Boltzmann constant
$k_{\textrm{B}}$. We thus have

$D_{\theta}=k_{\textrm{B}}T,$ (66)

which is consistent with the generalized fluctuating dissipation theorem for
non-equilibrium systems [68].

Of note, the non-equilibrium thermodynamic relations discussed in this
subsection are based on the dynamics in the shape space (46), or more
precisely, the dynamics of Eqs. (5) and (6) with additional noise. The entropy
production rate and the associated heat production for the full system should
account for the non-reciprocal swimming motion in the physical space and the
energy dissipation in the fluid.

Nonetheless, our analysis of the noisy limit cycle in shape space bridges the
entropy production and work done by the odd elasticity in the nonlinear
regime, being consistent with the physical interpretation that actuation
inside the elastic material is described by odd elasticity.

##  Appendix B Probability current with simple additive noise

In this appendix, to complement the results of the biologically relevant
multiplicative noise posited in Section A.1, we discuss the probabilistic
areal velocity for the additive Gaussian noise in Eq. (49). The form of the
noise (49) is given by the setting
$G_{\alpha\beta}=\sqrt{2D}\,\delta_{\alpha\beta}$.

The formal solution (54) then provides the steady probability distribution
[69]

$P_{\textrm{st}}(r)=N_{\textrm{st}}\exp\left[-\frac{k^{\textrm{e}}}{2D}r^{2}-\frac{k^{\textrm{ne}}}{4D}r^{4}\right],$
(67)

where the normalization factor $N_{\textrm{st}}$ is obtained by a direct
integral of
$\int_{-\infty}^{\infty}\int_{-\infty}^{\infty}P_{\textrm{st}}\,dq_{1}dq_{2}=1$
as

$N_{\textrm{st}}=\frac{1}{\pi}\exp\left(-\frac{(k^{\textrm{e}})^{2}}{4Dk^{\textrm{ne}}}\right)\bigg{/}\sqrt{\frac{\pi
D}{k^{\textrm{ne}}}}\,\textrm{erfc}\left(\frac{k^{\textrm{e}}}{2\sqrt{Dk^{\textrm{ne}}}}\right),$
(68)

where $\textrm{erfc}(x)$ is the complementary error function.

The exponential form in Eq. 67 represents the fourth-order potential
associated with the nonlinear odd elasticity, and the steady distribution
function is unimodal when $k^{\textrm{e}}>0$, whereas it becomes crater-shaped
with a degenerate maximum on a circle with radius
$r^{\ast}=(|k^{\textrm{e}}|/k^{\textrm{ne}})^{1/2}$ when $k^{\textrm{e}}<0$.

As in Eq. (60), we can calculate the probabilistic areal velocity via

$J=\frac{k^{\textrm{o}}}{2}\langle
r^{2}\rangle+\frac{k^{\textrm{no}}}{2}\langle
r^{4}\rangle=\mathrel{\mathop{\mathchar
58\relax}}J_{\textrm{o}}+J_{\textrm{no}},$ (69)

which we separated into two terms for subsequent analyses.

After some calculations, we can express the linear odd elastic contribution,
$J_{\textrm{o}}=(k^{\textrm{o}}/2)\langle r^{2}\rangle$, as

$J_{\textrm{o}}=\cfrac{k^{\textrm{o}}}{2}\left(\cfrac{2\pi
D}{k^{\textrm{ne}}}N_{\textrm{st}}-\cfrac{k^{\textrm{e}}}{k^{\textrm{ne}}}\right).$
(70)

Note that $N_{\textrm{st}}$ is a non-negative function of $D$. For example,
when $k^{\textrm{e}}=0$, it becomes
$N_{\textrm{st}}=\sqrt{k^{\textrm{ne}}/[\pi^{3}D]}$. To understand the effects
of noise, we need to examine the behavior of $DN_{\textrm{st}}$, and we found
that this quantity monotonically increases with $D$. Indeed, by the asymptotic
expression for a small $x=4Dk^{\textrm{ne}}/(k^{\textrm{e}})^{2}>0$,

$\frac{\sqrt{x}e^{-1/x}}{\textrm{erfc}(1/\sqrt{x})}=\sqrt{\pi}\left(1+\frac{x}{2}-\frac{x^{2}}{2}\right)+O(x^{3}),$
(71) $\frac{2\pi
D}{k^{\textrm{ne}}}N_{\textrm{st}}-\frac{k^{\textrm{e}}}{k^{\textrm{ne}}}=\frac{k^{\textrm{e}}}{2k^{\textrm{ne}}}\left(x-x^{2}\right)+O(x^{3}),$
(72)

and we obtain an asymptotic behavior for small noise as

$J_{\textrm{o}}=\frac{k^{\textrm{o}}}{k^{\textrm{e}}}D+O(D^{2})\leavevmode\nobreak\
\textrm{when}\leavevmode\nobreak\ k^{\textrm{e}}>0.$ (73)

When $k^{\textrm{e}}<0$, using the asymptotic behavior for small $x>0$,

$\frac{\sqrt{x}e^{-1/x}}{\textrm{erfc}(-1/\sqrt{x})}=\frac{1}{2}e^{-1/x}\left(\sqrt{x}+O(x^{7/2})\right),$
(74)

and we can obtain

$J_{\textrm{o}}\simeq\frac{k^{\textrm{no}}}{2}\frac{|k^{\textrm{e}}|}{k^{\textrm{e}}}\leavevmode\nobreak\
\textrm{when}\leavevmode\nobreak\ k^{\textrm{e}}<0$ (75)

with an exponentially small error.

We found that the average $\langle r^{2}\rangle$ monotonically increases with
$D$, irrespective of the sign of $k^{\textrm{e}}$. The magnitude of the
probabilistic areal velocity is therefore proportional to the odd elasticity
$k^{\textrm{e}}$ and increases as the size of the noise increases, although
the overall sign is determined by the sign of $k^{\textrm{e}}$. This noise
dependence is different from that for the multiplicative noise case, although
the noise effect is exponentially small for little noise.

By similar calculations, we can obtain the probabilistic areal velocity from
the nonlinear odd elasticity, $J_{\textrm{no}}=(k^{\textrm{no}}/2)\langle
r^{4}\rangle$, as

$J_{\textrm{no}}=\frac{k^{\textrm{no}}}{2}\left[\frac{2D}{k^{\textrm{ne}}}-\left(\frac{k^{\textrm{e}}}{k^{\textrm{ne}}}\right)\left(\frac{2\pi
D}{k^{\textrm{ne}}}N_{\textrm{st}}-\frac{k^{\textrm{e}}}{k^{\textrm{ne}}}\right)\right].$
(76)

By a similar asymptotic analysis, by using the expression (71), we obtain
$J_{\textrm{no}}$ for small $D$ when $k^{\textrm{e}}>0$ as

$J_{\textrm{no}}=\frac{2k^{\textrm{no}}}{(k^{\textrm{e}})^{2}}D^{2}+O(D^{3})\leavevmode\nobreak\
\textrm{when}\leavevmode\nobreak\ k^{\textrm{e}}>0,$ (77)

which is positive for small $D$ and monotonically increases (if
$k^{\textrm{no}}>0$) as diffusion is enhanced. When $k^{\textrm{e}}<0$, again
using the asymptotic (74), the probabilistic areal velocity is obtained for
small $D$ as

$J_{\textrm{no}}\simeq\frac{k^{\textrm{no}}}{2}\left[\frac{2D}{k^{\textrm{ne}}}+\left(\frac{k^{\textrm{e}}}{k^{\textrm{ne}}}\right)^{2}\right]\leavevmode\nobreak\
\textrm{when}\leavevmode\nobreak\ k^{\textrm{e}}<0,$ (78)

with an exponentially small error.

## References

  * Bernheim-Groswasser  _et al._ [2018] A. Bernheim-Groswasser, N. S. Gov, S. A. Safran, and S. Tzlil, Living matter: mesoscopic active materials, Advanced Materials 30, 1707028 (2018).
  * Das  _et al._ [2020] M. Das, C. F. Schmidt, and M. Murrell, Introduction to active matter, Soft Matter 16, 7185 (2020).
  * Corbetta and Toschi [2023] A. Corbetta and F. Toschi, Physics of human crowds, Annual Review of Condensed Matter Physics 14 (2023).
  * Shankar  _et al._ [2022] S. Shankar, A. Souslov, M. J. Bowick, M. C. Marchetti, and V. Vitelli, Topological active matter, Nature Reviews Physics 4, 380 (2022).
  * Klapp [2023] S. H. Klapp, Non-reciprocal interaction for living matter, Nature Nanotechnology 18, 8 (2023).
  * Scheibner  _et al._ [2020] C. Scheibner, A. Souslov, D. Banerjee, P. Surówka, W. T. Irvine, and V. Vitelli, Odd elasticity, Nature Physics 16, 475 (2020).
  * Tan  _et al._ [2022] T. H. Tan, A. Mietke, J. Li, Y. Chen, H. Higinbotham, P. J. Foster, S. Gokhale, J. Dunkel, and N. Fakhri, Odd dynamics of living chiral crystals, Nature 607, 287 (2022).
  * Fruchart  _et al._ [2023] M. Fruchart, C. Scheibner, and V. Vitelli, Odd viscosity and odd elasticity, Annual Review of Condensed Matter Physics 14, 471 (2023).
  * Brandenbourger  _et al._ [2021] M. Brandenbourger, C. Scheibner, J. Veenstra, V. Vitelli, and C. Coulais, Limit cycles turn active matter into robots, arXiv preprint arXiv:2108.08837 (2021).
  * Elgeti  _et al._ [2015] J. Elgeti, R. G. Winkler, and G. Gompper, Physics of microswimmers―single particle motion and collective behavior: a review, Reports on progress in physics 78, 056601 (2015).
  * Lauga [2020] E. Lauga,  _The fluid dynamics of cell motility_ , Vol. 62 (Cambridge University Press, 2020).
  * Purcell [1977] E. M. Purcell, Life at low reynolds number, American journal of physics 45, 3 (1977).
  * Ishimoto and Yamada [2012] K. Ishimoto and M. Yamada, A coordinate-based proof of the scallop theorem, SIAM Journal on Applied Mathematics 72, 1686 (2012).
  * Yasuda  _et al._ [2021] K. Yasuda, Y. Hosaka, I. Sou, and S. Komura, Odd microswimmer, Journal of the Physical Society of Japan 90, 075001 (2021).
  * Ishimoto  _et al._ [2022] K. Ishimoto, C. Moreau, and K. Yasuda, Self-organized swimming with odd elasticity, Physical Review E 105, 064603 (2022).
  * Gilpin  _et al._ [2020] W. Gilpin, M. S. Bull, and M. Prakash, The multiscale physics of cilia and flagella, Nature Reviews Physics 2, 74 (2020).
  * Camalet and Jülicher [2000] S. Camalet and F. Jülicher, Generic aspects of axonemal beating, New Journal of Physics 2, 324 (2000).
  * De Canio  _et al._ [2017] G. De Canio, E. Lauga, and R. E. Goldstein, Spontaneous oscillations of elastic filaments induced by molecular motors, Journal of The Royal Society Interface 14, 20170491 (2017).
  * Ling  _et al._ [2018] F. Ling, H. Guo, and E. Kanso, Instability-driven oscillations of elastic microfilaments, Journal of The Royal Society Interface 15, 20180594 (2018).
  * Chakrabarti and Saintillan [2019] B. Chakrabarti and D. Saintillan, Spontaneous oscillations, beating patterns, and hydrodynamics of active microfilaments, Physical Review Fluids 4, 043102 (2019).
  * Man and Kanso [2019] Y. Man and E. Kanso, Morphological transitions of axially-driven microfilaments, Soft Matter 15, 5163 (2019).
  * Ma  _et al._ [2014] R. Ma, G. S. Klindt, I. H. Riedel-Kruse, F. Jülicher, and B. M. Friedrich, Active phase and amplitude fluctuations of flagellar beating, Physical review letters 113, 048101 (2014).
  * Wan and Goldstein [2014] K. Y. Wan and R. E. Goldstein, Rhythmicity, recurrence, and recovery of flagellar beating, Physical review letters 113, 238103 (2014).
  * Saggiorato  _et al._ [2017] G. Saggiorato, L. Alvarez, J. F. Jikeli, U. B. Kaupp, G. Gompper, and J. Elgeti, Human sperm steer with second harmonics of the flagellar beat, Nature communications 8, 1415 (2017).
  * Ishimoto  _et al._ [2017] K. Ishimoto, H. Gadêlha, E. A. Gaffney, D. J. Smith, and J. Kirkman-Brown, Coarse-graining the fluid flow around a human sperm, Physical review letters 118, 124501 (2017).
  * Ishimoto  _et al._ [2018] K. Ishimoto, H. Gadêlha, E. A. Gaffney, D. J. Smith, and J. Kirkman-Brown, Human sperm swimming in a high viscosity mucus analogue, Journal of theoretical biology 446, 1 (2018).
  * Guasto  _et al._ [2020] J. S. Guasto, J. B. Estrada, F. Menolascina, L. J. Burton, M. Patel, C. Franck, A. Hosoi, R. K. Zimmer, and R. Stocker, Flagellar kinematics reveals the role of environment in shaping sperm motility, Journal of the Royal Society Interface 17, 20200525 (2020).
  * Lin  _et al._ [2023] L.-S. Lin, K. Yasuda, K. Ishimoto, Y. Hosaka, and S. Komura, Onsager’s variational principle for nonreciprocal systems with odd elasticity, Journal of the Physical Society of Japan 92, 033001 (2023).
  * Najafi and Golestanian [2004] A. Najafi and R. Golestanian, Simple swimmer at low reynolds number: Three linked spheres, Physical Review E 69, 062901 (2004).
  * Golestanian and Ajdari [2008] R. Golestanian and A. Ajdari, Analytic results for the three-sphere swimmer at low reynolds number, Physical Review E 77, 036308 (2008).
  * Yasuda  _et al._ [2017] K. Yasuda, Y. Hosaka, M. Kuroda, R. Okamoto, and S. Komura, Elastic three-sphere microswimmer in a viscous fluid, Journal of the Physical Society of Japan 86, 093801 (2017).
  * Becker  _et al._ [2003] L. E. Becker, S. A. Koehler, and H. A. Stone, On self-propulsion of micro-machines at low reynolds number: Purcell’s three-link swimmer, Journal of fluid mechanics 490, 15 (2003).
  * Moreau [2019] C. Moreau, Local controllability of a magnetized purcell’s swimmer, IEEE Control Systems Letters 3, 637 (2019).
  * Doi [2013] M. Doi,  _Soft matter physics_ (oxford university press, 2013).
  * Guckenheimer and Holmes [2013] J. Guckenheimer and P. Holmes,  _Nonlinear oscillations, dynamical systems, and bifurcations of vector fields_ , Vol. 42 (Springer Science & Business Media, 2013).
  * Shapere and Wilczek [1989] A. Shapere and F. Wilczek, Geometry of self-propulsion at low reynolds number, Journal of Fluid Mechanics 198, 557 (1989).
  * Koiller  _et al._ [1996] J. Koiller, K. Ehlers, and R. Montgomery, Problems and progress in microswimming, Journal of Nonlinear Science 6, 507 (1996).
  * Avron and Raz [2008] J. E. Avron and O. Raz, A geometric theory of swimming: Purcell’s swimmer and its symmetrized cousin, New Journal of Physics 10, 063016 (2008).
  * Machin [1958] K. Machin, Wave propagation along flagella, Journal of Experimental Biology 35, 796 (1958).
  * Fu  _et al._ [2008] H. C. Fu, C. W. Wolgemuth, and T. R. Powers, Beating patterns of filaments in viscoelastic fluids, Physical Review E 78, 041913 (2008).
  * Bayly and Wilson [2015] P. Bayly and K. Wilson, Analysis of unstable modes distinguishes mathematical models of flagellar motion, Journal of the Royal Society Interface 12, 20150124 (2015).
  * Gadêlha  _et al._ [2010] H. Gadêlha, E. Gaffney, D. Smith, and J. Kirkman-Brown, Nonlinear instability in flagellar dynamics: a novel modulation mechanism in sperm migration?, Journal of The Royal Society Interface 7, 1689 (2010).
  * Moreau  _et al._ [2018] C. Moreau, L. Giraldi, and H. Gadêlha, The asymptotic coarse-graining formulation of slender-rods, bio-filaments and flagella, Journal of the Royal Society Interface 15, 20180235 (2018).
  * Gong  _et al._ [2020] A. Gong, S. Rode, U. B. Kaupp, G. Gompper, J. Elgeti, B. M. Friedrich, and L. Alvarez, The steering gaits of sperm, Philosophical Transactions of the Royal Society B 375, 20190149 (2020).
  * Shiba  _et al._ [2008] K. Shiba, S. A. Baba, T. Inoue, and M. Yoshida, Ca2+ bursts occur around a local minimal concentration of attractant and trigger sperm chemotactic response, Proceedings of the National Academy of Sciences 105, 19312 (2008).
  * Geyer  _et al._ [2016] V. F. Geyer, P. Sartori, B. M. Friedrich, F. Jülicher, and J. Howard, Independent control of the static and dynamic components of the chlamydomonas flagellar beat, Current Biology 26, 1098 (2016).
  * Banerjee  _et al._ [2021] D. Banerjee, V. Vitelli, F. Jülicher, and P. Surówka, Active viscoelasticity of odd materials, Physical Review Letters 126, 138001 (2021).
  * Lier  _et al._ [2022] R. Lier, J. Armas, S. Bo, C. Duclut, F. Jülicher, and P. Surówka, Passive odd viscoelasticity, Physical Review E 105, 054607 (2022).
  * Schoeller  _et al._ [2020] S. F. Schoeller, W. V. Holt, and E. E. Keaveny, Collective dynamics of sperm cells, Philosophical Transactions of the Royal Society B 375, 20190384 (2020).
  * Taketoshi  _et al._ [2023] N. Taketoshi, T. Omori, and T. Ishikawa, Self-organisation of spermatozoa in bulk suspensions via unsteady elastohydrodynamic interactions, arXiv Preprint , 2302.10572 (2023).
  * Ishimoto and Gaffney [2016] K. Ishimoto and E. A. Gaffney, Mechanical tuning of mammalian sperm behaviour by hyperactivation, rheology and substrate adhesion: a numerical exploration, Journal of the Royal Society Interface 13, 20160633 (2016).
  * Simons and Olson [2018] J. E. Simons and S. D. Olson, Sperm motility: Models for dynamic behavior in complex environments, Cell Movement: Modeling and Applications , 169 (2018).
  * Lange and Friedrich [2021] S. Lange and B. M. Friedrich, Sperm chemotaxis in marine species is optimal at physiological flow rates according theory of filament surfing, PLoS computational biology 17, e1008826 (2021).
  * Bennett and Golestanian [2015] R. R. Bennett and R. Golestanian, A steering mechanism for phototaxis in chlamydomonas, Journal of The Royal Society Interface 12, 20141164 (2015).
  * De Maleprade  _et al._ [2020] H. De Maleprade, F. Moisy, T. Ishikawa, and R. E. Goldstein, Motility and phototaxis of gonium, the simplest differentiated colonial alga, Physical Review E 101, 022416 (2020).
  * Kunita  _et al._ [2016] I. Kunita, T. Yamaguchi, A. Tero, M. Akiyama, S. Kuroda, and T. Nakagaki, A ciliate memorizes the geometry of a swimming arena, Journal of the Royal Society interface 13, 20160155 (2016).
  * Ohmura  _et al._ [2018] T. Ohmura, Y. Nishigami, A. Taniguchi, S. Nonaka, J. Manabe, T. Ishikawa, and M. Ichikawa, Simple mechanosense and response of cilia motion reveal the intrinsic habits of ciliates, Proceedings of the National Academy of Sciences 115, 3231 (2018).
  * Echigoya  _et al._ [2022] S. Echigoya, K. Sato, O. Kishida, T. Nakagaki, and Y. Nishigami, Switching of behavioral modes and their modulation by a geometrical cue in the ciliate stentor coeruleus, Frontiers in Cell and Developmental Biology 10 (2022).
  * Risken and Risken [1996] H. Risken and H. Risken,  _Fokker-planck equation_ (Springer, 1996).
  * Tomita and Tomita [1974] K. Tomita and H. Tomita, Irreversible circulation of fluctuation, Progress of theoretical physics 51, 1731 (1974).
  * Wang  _et al._ [2008] J. Wang, L. Xu, and E. Wang, Potential landscape and flux framework of nonequilibrium networks: Robustness, dissipation, and coherence of biochemical oscillations, Proceedings of the National Academy of Sciences 105, 12271 (2008).
  * Fang  _et al._ [2019] X. Fang, K. Kruse, T. Lu, and J. Wang, Nonequilibrium physics in biology, Reviews of Modern Physics 91, 045004 (2019).
  * Yasuda and Komura [2021] K. Yasuda and S. Komura, Nonreciprocality of a micromachine driven by a catalytic chemical reaction, Physical Review E 103, 062113 (2021).
  * Kobayashi  _et al._ [2022] A. Kobayashi, K. Yasuda, K. Ishimoto, L.-S. Lin, I. Sou, Y. Hosaka, and S. Komura, Odd elasticity of a catalytic micromachine, arXiv preprint arXiv:2211.16089 (2022).
  * Graham [1982] R. Graham, Hopf bifurcation with fluctuating control parameter, Physical Review A 25, 3234 (1982).
  * Evans and Searles [2002] D. J. Evans and D. J. Searles, The fluctuation theorem, Advances in Physics 51, 1529 (2002).
  * Feng and Wang [2011] H. Feng and J. Wang, Potential and flux decomposition for dynamical systems and non-equilibrium thermodynamics: Curvature, gauge field, and generalized fluctuation-dissipation theorem, The Journal of chemical physics 135, 234511 (2011).
  * Seifert [2005] U. Seifert, Entropy production along a stochastic trajectory and an integral fluctuation theorem, Physical review letters 95, 040602 (2005).
  * San Miguel and Chaturvedi [1980] M. San Miguel and S. Chaturvedi, Limit cycles and detailed balance in fokker-planck equations, Zeitschrift für Physik B Condensed Matter 40, 167 (1980).