# A computational _medical XR_ discipline

George Papagiannakis
Department of Computer Science
University of Crete
Heraklion, GR 700113
<EMAIL_ADDRESS>also at Foundation for Research and Technology - Hellas,
ORamaVR SA-Switzerland, University of Geneva-Switzerland

###### Abstract

Computational medical XR (extended reality) unifies the computer science
applications of intelligent reality, medical virtual reality, medical
augmented reality and spatial computing for medical training, planning and
navigation content creation. It builds upon clinical XR by bringing on novel
low-code/no-code XR authoring platforms, suitable for medical professionals as
well as XR content creators.

Figure 1: A new computational medical XR (extended reality) concept, where
intelligent reality technologies are employed in the study of unified medical
training, preoperative planning and operative navigation, under integrated
authoring platforms and low-code/no-code metaphors.

_Keywords_ Computational medicine $\cdot$ Clinical XR $\cdot$ XR medicine

##  1 Introduction

This article is an extended version of the one previously published online at
Papagiannakis [2021a].

5 billion people today lack access to proper surgical and anesthesia care and
current training models cannot meet this demand. Furthermore, more than 1
billion jobs, almost one-third of all jobs worldwide, are likely to be
transformed by technology in the next decade, according to OECD estimates.
According to WHO, 40 million healthcare professionals are needed by 2030! This
growing need for continuous upskill and reskill becomes even more crucial in
the post COVID-19 pandemic era. Virtual Reality (VR) together with 5G spatial
computing technologies can pose as the next final computing frontier regarding
medical psychomotor/cognitive training, education, and empowerment.

##  2 VR in training

Virtual/augmented/mixed reality (VR/AR/MR)technologies (grouped by the
industry as XR) are now transforming medical training. Very recent articles
Slater [2020] and case study review articles Slater and Sanchez-Vives [2016]
in industry Dzyuba [2021] as well as dedicated academic special publication
collections Papagiannakis [2021b] highlight this fact. Only in 2020, more than
20 clinical trials were published in key scientific journals, measuring and
testifying the efficacy of medical VR training and skills transfer from
virtual to real, whereas only one in 2019: Hooper et al. [2019] . It has also
been recently reported Research and Innovation [2021] that XR technologies can
offer significant boost in experiential and collaborative learning of
healthcare professionals.

VR/AR can provide the means for remote qualitative education (knowledge) and
training (skills), using affordable technology with personalized, on-demand
and smooth learning curves. Based on recent major advances in the field of 5G
edge computing Kamarianakis et al. [2021], neuroscience Riva et al. [2019],
mixed-reality Zikas et al. [2020] and spatial computing:

_“VR/AR shares with our brain the same basic mechanism: embodied simulations”_
[Riva et al., 2019]

Such immersive technologies can facilitate continuous learning, provide
curriculum programs and self-improvement opportunities that both expand
medical professionals’ abilities as well as minimize skill gaps in their
training.

These spatial computing applications are being based on latest advances on
real-time 3D graphics, computer vision, novel algebraic representations and
simulation computational models.

_"After more than 30 years of intense R &D and intense validation by early
adopters, medical VR technologies are now moving into mainstream"_ [Greenleaf,
2021]

We know from neuroscientists as well as computer scientists and educators that
the feeling of Presence and the affordances of agency (body, hand interaction)
drive neuroplasticity and long-term potentiation, which refers to the
strengthening of connections between neurons in the brain based on recent
patterns of activity. With each repetition, cell-to-cell signaling improves
and the neural connections tighten, profoundly affecting learning, muscle
memory as well as spatial memory retention. Latest advances in affordable VR
hardware coupled with 5G-edge mobile and ubiquitous connectivity with massive
scale/outreach are further driving this proliferation. Therapeutic VR, is also
nowadays a field were the use of VR in the treatment of traumas, phobias,
stress, various mental illnesses as well as in rehabilitation, is now a fact
with clinically proven results Greenleaf [2021], Hassandra et al. [2021].

##  3 VR training in the post-pandemic era

Covid-19 has forever accelerated advances in this field and elevated these
spatial computing technologies from a mere “nice to have” to a “must have” for
the roughly 3000 medical schools worldwide, the 8400 medical device companies
and the 600 surgical training centers, that were all challenged by covid-19
lockdowns. The use of medical VR training is now a new reality and part of the
“quiet revolution”, as termed in a recent Accenture industry report accenture
[2020], where enterprise VR vastly outnumbers consumer VR applications.
Furthermore, healthcare professionals have now acknowledged that training
simulators are the number one immediate application of VR/AR in healthcare
Intelligence [2020]. The forefront of this revolution is the unexpected new
insight that “the rising demand for virtual training across industries drives
the VR market”, a $15Bmarketin2020,expectedtogrowto$57B, as was also recently
reported Research [2021].

However, the looming question is: _who is going to “VRify”(translate into VR)
the 400+ basic medical procedures and their countless variations?_

These very basic medical procedures wikipedia [2021] are constantly being
updated and/or being extended. Furthermore, as medical device equipment for
diagnostic or therapeutic procedures such as robotic surgery get updated, so
too do the variations of these basic medical procedures thereby exponentially
increasing the number of procedures requiring advanced and extensive training.
For healthcare professionals, it’s a constant challenge to upskill / reskill
following latest advances in the field as well as for the institutions to
train enough of them to meet the increasing demand. Every other
week/month/quarter, medical professionals must complete self-training courses
or attend lectures/workshops/conferences to provide the best treatment for
their patients and remain compliant with regulatory agencies. This is
especially true if you’re a med-school student, resident, young or experienced
surgeon or physician of any specialty or nurse. The traditional method of
training for physicians has been to train on cadavers or on-the-job training,
but of course, cadavers are scarce, expensive and logistically limited while
training on real patients poses other issues (limited opportunity,
repeatability not to mention possible adverse outcomes). With capped OR times
for trainees/residents and increasingly fewer of them entering the workforce,
the problem only gets worse over time. Similar needs/problems we also witness
in the adjacent verticals of vocational or corporate skills training.

Despite the above strides, the question of who will finally meet the vast
needs of “VRification” of medical training content, both in extensive depth as
well as exhaustive breadth, remains unanswered. This problem is not just a
“lack of VR content”, as many believe, but mainly a “lack of dedicated VR
content authoring tools” and metaphors [Zikas et al., 2020], specifically for
medical training. The excellent de facto 3D/XR content creation industrial
standards of Unity and Unreal Engines act as the essential enabling
technologies and general authoring platforms. However, their general-purpose
coverage (from games to automotive, engineering, construction to simulation
and film/broadcasting) requires highly skilled developer teams from various
technological disciplines (computer graphics, deep learning, computer
networking, human-computer-interaction, gamification, game development, UI/UX
design, affective computing) coupled with teams of education, medicine,
psychology professionals and ofcourse 3D artists/designers

##  4 The need for Computational Medical XR

I believe that new integrated software platforms and simulation-training
authoring tools are needed Ponder et al. [2003], Zikas et al. [2020]. They
will allow rapid prototyping of medical VR/XR training, as part of new “
computational medical XR” discipline that tackles medical education,
preoperative planning and real-time operative navigation under a new, holistic
integrated systems approach. Medical schools, dental schools, nursing
academies, medical device companies and surgical training centers should
employ teams in-house or outsource parts of the “VRification” of their
curricula but keeping ownership of the final source code that will enable
themselves to easily update/extend/modify this XR training content, hence
driving themselves further the adoption and customization of their medical VR
modules. Through the ability to control and develop their own XR training
material, they can ensure their medical professionals are properly and
(continuously) trained while ensuring optimal patient outcomes and fewer
medical errors/complications. I believe that such new tools will act as the
new ‘WordPress’, or ‘Powerpoint’ of medical VR/XR training, so that eventually
the subject-matter-experts themselves, can create their own VR/AR learning
modules, in the same manner they create their presentation slides today.
Ofcourse there is still work to be done from our side as technologists to
provide these low-code/no-code authoring platforms. An amazing recent
testament to this trend, is the in-house VR content creation department of the
International Committee of the Red Cross and the state-of-the-art content they
are creating ICRC [2021].

The example of aerospace and aviation, which is similar in regulatory
complexity as the healthcare industry, testifies to the unique benefits of
simulation-based training. Pilots must undergo mandatory simulation training
for several types of aircrafts from various aviation authorities worldwide and
this has led to increased safety in aviation worldwide. Several mandates for
simulation-based training have also been put forward in medicine, and it is
only a matter of time till VR simulation education labs for training and
competency assessment will be mandatory on a worldwide basis as well. Medical
education will be significantly boosted via experiential learning, since
learning by doing has proved to be much more effective than the traditional
master-apprentice model and the “see one, do one, teach one” concept.

##  5 Conclusions

For the above transformation to be realized, it is important that medical
institutions drive themselves the medical VR knowledge transformation via XR
content curriculum simulators and continuous education of their professionals,
for the immediate benefit of their patients with significant impact to
humanity in general.

Let’s make the world a better place by advancing XR computational medicine
with novel content authoring platforms.

##  6 Acknowledgements

Parts of this article have been presented in the 2nd Swiss Conference for
Telemedicine and Digital Health 2021, Bern, Switzerland. Special thanks to
Craig Patterson for reviewing this article.

## References

  * Papagiannakis [2021a] George Papagiannakis.  Computational xr medicine, 2021a.  URL https://george-papagiannakis.medium.com/a-new-xr-computational-medicine-discipline-9a48ae1550c.
  * Slater [2020] Mel Slater.  Virtual reality for training, 2020.  URL http://presence-thoughts.blogspot.com/2020/03/virtual-reality-for-training.html.
  * Slater and Sanchez-Vives [2016] Mel Slater and Maria V. Sanchez-Vives.  Enhancing our lives with immersive virtual reality.  _Frontiers in Robotics and AI_ , 3:74, 2016.  ISSN 2296-9144.  doi:10.3389/frobt.2016.00074.  URL https://www.frontiersin.org/article/10.3389/frobt.2016.00074.
  * Dzyuba [2021] Alex Dzyuba.  How immersive technology could enhance healthcare of the post-pandemic reality, 2021.  URL https://tinyurl.com/hrarczdm.
  * Papagiannakis [2021b] George Papagiannakis.  New virtual reality and spatial computing applications to empower, upskill and reskill medical professionals in a post-pandemic era, 2021b.  URL https://tinyurl.com/ys3u7jsy.
  * Hooper et al. [2019] Jessica Hooper, Eleftherios Tsiridis, James E. Feng, Ran Schwarzkopf, Daniel Waren, William J. Long, Lazaros Poultsides, William Macaulay, George Papagiannakis, Eustathios Kenanidis, Eduardo D. Rodriguez, James Slover, Kenneth A. Egol, Donna P. Phillips, Scott Friedlander, and Michael Collins.  Virtual reality simulation facilitates resident training in total hip arthroplasty: A randomized controlled trial.  _The Journal of Arthroplasty_ , 34(10):2278–2283, 2019.  ISSN 0883-5403.  doi:https://doi.org/10.1016/j.arth.2019.04.002.  URL https://www.sciencedirect.com/science/article/pii/S0883540319303341.
  * Research and Innovation [2021] UK Research and Innovation.  The growing value of xr in healthcare in the united kingdom, 2021.  URL https://www.xrhealthuk.org/the-growing-value-of-xr-in-healthcare#Contents.
  * Kamarianakis et al. [2021] Manos N. Kamarianakis, Nick Lydatakis, and George Papagiannakis.  Never ’drop the ball’ in the operating room: An efficient hand-based VR HMD controller interpolation algorithm, for collaborative, networked virtual environments.  _CoRR_ , abs/2107.04875, 2021.  URL https://arxiv.org/abs/2107.04875.
  * Riva et al. [2019] Giuseppe Riva, Brenda Wiederhold, and Fabrizia Mantovani.  Neuroscience of virtual reality: From virtual exposure to embodied medicine.  _Cyberpsychology, Behavior, and Social Networking_ , 22:82–96, 01 2019.  doi:10.1089/cyber.2017.29099.gri.
  * Zikas et al. [2020] Paul Zikas, George Papagiannakis, Nick Lydatakis, Steve Kateros, Stavroula Ntoa, Ilia Adami, and Constantine Stephanidis.  Immersive visual scripting based on VR software design patterns for experiential training.  _Vis. Comput._ , 36(10):1965–1977, 2020.  doi:10.1007/s00371-020-01919-0.  URL https://doi.org/10.1007/s00371-020-01919-0.
  * Greenleaf [2021] Walter Greenleaf.  How virtual reality technology will impact healthcare, 2021.  URL https://www.linkedin.com/pulse/how-virtual-reality-technology-impact-healthcare-greenleaf-phd/.
  * Hassandra et al. [2021] Mary Hassandra, Evangelos Galanis, Antonis Hatzigeorgiadis, Marios Goudas, Christos Mouzakidis, Eleni Maria Karathanasi, Niki Petridou, Magda Tsolaki, Paul Zikas, Giannis Evangelou, George Papagiannakis, George Bellis, Christos Kokkotis, Spyridon Rafail Panagiotopoulos, Giannis Giakas, and Yannis Theodorakis.  $\alpha$ virtual reality app for physical and cognitive training of older people with mild cognitive impairment: Mixed methods feasibility study.  _JMIR Serious Games_ , 9(1):e24170, Mar 2021.  ISSN 2291-9279.  doi:10.2196/24170.  URL https://games.jmir.org/2021/1/e24170.
  * accenture [2020] accenture.  Waking up to a new reality, 2020.  URL https://www.accenture.com/_acnmedia/Accenture/Redesign-Assets/DotCom/Documents/Global/1/Accenture-G20-YEA-report.pdf.
  * Intelligence [2020] Insider Intelligence.  What types of new applications/solutions do us ar/vr professionals expect immersive technologies to offer in the healthcare sector?, 2020.  URL https://tinyurl.com/sc8az5wh.
  * Research [2021] GrandView Research.  Virtual reality market size, share trends analysis report by technology, 2021.  URL https://www.grandviewresearch.com/industry-analysis/virtual-reality-vr-market.
  * wikipedia [2021] wikipedia.  Medical procedure, 2021.  URL https://en.wikipedia.org/wiki/Medical_procedure.
  * Ponder et al. [2003] Michal Ponder, Bruno Herbelin, Tom Molet, Sebastien Schertenlieb, Branislav Ulicny, George Papagiannakis, Nadia Magnenat-Thalmann, and Daniel Thalmann.  Immersive vr decision training: Telling interactive stories featuring advanced virtual human simulation technologies.  In _Proceedings of the Workshop on Virtual Environments 2003_ , EGVE ’03, page 97–106, New York, NY, USA, 2003. Association for Computing Machinery.  ISBN 1581136862.  doi:10.1145/769953.769965.  URL https://doi.org/10.1145/769953.769965.
  * ICRC [2021] ICRC.  Virtual reality innovation, 2021.  URL https://www.icrc.org/en/what-we-do/virtual-reality.


# Optimizing the Mean Swimming Velocity of a Model Two-sphere Swimmer

B. U. Felderhof<EMAIL_ADDRESS>Institut für Theorie der
Statistischen Physik
RWTH Aachen University
Templergraben 55
52056 Aachen
Germany

###### Abstract

The swimming of a two-sphere system oscillating in a viscous fluid is studied
on the basis of simplified equations of motion which take account of both
friction and inertial effects. In the model the friction follows from an Oseen
approximation to the mobility matrix, and the inertial effects follow from a
dipole approximation to the added mass matrix. The resulting mean swimming
velocity is evaluated analytically in a first harmonics approximation. For
specific choices of the parameters this is compared with the exact result
following from a numerical calculation including higher harmonics. The Oseen-
Dipole model is compared with the simpler Oseen* model, in which the added
mass effects are approximated by just the effective mass of the single spheres
and dipole interactions are neglected. The expression for the mean swimming
velocity can be reduced to a dimensionless scaling form. For given viscosity
and mass density of the fluid the frequency of the stroke and the ratio of
radii can be chosen such that the swimming velocity is optimized.

##  I Introduction

It was shown some years ago in experiment and computer simulation by Klotsa et
al. [2] that a system of two neutrally buoyant spheres, oscillating relative
to each other along the axis connecting the centers, will swim when immersed
in a viscous fluid, provided that the two radii differ. Subsequently we showed
that for a class of linear chain models with chosen hydrodynamic interactions,
the mean swimming velocity can be evaluated from an exact expression [3].
Unfortunately, in the expression the dependence on parameters is not explicit,
and this makes it hard to get qualitative insight. In the following we
consider two two-sphere models in more detail, the Oseen∗ model and the Oseen-
Dipole model.

In the Oseen model frictional hydrodynamic interactions are treated in Oseen
approximation, and added mass effects are neglected. In the Oseen∗ model the
mass of the spheres is replaced by the effective mass of the single spheres.
In the Oseen-Dipole model, in addition to the Oseen interaction, added mass
effects are treated in dipole approximation. We shall show that a first
harmonics approximation is accurate and provides qualitative information.

Due to the scallop theorem [4] the mean swimming velocity vanishes in the
Stokes limit of high viscosity in all three models. The behavior beyond this
limit in the Oseen model was studied recently by Hubert et al. [5] to second
order in the amplitude of the stroke. We commented that the mean swimming
velocity can be evaluated in a first harmonics approximation in close
agreement with the exact result, and that when expanded to second order in the
amplitude and for large center-to-center distance, this agrees with the second
order expression [6]. Hubert et al. [5] also performed experiments and lattice
Boltzmann simulations, and showed good agreement with their small amplitude
theory.

In the following we investigate the effect of inertia on the swimming of a
two-sphere on the basis of the Oseen∗ model and the Oseen-Dipole model. We
show that the analytic expression for the mean swimming velocity derived in
first harmonics approximation is practically useful for qualitative
exploration.In first harmonics approximation only zeroth and first harmonics
are taken into account. The results can be compared with the exact ones,
obtained by taking account of higher harmonics in a numerical calculation.

We consider two spheres of radii $a$ and $b$ immersed in a viscous
incompressible fluid of shear viscosity $\eta$ and mass density $\rho$. The
spheres are assumed to be uniform with mass densities $\rho_{a}$ and
$\rho_{b}$. We consider situations where the sphere centers perform harmonic
oscillations along the center-to-center line, which is taken to be the $x$
axis of a Cartesian system of coordinates. It suffices to consider the
coordinates $x_{1}(t)$ and $x_{2}(t)$ of the two centers. The model
calculations are based on postulated equations of motion [7] for the two
center positions ${\bf\sf R}=(x_{1},x_{2})$ and velocities ${\bf\sf
U}=(U_{1},U_{2})$. The equations involve interaction forces, actuating forces,
and kinetic forces deriving from the position dependence of the kinetic energy
of flow [8], as expressed in the mass matrix. We take a kinematic point of
view and assume that the relative position $x(t)=x_{2}(t)-x_{1}(t)$ is
prescribed as

$x(t)=d+f\sin(\omega t).$ (1)

This suffices to determine the corresponding periodic solution of the
equations of motion. If desired, the various forces can eventually be
calculated from the solution.

We are interested in the asymptotic periodic motion of the center
$C(t)=(x_{1}(t)+x_{2}(t))/2$. The average of the velocity $U(t)=dC/dt$ over a
period at long times defines the mean swimming velocity $\overline{U_{sw}}$.
We write this as

$\overline{U_{sw}}=Va\omega,$ (2)

with a dimensionless factor $V$, which is called the velocity function. We
note that $2\pi V$ is the fraction of $a$ traversed in a period of time
$T=2\pi/\omega$. We show that the first harmonics approximation $V^{(1)}$ is a
function of the length ratios $\xi=b/a,\;\delta=d/a,\;\varepsilon=f/a$, and
the mass density ratios $\rho_{b}/\rho_{a}$ and $\rho/\rho_{a}$. It depends in
addition on stroke frequency $\omega$, viscosity $\eta$, and mass density
$\rho$ via the dimensionless combination $R=a^{2}\omega\rho/\eta$. We focus in
particular on the question of how to optimize $V^{(1)}$ by variation of
parameters.

We expect‘that the Oseen-Dipole model catches the main features of the
physical two-sphere swimmer. It provides an approximate account of both
friction and added mass effects. However, the numerical solution of the
Navier-Stokes equations by Dombrowski and Klotsa [9] showed that at high
frequency the system shows a transition with a reversal of the direction of
swimming. Recently Derr et al. [10] analyzed the flow about an oscillating
dimer at small amplitudes on the basis of the Navier-Stokes equations. They
found a reversal of swimming velocity at high frequency due to a boundary
layer effect, similar to that found by Riley for steady streaming about a
single sphere [11]. Hubert et al. [5] did not see the transition in their
lattice Boltzmann simulations. It would clearly be of interest to modify the
hydrodynamic interactions of our mechanical model in such a way as to
reproduce the reversal of velocity. Such an investigation is beyond the scope
of the present work.

##  II Oseen-Dipole two-sphere swimmer

The so-called Oseen-Dipole two-sphere swimmer is a simple model of a physical
two-sphere swimmer consisting of two spheres immersed in an incompressible
viscous fluid oscillating relative to each other along a common axis. In the
model the hydrodynamic interactions between the two spheres are approximated
by instantaneous point interactions calculated from an Oseen approximation to
the $2\times 2$ mobility matrix $\mu$ and a dipole approximation to the
$2\times 2$ mass matrix ${\bf\sf m}$. Higher order multipoles and forces of
the Basset type are neglected.

The inverse of the Oseen mobility matrix yields a corresponding approximation
to the friction matrix $\zeta$. Explicitly the friction matrix is given by

$\mbox{\boldmath$\zeta$}=\frac{12\pi\eta
x}{4x^{2}-9ab}\left(\begin{array}[]{cc}2ax&-3ab\\\
-3ab&2bx\end{array}\right).$ (3)

Here the time-dependent distance $x$ is given by Eq. (1.1).

The approximation to the mass matrix reads

${\bf\sf
m}=\frac{2\pi}{3x^{6}-3a^{3}b^{3}}\left(\begin{array}[]{cc}a^{3}[x^{6}(\rho+2\rho_{a})+2a^{3}b^{3}(\rho-\rho_{a})]&-3a^{3}b^{3}x^{3}\rho\\\
-3a^{3}b^{3}x^{3}\rho&b^{3}[x^{6}(\rho+2\rho_{b})+2a^{3}b^{3}(\rho-\rho_{b})]\end{array}\right).$
(4)

This also depends on time via Eq. (1.1). At large distance the matrix tends to
the time-independent form

${\bf\sf m}_{\infty}=\left(\begin{array}[]{cc}m_{1}^{*}&0\\\
0&m_{2}^{*}\end{array}\right)$ (5)

with effective masses

$m_{1}^{*}=\frac{2\pi}{3}\;a^{3}(2\rho_{a}+\rho),\qquad
m_{2}^{*}=\frac{2\pi}{3}\;b^{3}(2\rho_{b}+\rho).$ (6)

Here the second term in each case is called the added mass. This represents
the kinetic energy of dipolar flow about a moving sphere [8].

The center velocity $U(t)=dC/dt$ satisfies the equation of motion [7]

$\frac{d}{dt}\;(MU)+ZU=\mathcal{I}$ (7)

with time-dependent mass $M(t)$ and friction coefficient $Z(t)$ given by

$M={\bf\sf u}\cdot{\bf\sf m}\cdot{\bf\sf u},\qquad Z={\bf\sf
u}\cdot\mbox{\boldmath$\zeta$}\cdot{\bf\sf u},\qquad{\bf\sf u}=(1,1),$ (8)

and impetus $\mathcal{I}(t)$ given by

$\mathcal{I}(t)=-\frac{d}{dt}({\bf\sf u}\cdot{\bf\sf m}\cdot\dot{{\bf\sf
d}})-{\bf\sf u}\cdot\mbox{\boldmath$\zeta$}\cdot\dot{{\bf\sf d}},$ (9)

where $\dot{{\bf\sf d}}$ is the time-derivative of the displacement vector
${\bf\sf d}(t)$,

${\bf\sf d}=\frac{1}{2}(-x+d,x-d).$ (10)

The time-dependent mass is found to be

$M(t)=m_{1}^{*}+m_{2}^{*}-2\pi\rho
a^{3}b^{3}\frac{2x^{3}-a^{3}-b^{3}}{x^{6}-a^{3}b^{3}},$ (11)

with $x(t)$ given by Eq. (1.1). The friction coefficient is found as

$Z(t)=24\pi\eta x\frac{(a+b)x-3ab}{4x^{2}-9ab}.$ (12)

The impetus $\mathcal{I}(t)$ is a sum of two terms,

$\mathcal{I}(t)=\mathcal{I}_{m}(t)+\mathcal{I}_{f}(t),$ (13)

with mass contribution $\mathcal{I}_{m}(t)$ and friction contribution
$\mathcal{I}_{f}(t)$. The mass contribution is written conveniently in the
form

$\mathcal{I}_{m}(t)=\frac{dG_{m}}{dt},$ (14)

with $G_{m}(t)$ given by

$G_{m}(t)=\frac{1}{2}f\omega(m_{1}^{*}-m_{2}^{*})\cos(\omega t)-\pi
f\rho\omega\frac{a^{3}b^{3}(a^{3}-b^{3})}{a^{3}b^{3}-x^{6}}\cos(\omega t).$
(15)

The friction contribution is

$\mathcal{I}_{f}(t)=12\pi\eta\omega(a-b)f\cos(\omega
t)\frac{x^{2}}{4x^{2}-9ab}.$ (16)

We note that this can be expressed as

$\mathcal{I}_{f}(t)=\frac{dG_{f}}{dt}$ (17)

with periodic function $G_{f}(t)$. Explicitly

$G_{f}(t)=3\pi\eta(a-b)\bigg{[}x+\frac{3}{4}\sqrt{ab}\log\frac{(3\sqrt{ab}+2d)(3\sqrt{ab}-2x)}{(3\sqrt{ab}-2d)(3\sqrt{ab}+2x)}\bigg{]}.$
(18)

It follows from Eqs. (2.12) and (2.15) that the time average of
$\mathcal{I}(t)$ over a period vanishes. It follows from this property and
from Eq. (2.5) that for the periodic solution $U(t)$ the time-average of the
drag force $Z(t)U(t)$ over a period vanishes: $\overline{ZU}=0$, in agreement
with general theory [7]. In the following we discuss the solution of Eq. (2.5)
and the corresponding calculation of the mean swimming velocity
$\overline{U_{sw}}$.

##  III First harmonics approximation

In earlier work [3] we derived an integral expression for the swim velocity
$U_{sw}(t)$ corresponding to the asymptotic periodic solution of the swim
equation Eq. (2.5). Hence we derived an expression for the mean swimming
velocity $\overline{U_{sw}}$, the average of $U_{sw}(t)$ over a period
$T=2\pi/\omega$.

The expressions allow exact calculation, but do not provide easy information
on the dependence on parameters. Here we describe an alternative method which
provides more insight and leads to a faster numerical scheme. In the new
scheme we solve Eq. (2.5) by expansion in harmonics. The periodic swim
velocity $U_{sw}(t)$ is expanded as

$U_{sw}(t)=U_{0}+\sum^{\infty}_{n=1}\big{[}U_{nc}\cos(n\omega
t)+U_{ns}\sin(n\omega t)\big{]}.$ (19)

Clearly $\overline{U_{sw}}=U_{0}$. In $n$th harmonics approximation we include
harmonics of orders 0 to $n$, and ignore higher orders. The calculation
converges rapidly with increasing $n$. Moreover, it turns out that already the
first harmonics approximation with $n=1$ is numerically close to the exact
result.

We examine the first harmonics approximation in analytic form. In this
approximation we have

$U_{sw}^{(1)}(t)=U_{0}^{(1)}+U_{1c}^{(1)}\cos(\omega
t)+U_{1s}^{(1)}\sin(\omega t).$ (20)

Substituting into Eq. (2.5) and expanding the time-dependent factors as in Eq.
(3.2) we obtain a set of three linear equations for the three coefficients
$U_{0}^{(1)},\;U_{1c}^{(1)},\;U_{1s}^{(1)}$ reading

$\displaystyle 2Z_{0}U_{0}^{(1)}+Z_{1c}U_{1c}^{(1)}+Z_{1s}U_{1s}^{(1)}$
$\displaystyle=$ $\displaystyle 2\mathcal{I}_{0},$
$\displaystyle(Z_{1c}+\omega M_{1s})U_{0}^{(1)}+Z_{0}U_{1c}^{(1)}+\omega
M_{0}U_{1s}^{(1)}$ $\displaystyle=$ $\displaystyle\mathcal{I}_{1c},$
$\displaystyle(Z_{1s}-\omega M_{1c}U_{0}^{(1)}-\omega
M_{0}U_{1c}^{(1)}+Z_{0}U_{1s}^{(1)}$ $\displaystyle=$
$\displaystyle\mathcal{I}_{1s},$ (21)

with coefficients which can be calculated from the above. The equations
simplify by use of the properties

$\mathcal{I}_{0}=0,\qquad Z_{1c}=0,\qquad M_{1c}=0.$ (22)

From the solution of Eq. (3.3) we find in particular

$\overline{U_{sw}}^{(1)}=\frac{-Z_{1s}(Z_{0}\mathcal{I}_{1s}+\omega
M_{0}\mathcal{I}_{1c})}{2Z_{0}^{3}-Z_{0}Z_{1s}^{2}+2\omega^{2}M_{0}^{2}Z_{0}-\omega^{2}M_{0}M_{1s}Z_{1s}}.$
(23)

We also find the solution for $U_{1c}^{(1)}$ and $U_{1s}^{(1)}$.

For the mean friction coefficient $Z_{0}$ we find

$Z_{0}=3\pi\eta\bigg{[}2(a+b)-3\alpha\beta\frac{(\alpha+\beta)^{2}}{q_{+}}+3\alpha\beta\frac{(\alpha-\beta)^{2}}{q_{-}}\bigg{]},$
(24)

with the abbreviations $\alpha=\sqrt{a},\;\beta=\sqrt{b}$, and

$q_{\pm}=\sqrt{4d^{2}-4f^{2}+9ab\pm 12d\sqrt{ab}}.$ (25)

For the coefficient $Z_{1s}$ we find

$Z_{1s}=\frac{9\pi\eta\alpha\beta}{f}\bigg{[}\frac{(\alpha+\beta)^{2}}{q_{+}}(2d+3\alpha\beta)-\frac{(\alpha-\beta)^{2}}{q_{-}}(2d-3\alpha\beta)-4\alpha\beta\bigg{]}.$
(26)

The coefficients $\mathcal{I}_{1c}$ and $\mathcal{I}_{1s}$ can be expressed as
a sum of mass and friction contributions

$\mathcal{I}_{1c}=\mathcal{I}_{1cm}+\mathcal{I}_{1cf},\qquad\mathcal{I}_{1s}=\mathcal{I}_{1sm}+\mathcal{I}_{1sf},$
(27)

corresponding to the two terms in Eq. (2.7). The coefficients
$\mathcal{I}_{1cm}$ and $\mathcal{I}_{1sf}$ vanish. The coefficient
$\mathcal{I}_{1c}$ is given by

$\mathcal{I}_{1c}=\frac{3\pi\eta\omega}{4f}(a-b)\big{[}4f^{2}-18ab+3\sqrt{ab}(q_{+}-q_{-})\big{]}.$
(28)

The coefficients $\mathcal{I}_{1sm}$, $M_{0}$ and $M_{1s}$ must be calculated
from integrals over a period after substitution of $x(t)$ from Eq. (1.1). We
evaluate the required integrals by regarding $x$ as a complex variable and
decomposing the expression into partial fractions. With the definition

$w_{j}=\exp(2j\pi i/6),\qquad j=(1,...,6)$ (29)

we have

$\frac{1}{g^{6}-x^{6}}=\frac{1}{6g^{5}}\sum^{6}_{j=1}\frac{w_{j}}{w_{j}g-x}.$
(30)

Correspondingly we find in particular for the coefficient $\mathcal{I}_{1sm}$

$\mathcal{I}_{1sm}=-\frac{1}{2}M_{d}\omega^{2}f+J_{1sm},$ (31)

with $M_{d}=m_{1}^{*}-m_{2}^{*}$ and last term

$J_{1sm}=\frac{\pi}{3}\rho\omega^{2}\;\alpha\beta(\alpha^{6}-\beta^{6})\sum^{6}_{j=1}w_{j}Q(f,d+w_{j}\alpha\beta),$
(32)

with the function

$Q(f,x)=\frac{1}{f}\big{[}x-\sqrt{x^{2}-f^{2}}\big{]}.$ (33)

We note that $Q(f,x)\approx f/(2x)$ for $f<<x$, and that it is odd in $f$.

By the same method we find for the mean total mass

$M_{0}=m_{1}^{*}+m_{2}^{*}+M_{0int},$ (34)

where $M_{0int}$ is due to hydrodynamic interactions between the spheres. We
find

$M_{0int}=\frac{-\pi}{3}\;\rho\alpha\beta\sum^{6}_{j=1}\big{(}\alpha^{3}+(-1)^{j}\beta^{3}\big{)}^{2}\;\frac{w_{j}}{W_{j}},$
(35)

with

$W_{j}=\big{[}d^{2}-f^{2}+2w_{j}d\alpha\beta+w_{j}^{2}\alpha^{2}\beta^{2}\big{]}^{1/2}.$
(36)

Similarly

$\displaystyle
M_{1s}=\frac{2\pi}{3f}\rho\alpha\beta\sum^{6}_{j=1}\big{(}\alpha^{3}+(-1)^{j}\beta^{3}\big{)}^{2}\;\frac{w_{j}}{W_{j}}\big{(}d+w_{j}\alpha\beta\big{)}.$
(37)

We note that the coefficients $M_{0int},M_{1s}$ are independent of the mass
densities $\rho_{a},\rho_{b}$ of the spheres, and the dependence of
$\mathcal{I}_{1s}$ is simple. This can be used to advantage in the
optimization of the mean swimming velocity.

##  IV Scaling in first harmonics approximation

The added mass effect significantly modifies the swimming performance for
given stroke. In order to get some qualitative understanding it is worthwhile
to consider the first harmonics approximation. In Eq. (3.5) for the mean
swimming velocity we put

$\displaystyle Z_{0}$ $\displaystyle=$ $\displaystyle\eta az_{0},\qquad
Z_{1s}=\eta az_{1s},\qquad M_{0}=\rho a^{3}m_{0},$ $\displaystyle M_{d}$
$\displaystyle=$ $\displaystyle\rho a^{3}m_{d},\qquad M_{1s}=\rho
a^{3}m_{1s},\qquad M_{0int}=\rho a^{3}m_{0i},$
$\displaystyle\mathcal{I}_{1cf}$ $\displaystyle=$ $\displaystyle\eta\omega
a^{2}j_{1cf},\qquad\mathcal{I}_{1sm}=-\frac{1}{2}f(m_{1}^{*}-m_{2}^{*})\omega^{2}+R\eta\omega
a^{2}j_{1sm},$ (38)

where $R=a^{2}\omega\rho/\eta$. By comparison with Eq. (3.13)

$j_{1sm}=\frac{\pi}{3}\frac{\alpha\beta(\alpha^{6}-\beta^{6})}{a^{4}}\sum^{6}_{j=1}w_{j}Q(f,d+w_{j}\alpha\beta).$
(39)

The lower case coefficients $z_{0},\;z_{1s},\;j_{1cf},\;j_{1sm},\;m_{1s}$ are
complicated dimensionless functions of the length ratios $b/a,\;d/a,\;f/a$.
The mass coefficient $m_{0}$ involves the single sphere bare and added masses,
as well as the added mass corresponding to $M_{0int}$. The dimensionless
velocity function, defined by $V=\overline{U_{sw}}/(a\omega)$, is given in
first harmonics approximation by

$V^{(1)}=\frac{1}{3}z_{1s}\;\frac{\pi\varepsilon
z_{0}(1-\xi^{3}+2\sigma_{a}-2\xi^{3}\sigma_{b})-3z_{0}j_{1sm}-3m_{0}j_{1cf}}{z_{0}(2z_{0}^{2}-z_{1s}^{2})+m_{0}(2m_{0}z_{0}-m_{1s}z_{1s})R^{2}}\;R,$
(40)

with dimensionless coefficients

$\varepsilon=f/a,\qquad\xi=b/a,\qquad\sigma_{a}=\rho_{a}/\rho=1/S_{a},\qquad\sigma_{b}=\rho_{b}/\rho.$
(41)

The velocity function $V^{(1)}$ depends in simple fashion on viscosity,
frequency and fluid mass density via the dimensionless scaling variable
$R=a^{2}\omega\rho/\eta$. The dependence on $\sigma_{a},\sigma_{b}$ is also
simple. The dependence on the four lengths $a,b,d,f$ is complicated. The
expression Eq. (4.3) shows that the mean swimming velocity is generated by an
intricate interplay of the effects of friction, mass, and impetus.

We write Eq. (4.3) in the abbreviated form

$V^{(1)}=\frac{AR}{B+CR^{2}}.$ (42)

As a function of the number $R$ this takes its maximum value $V_{x}^{(1)}$ at

$R^{(1)}_{x}=\sqrt{\frac{B}{C}},\qquad V^{(1)}_{x}=\frac{A}{2B}R^{(1)}_{x}.$
(43)

The half-width of the curve $V^{(1)}(R)$ is

$\Delta V^{(1)}=2\sqrt{3}R^{(1)}_{x}.$ (44)

We can rewrite Eq. (4.5) as

$\qquad V^{(1)}=2\frac{R^{(1)}_{x}V^{(1)}_{x}R}{{R^{(1)}_{x}}^{2}+R^{2}}.$
(45)

The two ratios $A/B$ and $A/C$ correspond to the behavior of the mean swimming
velocity $\overline{U_{sw}}$ of a two-sphere as a function of frequency
$\omega$ in a fluid of kinematic viscosity $\eta/\rho$ according to

$\displaystyle\overline{U_{sw}}$ $\displaystyle\approx$ $\displaystyle
2\frac{a^{3}\omega^{2}\rho}{\eta}\;\frac{V^{(1)}_{x}}{R^{(1)}_{x}}\;\;\mathrm{as}\;\;\omega\rightarrow
0,$ $\displaystyle\overline{U_{sw}}$ $\displaystyle\approx$ $\displaystyle
2\frac{\eta}{a\rho}R^{(1)}_{x}V^{(1)}_{x}\;\;\mathrm{as}\;\;\omega\rightarrow\infty.$
(46)

It is interesting to observe that the coefficient $B$ depends only on friction
coefficients, that the coefficient $C$ arises from a combination of mass and
friction, and that the coefficient $A$ involves a combination of friction,
mass, and impetus. The first term in the numerator of the expression in Eq.
(4.3) is a measure of the asymmetry of the two-sphere. In addition there are
the two terms with moments of the impetus. The precise combination of terms in
Eq. (4.3) follows from the details of the equations of motion. It would be
difficult to write down the expression on intuitive grounds.

We comment on the simple frequency-dependence found in Eq. (4.3). Consider the
linear susceptibility $\chi(\omega)$ corresponding to a single exponential
relaxation process,

$\chi(\omega)=\chi_{0}\int^{\infty}_{0}e^{i\omega t}\gamma e^{-\gamma t}dt.$
(47)

This is given by

$\chi(\omega)=\chi_{0}\frac{\gamma}{\gamma-i\omega}.$ (48)

The absorption is given by the imaginary part

$\chi^{\prime\prime}(\omega)=\chi_{0}\frac{\gamma\omega}{\gamma^{2}+\omega^{2}}.$
(49)

This is precisely the frequency-dependence seen in Eq. (4.5). The mathematical
equivalence does not allow us to conclude that there is a relaxation process
involved in the mechanism of swimming. However, it does suggest how to modify
the hydromechanical model to get agreement with the behavior seen in the
numerical calculations by Dombrowski and Klotsa [9].

We can generalize Eq. (4.11) by adding a second pole on the imaginary
$\omega$-axis to get the susceptibility

$\chi_{2}(\omega)=\chi_{0}\frac{\gamma}{\gamma-i\omega}+\chi_{02}\frac{\gamma_{2}}{\gamma_{2}-i\omega}.$
(50)

In Fig. 1 we plot $\chi^{\prime\prime}_{2}(\omega)$ and compare with
$\chi^{\prime\prime}(\omega)$ as functions of $\omega$ for the numerical
example $\chi_{0}=1,\;\gamma=1,\;\chi_{02}=-1,\;\gamma_{2}=10$. This shows
that the addition of the second pole with negative residue leads to the type
of behavior seen by Dombrowski and Klotsa [9]. The hydrodynamic interactions
in our hydromechanical model must be modified to reproduce this behavior in
the high frequency regime.

##  V Expansion in amplitude

The dimensionless amplitude of the stroke is defined as the ratio
$\varepsilon=f/a$. Since the mean swimming velocity cannot depend on the phase
of the stroke the velocity function $V$ must be an even function of
$\varepsilon$. The first few terms of the Taylor expansion read

$V=v_{2}\varepsilon^{2}+v_{4}\varepsilon^{4}+O(\varepsilon^{6}).$ (51)

However, in first harmonic approximation it is clearly advantageous to keep
the analytic structure in the variable $R$, as given by Eq. (4.5). The number
$R=a^{2}\omega\rho/\eta$ contains the properties of the fluid, whereas the
other variables refer to the spheres. In the following we derive expressions
for the first two terms in the expansion of the various coefficients in Eq.
(4.3) in powers of the amplitude $\varepsilon$.

First we consider the friction coefficient. From Eq. (3.6) we find

$z_{0}=z_{0,0}+\varepsilon^{2}z_{0,2}+O(\varepsilon^{4}),$ (52)

with first term

$z_{0,0}=6\pi\frac{a+b}{a}-18\pi\frac{b(4d-3a-3b)}{4d^{2}-9ab},$ (53)

and second term

$z_{0,2}=36\pi\frac{a^{2}b}{(4d^{2}-9ab)^{3}}\big{[}(36d^{2}+27ab)(a+b)-108abd-16d^{3}\big{]}.$
(54)

For the friction moment $z_{1s}$ we find

$z_{1s}=\varepsilon z_{1s,1}+\varepsilon^{3}z_{1s,3}+O(\varepsilon^{5}),$ (55)

with coefficients

$\displaystyle z_{1s,1}$ $\displaystyle=$ $\displaystyle 72\pi
ab\frac{(2d-3a)(2d-3b)}{(4d^{2}-9ab)^{2}},$ $\displaystyle z_{1s,3}$
$\displaystyle=$ $\displaystyle
216\pi\frac{a^{3}b}{(4d^{2}-9ab)^{4}}\big{[}81a^{2}b^{2}-108abd(a+b)+216abd^{2}-48(a+b)d^{3}+16d^{4}\big{]}.$

Next we consider the mass. For the added mass in Eq. (3.17) we find

$m_{0int}=m_{0i,0}+\varepsilon^{2}m_{0i,2}+O(\varepsilon^{4}),$ (57)

with the coefficient

$m_{0i,0}=-2\pi\frac{b^{3}(a^{3}+b^{3}-2d^{3})}{a^{3}b^{3}-d^{6}},$ (58)

and with

$m_{0i,2}=3\pi\frac{a^{2}b^{3}d}{(a^{3}b^{3}-d^{6})^{3}}\big{[}2a^{6}b^{6}-5a^{3}b^{3}(a^{3}+b^{3})d^{3}+18a^{3}b^{3}d^{6}-7(a^{3}+b^{3})d^{9}+4d^{12}\big{]}.$
(59)

For the total mass we have

$M_{0,0}=m_{1}+\frac{1}{2}m_{1f}+m_{2}+\frac{1}{2}m_{2f}+\rho
a^{3}m_{0i,0},\qquad m_{0,2}=m_{0i,2}.$ (60)

For the moment $M_{1s}=\rho a^{3}m_{1s}$ in Eq. (3.19) we find the expansion

$m_{1s}=m_{1s,1}\varepsilon+m_{1s,3}\varepsilon^{3}+O(\varepsilon^{5}),$ (61)

with coefficient of $\varepsilon$ given by

$m_{1s,1}=\frac{\pi}{3a^{2}}\alpha\beta\big{[}(\alpha^{6}+\beta^{6})\mu_{1}+2\alpha^{3}\beta^{3}\mu_{2}\big{]},$
(62)

with factors

$\displaystyle\mu_{1}$ $\displaystyle=$
$\displaystyle\sum^{6}_{j=1}\big{[}w_{j}(d-(-1)^{j}w_{j}^{2}\alpha\beta)^{2}\big{]}^{-1}=-36\frac{d^{5}\alpha^{5}\beta^{5}}{(d^{6}-a^{3}b^{3})^{2}},$
$\displaystyle\mu_{2}$ $\displaystyle=$
$\displaystyle\sum^{6}_{j=1}\big{(}-1)^{j}[w_{j}(d-(-1)^{j}w_{j}^{2}\alpha\beta)^{2}\big{]}^{-1}=18d^{2}ab\frac{d^{6}+a^{3}b^{3}}{(d^{6}-a^{3}b^{3})^{2}}.$
(63)

The coefficient of $\varepsilon^{3}$ in Eq. (5.11) is given by

$m_{1s,3}=\frac{\pi}{3}\alpha\beta\big{[}(\alpha^{6}+\beta^{6})\mu_{3}+2\alpha^{3}\beta^{3}\mu_{4}\big{]},$
(64)

with factors

$\displaystyle\mu_{3}$ $\displaystyle=$
$\displaystyle\frac{3}{4}\sum^{6}_{j=1}\big{[}-w_{j}(d-(-1)^{j}w_{j}^{2}\alpha\beta)^{4}\big{]}^{-1}$
$\displaystyle=$ $\displaystyle-18d^{3}\alpha\beta
a^{2}b^{2}\frac{14d^{12}+35d^{6}a^{3}b^{3}+5a^{6}b^{6}}{(d^{6}-a^{3}b^{3})^{4}},$
$\displaystyle\mu_{4}$ $\displaystyle=$
$\displaystyle\frac{3}{4}\sum^{6}_{j=1}(-1)^{j}\big{[}w_{j}(d-(-1)^{j}w_{j}^{2}\alpha\beta)^{4}\big{]}^{-1}$
(65) $\displaystyle=$
$\displaystyle\frac{9}{2}ab\big{[}10d^{18}+125d^{12}a^{3}b^{3}+80d^{6}a^{6}b^{6}+a^{9}b^{9}\big{]}\big{/}(d^{6}-a^{3}b^{3})^{4}.$

Finally we consider the impetus. Expansion of the function $j_{1sm}$ in Eq.
(4.2) in powers of $\varepsilon$ yields

$j_{1sm}=j_{1sm,1}\varepsilon+j_{1sm,3}\varepsilon^{3}+O(\varepsilon^{5}),$
(66)

with coefficients

$\displaystyle j_{1sm,1}$ $\displaystyle=$
$\displaystyle\pi\frac{b^{3}(a^{3}-b^{3})}{a^{3}b^{3}-d^{6}},$ $\displaystyle
j_{1sm,3}$ $\displaystyle=$
$\displaystyle\frac{3\pi}{4}\frac{a^{2}b^{3}d^{4}(a^{3}-b^{3})(5a^{3}b^{3}+7d^{6})}{(a^{3}b^{3}-d^{6})^{3}}.$
(67)

Similarly we find for the function $\mathcal{I}_{1c}$ in Eq. (2.18)

$j_{1c}=j_{1c,1}\varepsilon+j_{1c,3}\varepsilon^{3}+O(\varepsilon^{5}),$ (68)

with the expressions

$\displaystyle j_{1c,1}$ $\displaystyle=$ $\displaystyle
12\pi\frac{d^{2}(a-b)}{a(4d^{2}-9ab)},$ $\displaystyle j_{1c,3}$
$\displaystyle=$ $\displaystyle
81\pi\frac{a^{2}b(a-b)(3ab+4d^{2})}{(4d^{2}-9ab)^{3}}.$ (69)

Substituting the above expansions into Eq. (4.3) we can compare with the
complete result and with the first term in Eq. (5.1). As an example we
consider the case
$\eta=0.5,\;\rho=1,\;a=1,\;b=0.5,\;d=3.5,\;\rho_{a}=\rho_{b}=1,\;\omega=1$.
The comparison in Fig. 2 shows that the perturbation expansion yields valid
results only up to $\varepsilon\approx 0.5$. For small $\varepsilon$ one finds
in this case $V^{(1)}=0.00201\;\varepsilon^{2}+O(\varepsilon^{4})$. The
velocity function $V^{(1)}$ increases monotonically with $\varepsilon$.

##  VI Small amplitude

In this section we study the velocity function to second order on the
amplitude $\varepsilon=f/a$. First we formulate the asymptotic results for
small amplitude $f$ and large distance $d$. We find by substitution of the
above results into Eq. (4.3) for the velocity function to order $f^{2}/d^{2}$

$V_{as}=\frac{27\xi^{2}}{1+\xi}\;\frac{\sigma_{a}^{*}-\xi^{2}\sigma_{b}^{*}}{81(1+\xi)^{2}+4(\sigma_{a}^{*}+\xi^{3}\sigma_{b}^{*})^{2}R^{2}}\;\frac{f^{2}}{d^{2}}R,$
(70)

with dimensionless variables

$\xi=b/a,\qquad\sigma_{a}^{*}=\sigma_{a}+\frac{1}{2},\qquad\sigma_{b}^{*}=\sigma_{b}+\frac{1}{2}.$
(71)

The result Eq. (6.1) holds for all $R=a^{2}\omega\rho/\eta$, small $f$ and
large $d$.

For large viscosity $\eta$ or low frequency $\omega$ the behavior becomes

$V_{as}\approx\frac{1}{3}\;\frac{ab^{2}}{(a+b)^{3}}\;(a^{2}\rho_{a}^{*}-b^{2}\rho_{b}^{*})\;\frac{\omega
f^{2}}{\eta
d^{2}}\qquad(\mathrm{large}\;\eta,\mathrm{small}\;f\;\mathrm{and}\;\mathrm{large}\;d).$
(72)

For small viscosity or high frequency

$V_{as}\approx\frac{27}{4}\;\frac{ab^{2}}{a+b}\;\frac{a^{2}\rho_{a}^{*}-b^{2}\rho_{b}^{*}}{(a^{3}\rho_{a}^{*}+b^{3}\rho_{b}^{*})^{2}}\;\frac{\eta
f^{2}}{\omega
d^{2}}\qquad(\mathrm{small}\;\eta,\mathrm{small}\;f\;\mathrm{and}\;\mathrm{large}\;d).$
(73)

The expressions in Eqs. (6.1-4) become identical to those derived earlier
[5],[6] when $(\rho_{a}^{*},\rho_{b}^{*})$ is replaced by
$(\rho_{a},\rho_{b})$, i.e. if the single sphere added mass is ignored. The
expressions similar to Eqs. (6.3) and (6.4) derived by Derr et al.[10] are the
same apart from numerical factors. The potential flow appears only in the
single sphere added mass. The dipolar interaction does not contribute in Eq.
(6.1). We refer to the Oseen model with single sphere added mass included as
the Oseen* model. The velocity function $V^{(1)}_{O*}$ of this model is found
from Eq. (4.3) with $j_{1sm},\;m_{0i},\;m_{1s}$ put equal to zero. The
asymptotic behavior of the Oseen-Dipole model is the same as that of the
Oseen* model.

The analytic form of Eq. (6.1) is the same as that encountered earlier in Eq.
(4.5). Hence we conclude that for small amplitude $\varepsilon=f/a$ and long
distance $d$ the velocity function is maximal at frequency given by

$R_{asx}=\frac{9}{2}\;\frac{1+\xi}{\sigma_{a}^{*}+\xi^{3}\sigma_{b}^{*}}.$
(74)

The corresponding value of the velocity function is

$V_{asx}=\frac{3\xi^{2}}{4(1+\xi)^{2}}\;\frac{\sigma_{a}^{*}-\xi^{2}\sigma_{b}^{*}}{\sigma_{a}^{*}+\xi^{3}\sigma_{b}^{*}}\;\frac{f^{2}}{d^{2}}.$
(75)

With these values the asymptotic results of Eq. (4.9) correspond to Eqs. (6.5)
and (6.6).

One can choose $\xi$ such that $V_{asx}$ is maximized. This leads to a value
$\xi_{asx}$ in the interval $0<\xi<1$ which is independent of $f/d$. In
particular in case $\rho_{a}=\rho_{b}$ this becomes a root of the quartic
equation $\xi^{4}-3\xi^{3}+2\xi^{2}-4\xi+2=0$ with the numerical value
$\xi_{asx}=0.549$. In this case $\sigma_{a}^{*}=\sigma_{b}^{*}$ drops out in
Eq. (6.6), so that the value of the velocity function maximum $V_{asx}^{(1)}$
is independent of the mass density $\rho_{a}=\rho_{b}$. The value of the
maximum at $\xi_{asx}=0.549$ is
$V_{asxx}^{(1)}=0.05647\varepsilon^{2}/\delta^{2}$, also independent of the
mass density. The corresponding $R$-number is
$R_{asxx}^{(1)}=5.981/\sigma_{a}^{*}$. For neutrally buoyant spheres
$\sigma_{a}^{*}=3/2$ and $R_{asxx}^{(1)}=3.987$.

More generally we consider small amplitude $\varepsilon$, but any distance
$d$. It is clear from Eq. (4.3) that the numerator is of order
$\varepsilon^{2}$, whereas the denominator remains positive in the small
$\varepsilon$ limit. To second order in $\varepsilon$ Eq. (4.5) becomes

$V_{2}=v_{2}\varepsilon^{2}=\frac{A_{2}R}{B_{0}+C_{0}R^{2}},$ (76)

with $A_{2}$ given by

$A_{2}=\frac{1}{6}z_{1s,1}\big{[}2\pi
z_{0,0}(\sigma_{a}^{*}-2\xi^{3}\sigma_{b}^{*})-3z_{0,0}j_{1sm,1}-3m_{0,0}j_{1c,1}\big{]}\varepsilon^{2},$
(77)

with coefficients as given in the preceding section, and with $B_{0}$ and
$C_{0}$ given by

$B_{0}=z_{0,0}^{3},\qquad C_{0}=m_{0,0}^{2}z_{0,0}.$ (78)

We have dropped the superscript $(1)$, because the same result is obtained
when higher order harmonics are included. Eq. (6.7) is again of the form Eq.
(4.5), so that we find values $R_{2x}$ and $V_{2x}$ as in Eq. (4.6),

$R_{2x}=\frac{z_{0,0}}{m_{0,0}},\qquad
V_{2x}=\frac{A_{2}}{2z_{0,0}^{2}m_{0,0}}.$ (79)

From Eq. (5.3)

$z_{0,0}=6\pi(1+\xi)-18\pi\xi\frac{4\delta-3-3\xi}{4\delta^{2}-9\xi},$ (80)

and from Eqs. (3.16) and (5.8)

$\displaystyle m_{0,0}$ $\displaystyle=$
$\displaystyle\frac{4\pi}{3}(\sigma^{*}_{a}+\xi^{3}\sigma_{b}^{*}),\qquad(Oseen^{*})$
$\displaystyle m_{0,0}$ $\displaystyle=$
$\displaystyle\frac{4\pi}{3}(\sigma^{*}_{a}+\xi^{3}\sigma_{b}^{*})-2\pi\xi^{3}\frac{2\delta^{3}-1-\xi^{3}}{\delta^{6}-\xi^{3}},\qquad(Oseen-
Dip).$ (81)

The numerator $A_{2}$ in Eq. (6.10) reads for the Oseen*-model

$A_{2}=\varepsilon^{2}X[(2\delta-3)\sigma_{a}^{*}-(2\delta-3\xi)\xi^{2}\sigma_{b}^{*}],\qquad(Oseen^{*}),$
(82)

with prefactor

$X=576\pi^{3}\delta\xi^{2}\frac{(2\delta-3)(2\delta-3\xi)}{(4\delta^{2}-9\xi)^{3}}.$
(83)

Similarly for the Oseen-Dipole model

$A_{2}=\varepsilon^{2}\frac{X}{2(\delta^{6}-\xi^{3})}[r_{a}\sigma_{a}^{*}+r_{b}\sigma_{b}^{*}+r_{0}],\qquad(Oseen-
Dip),$ (84)

with coefficients

$\displaystyle r_{a}$ $\displaystyle=$ $\displaystyle
2(2\delta-3)(\delta^{6}-\xi^{3}),$ $\displaystyle r_{b}$ $\displaystyle=$
$\displaystyle
2\xi^{2}(\xi^{2}-\delta^{2})(2\delta-3\xi)(\delta^{4}+\delta^{2}\xi+\xi^{2}),$
$\displaystyle r_{0}$ $\displaystyle=$ $\displaystyle
3\xi^{2}(1-\xi)(2\delta^{4}-3\xi+2\delta\xi-3\xi^{2}+2\delta\xi^{2}-3\xi^{3}).$
(85)

We seek to optimize by finding the maximum value of $V_{2x}$ for fixed
$\delta$. This leads to the optimum $\xi_{x}(\delta)$, i.e. the optimal radius
of the second sphere for given length $d=\delta a$, and the corresponding
maximum $V_{2xx}(\delta)$. For both models the value $\xi_{x}(\delta)$ is
found as zero on the interval $0<\xi<1$ of a polynomial in $\xi$ of degree
eight. In Fig. 3 we plot the function $\xi_{x}(\delta)$ for both models for
the case of neutrally buoyant spheres.

In Fig. 4 we plot the corresponding functions
$V_{2xx}(\delta)/\varepsilon^{2}$. The plots nearly coincide with the function
$V_{asxx}/\varepsilon^{2}=0.05647/\delta^{2}$ found from the asymptotic
calculation below Eq. (6.6).

As we noted below Eq. (6.9), the first harmonics approximation yields the
exact small amplitude velocity function $V_{2}$ for the model considered. The
higher order correction terms in Sec. V can be used to calculate the first
harmonics approximation to the next order term $V_{4}=v_{4}\varepsilon^{4}$.
We would need to include the second harmonics in order to find the exact
$V_{4}$ for the model considered.

##  VII Large amplitude

The first harmonics approximation can be used also to perform model
calculations for large amplitude of stroke. In this section we examine some
situations of interest.

As a first example we show in Fig. 5 the velocity function $V$ of the Oseen-
Dipole model for the case
$\omega=1,\;a=1,\;b=0.5,\;d=3.5,\;\rho=\rho_{a}=\rho_{b}=1,\;f=2$ as a
function of the scaling number $R=a^{2}\omega\rho/\eta$, as calculated from
the exact expression Eq. (3.10) of Ref. 2, or alternatively, from the
expansion in harmonics with a sufficient number of higher order harmonics
(solid curve). The latter procedure is faster. We have included up to order
six.

We compare first with the asymptotic result Eq. (6.1) (short dashes), as well
as with this result with $\sigma_{a}^{*}$ replaced by $\sigma_{a}$
and$\sigma_{b}^{*}$ replaced by $\sigma_{b}$ (dot-dashed curve). The latter
result is identical with that derived by Hubert et al.[5]. For these values of
$f$ and $d$ the asymptotic result is rather different from the exact one.

We compare with the result in first harmonics approximation, found from Eq.
(4.3) (long dashes). On the scale of the figure the latter curve also
describes the result for the Oseen∗ model, suggesting that for this value of
$d$ the effect of the dipolar interaction is small. The long-dashed curve is
close to the solid one, showing that the first harmonics approximation is
nearly exact.

The velocity function $V^{(1)}$ as a function of
$R_{a}=a^{2}\omega\rho_{a}/\eta$ and $S_{a}=\rho/\rho_{a}$ is also found from
Eq. (4.3). The expression is of the form Eq. (4.5) with $R$ replaced by
$R_{a}$, with $A$ linear in $S_{a}$, $B$ independent of $S_{a}$, and $C$
quadratic in $S_{a}$.

Next we consider a swimmer corresponding to the experiments and lattice
Boltzmann simulations of Hubert et al. [5] with parameters
$a=8,\;b=5,\;d=28,\;\eta=1/6,\;\rho=1,\;\rho_{a}=\rho_{b}=8$, and determine
the velocity function maximum $V^{(1)}_{x}$ and corresponding number
$R^{(1)}_{x}$ as functions of the stroke amplitude $\varepsilon=f/a$ for the
Oseen-Dipole and the Oseen* model. From the asymptotic expressions Eqs. (6.5)
and (6.6) we find for both models $R^{(1)}_{asx}=0.691$ and
$V^{(1)}_{asx}=0.00444\;\varepsilon^{2}$. In Figs. 6 and 7 we plot the
functions $R^{(1)}_{x}(\varepsilon)$ and $V^{(1)}_{x}(\varepsilon)$ for the
two models. The models yield nearly identical results, but there is a
significant difference with the asymptotic values. In the limit
$\varepsilon\rightarrow 0$ one finds from Eq. (6.10) for the Oseen-Dipole
model $R_{2x}=0.5245$ and $V_{2x}=0.002194\;\varepsilon^{2}$, and for the
Oseen* model $R_{2x}=0.5236$ and $V_{2x}=0.002186\;\varepsilon^{2}$.

Returning to the general situation we seek to optimize the swimming speed by
appropriate choice of the ratio $\xi=b/a$ and the frequency $\omega$. For a
start we consider three cases with $d=3.5\;a$ fixed, but with different mass
densities of the spheres. In the first case $\rho_{a}=\rho_{b}=\rho=1$. In the
second case $\rho_{a}=\rho$, but $\rho_{b}=0$. A consideration of the contours
of $V^{(1)}$ in the $b\rho_{a}$-plane suggests that the latter choice
optimizes the velocity for fixed values of the other parameters. In the third
case both $\rho_{a}=0$ and $\rho_{b}=0$. We choose amplitude $f=d-a-b$, the
maximum amplitude for which there is no overlap.

We use the notation $\sigma=\rho_{b}/\rho_{a}$. In the first case with
$b=0.5\;a$ and $\sigma=1$ the velocity function takes the value
$V_{x}^{(1)}=0.00885$ at the maximum at $R_{x}^{(1)}=3.117$. The exact value
is $V_{x}=0.00859$. In the second case with $b=0.5\;a$ and $\sigma=0$ we find
the value $V_{x}^{(1)}=0.0133$ at the maximum at $R_{x}^{(1)}=3.371$. The
exact value is $V_{x}=0.0131$. In the third case with $b=a$ and $\sigma=0$ we
find the value $V_{x}^{(1)}=0.0141$ at the maximum at $R_{x}^{(1)}=3.270$. The
exact value is $V_{x}=0.0144$.

Next we consider a range of $d$ values. For each value of $\delta=d/a$ and
given mass densities we choose the ratio $\xi=b/a$ and the number $R$ such
that $V^{(1)}$ is maximized. For these values $\xi_{x}(\delta)$ we then
calculate $V^{(1)}$ in the form Eq. (4.5) and determine the values
$R_{xx}^{(1)}$ and $V_{xx}^{(1)}$ as functions of $\delta$. In the figures
below we compare these functions for the Oseen-Dipole model with those for the
Oseen∗ model.

We show first that for the Oseen∗ model with $\rho_{a}=\rho_{b}$ the maximum
value $V^{[1]}_{x}$ has the simple property that it is independent of the mass
density $\rho_{a}=\rho_{b}$ of the spheres, as we found below Eq. (6.6) for
the asymptotic situation. We note from Eq. (4.6) that
$V^{[1]}_{x}=A/(2\sqrt{BC})$ with $A,B,C$ read off from Eq. (4.3). In the
Oseen* model the moments $j_{1sm}$, $m_{0i}$, and $m_{1s}$ vanish. Hence in
this model $B=z_{0}^{3}$, $C=m_{0}^{2}z_{0}$ and

$V_{x}^{(1)}=\frac{A}{2m_{0}z_{0}^{2}},\qquad(Oseen^{*})$ (86)

With the symmetry $\rho_{a}=\rho_{b}$ one has

$\displaystyle A=\frac{2\pi}{3}z_{1s}\sigma^{*}_{a}[(1-\xi^{3})\varepsilon
z_{0}-2(1+\xi^{3})j_{1cf}],$ $\displaystyle
m_{0}=\frac{4\pi}{3}\sigma_{a}^{*}(1+\xi^{3}),\qquad(Oseen^{*},\rho_{a}=\rho_{b}).$
(87)

In Eq. (7.1) the factor $\sigma_{a}^{*}$ cancels and the other coefficients do
not depend on the mass densities, so that the velocity function maximum
$V_{x}^{(1)}$ is independent of the mass density $\rho_{a}=\rho_{b}$ of the
spheres. It follows that then also $\xi_{x}(\delta)$ is independent of
$\rho_{a}=\rho_{b}$. The above properties hold also for the Oseen model.

We choose amplitude $f=d-2a$. In Fig. 8 we show the function $\xi_{x}(\delta)$
for the case $\rho_{a}=\rho_{b}=\rho$ for both the Oseen-Dipole model and the
Oseen* model. In Fig. 9 we show the corresponding functions
$R_{xx}^{(1)}(\delta)$ and in Fig. 10 we show the functions
$V_{xx}^{(1)}(\delta)$. In particular, at $\delta=3.5$ we have
$\xi_{x}(3.5)=0.570$, $R_{xx}^{(1)}(3.5)=3.054$ , $V_{xx}^{(1)}(3.5)=0.00533$
for the Oseen-Dipole model and $\xi_{x}(3.5)=0.561$, $R_{xx}^{(1)}(3.5)=3.032$
, $V_{xx}^{(1)}(3.5)=0.00516$ for the Oseen* model.

According to the above argument for other values of $\rho_{a}=\rho_{b}$ the
curves for $\xi_{x}(\delta)$ and $V_{xx}^{(1)}(\delta)$ for the Oseen* model
are the same as shown in Figs. 9 and 10. The corresponding curves for the
Oseen-Dipole model are not much different. In Figs. 11 and 12 we show the
curves for $\xi_{x}(\delta)$ and $R_{xx}^{(1)}(\delta)$ at
$\rho_{a}=\rho_{b}=0.5\rho$ for both models. The curve for
$V_{xx}^{(1)}(\delta)$ for the Oseen* model is the same as the one in Fig. 10,
and the one for the Oseen-Dipole model is quite similar. In particular at
$\delta=3.5$ we have
$\xi_{x}(3.5)=0.574,\;R_{xx}^{(1)}(3.5)=4.555,\;V_{xx}^{(1)}(3.5)=0.00542$ for
the Oseen-Dipole model and
$\xi_{x}(3.5)=0.561,\;R_{xx}^{(1)}(3.5)=4.529,\;V_{xx}^{(1)}(3.5)=0.00516$ for
the Oseen* model.

Consider again the Oseen-Dipole swimmer corresponding to Fig. 5,
$a=1,\;b=0.5,\;d=3.5$, but with amplitude of stroke $f=1.5$. We choose units
such that the fluid has viscosity $\eta=1$ and mass density $\rho=1$. If the
mass density of both spheres is $\rho_{a}=\rho_{b}=1$, then the two-sphere in
first harmonics approximation is found from Eq. (4.6) to have velocity
function maximum $V^{(1)}_{x}=0.00513$, achieved at $R$-number $R_{x}=3.129$.
As we showed above, if we increase the radius of the second sphere to
$\xi_{x1}=0.570$ then the velocity function maximum increases to
$V^{(1)}_{xx}=0.00533$, achieved at $R$-number $R_{xx}=3.054$. If instead we
first decrease the mass density of the spheres to $\rho_{a}=\rho_{b}=0.5$, and
increase the radius of the second sphere to $\xi_{x2}=0.574$, then we achieve
velocity function maximum $V^{(1)}_{xx}=0.00542$, at $R$-number
$R_{xx}=4.555$.

##  VIII Power and efficiency

In the preceding sections we concentrated on the calculation of the mean
swimming velocity. Other properties of both the Oseen* model and the Oseen-
Dipole model can also be investigated in full detail. In this section we
consider the power dissipated per period and the efficiency of swimming. The
latter is defined as the ratio of speed and power.

In both models the relative distance between both centers is prescribed as
$x_{2}(t)-x_{1}(t)=d+f\sin(\omega t)$. The center position
$C(t)=\frac{1}{2}(x_{1}(t)+x_{2}(t))$ is found by integration of the swimming
velocity $U_{sw}(t)=dC(t)/dt$. Hence the two sphere centers move as
$x_{1}(t)=C(t)-(d+f\sin(\omega t))/2$ and $x_{2}(t)=C(t)+(d+f\sin(\omega
t))/2$ The two sphere velocities are $U_{1}(t)=U_{sw}(t)-\frac{1}{2}\omega
f\cos(\omega t)$ and $U_{2}(t)=U_{sw}(t)+\frac{1}{2}\omega f\cos(\omega t)$.

The time-dependent rate of dissipation is given by

$\mathcal{D}={\bf\sf U}\cdot\mbox{\boldmath$\zeta$}\cdot{\bf\sf U},$ (88)

where ${\bf\sf U}=(U_{1},U_{2})$. The power used for a stroke in periodic
swimming is the average dissipation in a period,

$P=\overline{\mathcal{D}}.$ (89)

The efficiency is defined as [13]

$L=8\pi\eta\omega a^{2}\frac{\overline{U_{sw}}}{\overline{\mathcal{D}}}.$ (90)

The dimensionless mean rate of dissipation $\hat{\mathcal{D}}$ is defined by

$\overline{\mathcal{D}}=8\pi\eta\omega^{2}a^{3}\hat{\mathcal{D}},$ (91)

so that $L=V/\hat{\mathcal{D}}$.

The mean rate of dissipation $\hat{\mathcal{D}}$ can be calculated in first
harmonics approximation in analytic form, but the expression for
$\hat{\mathcal{D}}^{(1)}$ becomes too complicated to be practically useful.
The required values can be calculated from the numerical calculation as before
[12].

Depending on the starting point and the size of the shift the change of
$V^{(1)}$ upon doubling the viscosity or the mass density can be positive,
negative, or zero. In Table I we list values of $V^{(1)}$ and $L^{(1)}$ for
selected values of $R_{a}=a^{2}\omega\rho_{a}/\eta$ and $S_{a}=\rho/\rho_{a}$.
These can be compared with the corresponding values of $V$ and $L$ in Table I
of Ref. 11, calculated by inclusion of higher order harmonics. In the latter
Table some copying errors need to be corrected. The last but third column of
that Table should have heading $10^{5}L$ and entries $631,1000,873$.
Comparison of the two Tables shows that for this large amplitude the
approximate value $V^{(1)}$ differs from the exact $V$ by about ten percent.
Table I shows that the efficiency varies significantly over the
$R_{a}S_{a}$-plane. We note that $R=R_{a}S_{a}$.

Table 1:  $R_{a},S_{a}$ | $10^{6}V^{(1)}$ | $10^{6}L^{(1)}$ | $R_{a},S_{a}$ | $10^{5}V^{(1)}$ | $10^{5}L^{(1)}$ | $R_{a},S_{a}$ | $10^{5}V^{(1)}$ | $10^{5}L^{(1)}$ | $R_{a},S_{a}$ | $10^{6}V^{(1)}$ |  $10^{5}L^{(1)}$
---|---|---|---|---|---|---|---|---|---|---|---
$0.1,1$ | $569$ | $752$ | $1,1$ | $516$ | $679$ | $10,1$ | $504$ | $638$ | $100,1$ | $553$ |  $70$
$0.05,1$ | $284$ | $376$ | $0.5,1$ | $277$ | $367$ | $5,1$ | $797$ | $1017$ | $50,1$ | $1102$ |  $139$
$0.05,2$ | $386$ | $508$ | $0.5,2$ | $367$ | $486$ | $5,2$ | $702$ | $886$ | $50,2$ | $856$ |  $107$

-

###

List of values of the velocity function $V^{(1)}(R_{a},S_{a})$ and the
efficiency function $L^{(1)}(R_{a},S_{a})$ for the Oseen-Dipole model of the
two-sphere with $\xi=0.5,\;\delta=3.5,\;\varepsilon=2,\;\rho_{a}=\rho_{b}$.

In Ref. 11 we found for a particular swimmer that the surface $L(R_{a},S_{a})$
looks quite similar to the surface $V(R_{a},S_{a})$. This implies that
$\hat{\mathcal{D}}$ is nearly constant in the $R_{a}S_{a}$-plane. In Fig. 13
we show the surface $\hat{\mathcal{D}}^{(1)}(R_{a},S_{a})$ for the case
$\xi=0.5,\;\delta=3.5,\;\varepsilon=2,\;\rho_{a}=\rho_{b}=\rho$.

##  IX Discussion

The study shows that the first harmonics approximation provides a useful tool
for the analysis of the two-sphere Oseen-Dipole model and of the Oseen* model.
The approximation yields an analytic expression for the mean swimming velocity
which can be studied in its dependence on parameters. The swimming speed can
be compared with the exact one as found from another expression derived
earlier [3], or as found numerically from an expansion including higher order
harmonics. The comparison in special cases with large amplitude of stroke
shows that the first harmonics approximation leads to results in good
agreement with the exact ones. For small amplitude the first harmonics
approximation becomes exact. In Sec. VI we derived the explicit expression for
the small amplitude mean swimming velocity for both models.

The first harmonics expression for the mean swimming velocity Eq. (4.3) shows
a simple dependence on the scaling numbers $R=a^{2}\omega\rho/\eta$ and
$S_{a}=\rho/\rho_{a}$, consisting of a Padé type ratio of two simple
polynomials in terms of the two numbers. The dependence on $R$ shows that the
mean swimming velocity in any chosen fluid varies slowly over a wide range of
frequency of the stroke. The expression allows us to optimize the swimming
velocity by adjusting the frequency and the ratio of radii.

The Oseen-Dipole model is more sophisticated than the Oseen∗ model since it
includes the dipolar interactions of potential flow. However, it fails to
describe the interesting reversal of swimming direction at high frequency seen
by Dombrowski and Klotsa [9] in a numerical analysis based on the Navier-
Stokes equations. We view the present study as a necessary prologue to further
investigation. We expect that the model can be adapted with more realistic
hydrodynamic interactions to describe the reversal phenomenon.

The author has no conflicts of interest to disclose.

## References

  * [1]
  * [2] D. Klotsa, K. A. Baldwin, R. J. A. Hill, R. M. Bowley, and M. R. Swift, Phys. Rev. Lett. 115, 248102 (2015).
  * [3] B. U. Felderhof, Phys. Rev. E 94, 063114 (2016).
  * [4] E. M. Purcell, Am. J. Phys. 45, 3 (1977).
  * [5] M. Hubert, O. Trosman, Y. Collard, A. Sukhov, J. Harting, N. Vandewalle, and A.-S. Smith, Phys. Rev. Lett. 126, 224501 (2021).
  * [6] B. U. Felderhof, arXiv:2107.10629[physics.flu-dyn.].
  * [7] B. U. Felderhof, Phys. Rev. E 92, 053011 (2015).
  * [8] H. Lamb, Hydrodynamics (Dover Publications, New York, 1945).
  * [9] T. Dombrowski and D. Klotsa, Phys. Rev. Fluids 5, 063103 (2020).
  * [10] N. J. Derr, T. Dombrowski, C. H. Rycroft, and D. Klotsa, arXiv:2202.03669 [physics.flu-dyn].
  * [11] N. Riley, Quart. J. Mech. and Applied Math. 19, 461 (1966).
  * [12] B. U. Felderhof, Phys. Fluids 34, 011903 (2022).
  * [13] B. U. Felderhof and R. B. Jones, Eur. J. Mech./B Fluids 85, 58 (2021).

Figure 1: Plot of $\chi^{\prime\prime}(\omega)$ given by Eq. (4.12) (dashed
curve) and of $\chi^{\prime\prime}_{2}(\omega)$ given by Eq. (4.13) (solid
curve) for the numerical example
$\chi_{0}=1,\;\gamma=1,\;\chi_{02}=-1,\;\gamma_{2}=10$.

\put(-9.1,3.1){}

Figure 2: Plot of the velocity function $V^{(1)}$ as a function of amplitude
$\varepsilon=f/a$ for the case
$\eta=0.5,\;\rho=1,\;a=1,\;b=0.5,\;d=3.5,\;\rho_{a}=\rho_{b}=1,\;\omega=1$
(solid curve) and of the approximation given by the expansion of coefficients
in powers of $\varepsilon$ as derived in Sec. V (dashed curve). On the scale
of the figure the dashed curve cannot be distinguished from the expansion up
to the quartic term, as given by Eq. (5.1).

\put(-9.1,3.1){}

Figure 3: Plot of the function $\xi_{x}(\delta)$ in the small amplitude limit
for a two-sphere with mass densities $\rho_{a}=\rho_{b}=\rho$ for the Oseen-
Dipole model (solid curve) and the Oseen*-model (long dashes). The horizontal
line (short dashes) is the asymptotic limit $\xi_{asx}=0.549$.

\put(-9.1,3.1){}

Figure 4: Plot of the optimal velocity function $V_{2xx}(\delta)$ for a two-
sphere with mass densities $\rho_{a}=\rho_{b}=\rho$ for the Oseen-Dipole model
(solid curve) and the Oseen*-model (long dashes). The plots nearly coincide
with the function $V_{asxx}/\varepsilon^{2}=0.05647/\delta^{2}$ found from the
asymptotic calculation below Eq. (6.6).

\put(-9.1,3.1){}

Figure 5: Plot of the velocity function $V$ for parameters
$\xi=b/a=0.5,\;\delta=d/a=3.5,\;\varepsilon=f/a=2,\;\rho_{a}=\rho_{b}=\rho=1$
as a function of $R=a^{2}\omega\rho/\eta$ (solid curve). We compare with the
velocity function $V^{(1)}$ calculated in first harmonics approximation (short
dashes) and with the asymptotic result Eq. (6.1) (long dashes), and the
analogous result for the Oseen model [5] (dot-dashed curve).

\put(-9.1,3.1){}

Figure 6: Plot of the function $R^{(1)}_{x}(\varepsilon)$ for a two-sphere
with parameters $a=8,\;b=5,\;d=28,\;\eta=1/6,\;\rho=1,\;\rho_{a}=\rho_{b}=8$,
as a function of $\varepsilon=f/a$ for the Oseen-Dipole model (solid curve)
and the Oseen∗ model (dashed curve).

\put(-9.1,3.1){}

Figure 7: As in Fig. 6 for the function $V^{(1)}_{x}(\varepsilon)$.

\put(-9.1,3.1){}

Figure 8: Plot of the function $\xi_{x}(\delta)$ representing the ratio $b/a$
with optimal $V^{(1)}_{x}$ for a two-sphere with mean distance $d$ between
centers and stroke of amplitude $\varepsilon=\delta-1-\xi_{x}$, as a function
of $\delta=d/a$ for the Oseen-Dipole model (solid curve) and the Oseen∗ model
(dashed curve) with mass densities $\rho_{a}=\rho_{b}=\rho$.

\put(-9.1,3.1){}

Figure 9: Plot of the function $R^{(1)}_{xx}(\delta)$ corresponding to the
optimum specified in caption 8.

\put(-9.1,3.1){}

Figure 10: Plot of the function $V^{(1)}_{xx}(\delta)$ corresponding to the
optimum specified in caption 8.

\put(-9.1,3.1){}

Figure 11: As in Fig. 8 for a two-sphere with mass densities
$\rho_{a}=\rho_{b}=0.5\rho$.

\put(-9.1,3.1){}

Figure 12: As in Fig. 9 for a two-sphere with mass densities
$\rho_{a}=\rho_{b}=0.5\rho$.

Figure 13: Plot of the surface $\hat{\mathcal{D}}^{(1)}(R_{a},S_{a})$, as
defined by Eqs. (8.1-4), for the case
$\xi=0.5,\;\delta=3.5,\;\varepsilon=2,\;\rho_{a}=\rho_{b}=\rho$.