# Self-propulsion of an elliptic phoretic disk emitting solute uniformly

Guangpu Zhu  Department of Mechanical Engineering, National University of
Singapore, 117575, Singapore  Lailai Zhu<EMAIL_ADDRESS>Department of
Mechanical Engineering, National University of Singapore, 117575, Singapore

###### Abstract

Self-propulsion of chemically active droplet and phoretic disk has been widely
studied; however, most research overlooks the influence of disk shape on
swimming dynamics. Inspired by the experimentally observed prolate composite
droplets and elliptic camphor disks, we employ simulations to investigate the
phoretic dynamics of an elliptic disk that uniformly emits solutes in the
creeping flow regime. By varying the disk’s eccentricity $e$ and the P’eclet
number $\mathrm{Pe}$, we distinguish five disk behaviors: stationary, steady,
orbiting, periodic, and chaotic. We perform a global linear stability analysis
(LSA) to predict the onset of instability and the most unstable eigenmode when
a stationary disk spontaneously transitions to steady self-propulsion. In
addition to the LSA, we use an alternative approach to determine the
perturbation growth rate, offering valuable insights into the competing roles
of advection and diffusion. The steady motion features a transition from a
puller-type to a neutral-type swimmer as $\mathrm{Pe}$ increases, which occurs
as a bimodal concentration profile at the disk surface shifts to a polarized
solute distribution, driven by convective solute transport. An elliptic disk
achieves an orbiting motion through a chiral symmetry-breaking instability,
wherein it repeatedly follows a circular path while simultaneously rotating.
The swinging periodic motion, emerging from a steady motion via a
supercritical Hopf bifurcation, is characterized by a wave-like trajectory. We
uncover a transition from normal diffusion to superdiffusion as eccentricity
$e$ increases, corresponding to a random-walking circular disk and a
ballistically swimming elliptic counterpart, respectively.

##  I Introduction

Synthetic micro-swimmers have attracted much attention owing to their
promising potential in biomedical and bioengineering applications [70], e.g.,
detection and collection of metal ions [2], targeted controlled drug delivery
[19, 76], and cancer-cells microsurgery [80]. Drawing inspiration from the
propulsion strategies of microorganisms in nature, various bio-mimetic
swimmers that propel in viscous fluids have been developed [21, 20, 38, 71].
Unlike their biological counterparts, these synthetic micro-swimmers are
commonly powered by external forces or torques coming from the electric,
optic, acoustic, or magnetic fields [10, 7, 66, 38, 23], for example, a sperm-
mimicking micro-swimmer with a flexible filament actuated magnetically [15].
Despite the rapid development of externally actuated micro-swimmers, some
practical difficulties, such as miniaturization and manufacturing of moving
parts for certain swimmers, have limited their applications in realistic
scenarios [18, 33, 43, 64].

Unlike externally actuated swimmers, chemically active swimmers covert
chemical energy stored internally or extracted from their surroundings into
motion [53]. They can be broadly classified by whether their surface
properties, e.g., surface activity and mobility, are anisotropic or isotropic.
A classical anisotropic swimmer is the Janus colloid, e.g., the autophoretic
Au-Pt Janus colloid [60]. The chemically patterned asymmetric colloid
typically features two distinct compartments, each composed of a different
material or bearing diverse functional groups [40], which enables asymmetric
chemical reactions at the surface. The inherent asymmetry allows it to self-
generate a concentration gradient, which drives a slip flow inducing net
phoretic propulsion, as revealed by experimental [61, 16, 8], theoretical [22,
6, 57, 11] and numerical [65, 69, 37] studies. The Janus swimmer is typically
micro-scale or even smaller, whose phoretic motion is Brownian [51]. Its self-
propulsion requires a built-in asymmetry in the surface properties. This
requirement presents a challenge to the controlled and reproducible
manufacturing of Janus colloids, hence hindering their high-throughput
production [52].

Chemically isotropic swimmers are much easier to manufacture compared to their
anisotropic counterparts. A simple and typical representative of such swimmers
is a chemically active droplet, e.g., a water droplet slowly dissolving in a
surfactant-saturated oil phase [4], which have been extensively researched
since their first experimental realization [30]. These active droplets are
generally larger than the Janus microswimmers and have typical radii ranging
from 10 to 100 $\mu\mathrm{m}$. Active droplets mainly consist of reacting
droplets [78, 35, 74] and solubilizing droplets [62, 68, 26]. The former
involve chemical reactions producing or changing surfactant molecules at their
surface, while the latter feature a micellar dissolution into the surfactant-
saturated ambient phase. In both instances, spatial modulation of the surface
tension at the droplet interface may potentially induce Marangoni flows. Such
droplets do not rely on built-in asymmetry like the Janus colloids, but
instead attain self-propulsion via an instability spontaneously breaking the
spatial symmetry. This instability arises from the nonlinear convective
transport of solute species by fluid flow, resulting from Marangoni and/or
phoretic effects produced by local chemical gradients surrounding the droplet
[54, 64]. Besides the active droplet, another class of chemically isotropic
swimmers is camphor disks that surf at a liquid-air interface [79, 56, 73]. In
this scenario, camphor molecules dissolved from the disk diffuse into the
interface and further the subsurface liquid, and the Marangoni flows resulting
from the solute gradient drive the disk to propel [48, 5]. Notably, the
Marangoni flows generated by active droplets are solely at their surface,
while those triggered by camphor disks are along the air-liquid interface and
depend significantly on the depth of the subsurface liquid [48, 51].

Active droplets exhibit complex and tunable motion as a result of the
nonlinear physico-chemical hydrodynamics [26, 42], characterized by the Péclet
number $\mathrm{Pe}$ as the ratio of flow advection to solute diffusion. At
low $\mathrm{Pe}$, an isolated droplet remains stationary. Morozov & Michelin
[55] identified the critical Péclet number for an undeformable droplet through
a stability analysis, beyond which an unstable dipolar mode of hydrodynamics
emerges, driving the spontaneous propulsion of the droplet. This critical
$\mathrm{Pe}$ exactly matches that for a spherical phoretic particle [52]. In
fact, the spherical particle can be considered as a droplet with infinite
viscosity. Besides, the critical $\mathrm{Pe}$ determined for a two-
dimensional (2D) undeformable droplet [42] is also consistent with that for a
phoretic disk [27]. Near the critical $\mathrm{Pe}$, the dipolar mode is the
only unstable one [67, 63]. However, higher-order modes, e.g., the quadrupolar
mode, becomes successively unstable with increasing $\mathrm{Pe}$, leading to
the possible coexistence of multiple unstable modes with different polar
symmetries [54, 26]. Accordingly, the active droplet sequentially exhibits
quasi-ballistic, unsteady curvilinear, and even chaotic motions [39, 75, 26,
42] as $\mathrm{Pe}$ grows. Analogous behaviors of 2D [27, 28] and three-
dimensional (3D) isotropic phoretic particles have also been observed [9, 28,
34]. Besides an isolated unbounded droplet/particle, the effect of nearby
boundaries [77, 31, 4, 46, 12, 14, 64], that of an ambient flow [83, 13, 17],
and interaction among multiple droplets/particles [32, 45, 49, 58, 81, 24, 82]
have been investigated.

Most of the active droplets observed in experiments were weakly deformed and
remained spherical, one exception is the very recent work [25] reporting the
self-propulsion of a prolate composite droplet along its minor axis. That oil
droplet trapping two aqueous daughter droplets at the opposing poles of its
major axis, suspended in an aqueous surfactant solution. The daughter droplets
are submicellar, whereas the external aqueous phase is supramicellar. The
micellar dissolution at the external oil-water interface induces a self-
sustaining surface tension gradient, driving the droplet motion. Compared to
droplets, solid self-propelling swimmers relying on a similar symmetry-
breaking mechanism exhibit greater flexibility in shape. Kitahata et al. [36]
and [29] experimentally and theoretically investigated the spontaneous motion
of an elliptic camphor disk at the air-liquid interface, and found that the
disk swam along its minor axis resembling the swimming prolate droplet [25].
Motivated by these non-spherical phoretic swimmers with uniform chemical
reactions at their surface, especially Kitahata et al. [36] and Hokmabad et
al. [25], here, we theoretically and numerically explore, in the creeping flow
regime, the instability-driven spontaneous propulsion of an elliptic phoretic
disk that releases chemical species uniformly. We perform a global linear
stability analysis to investigate the onset of instability, and direct
numerical simulations to explore the swimming behavior of the phoretic disk.

This paper is organized as follows. We describe the problem setup,
assumptions, and governing equations in Section II. The implementation for the
global LSA is introduced in Section III, followed by Section IV demonstrating
numerical and theoretical results. Finally, we conclude our observations and
provide some discussion in Section V.

##  II Problem setup, governing equations and methodology

###  II.1 Problem setup and governing equations

We consider a chemically active elliptic disk emitting or absorbing solute
molecules uniformly in an incompressible Newtonian fluid of dynamic viscosity
$\bar{\eta}$ (see figure 1). From here on, $\bar{}$ indicates dimensional
variables unless otherwise mentioned. The semi-major and semi-minor axes of
the disk are $\bar{a}$ and $\bar{b}\leq\bar{a}$, respectively, and
$\bar{c}_{f}=\sqrt{\bar{a}^{2}-\bar{b}^{2}}$ denotes the half of its focal
length. Hence, the disk shape can be characterized by the eccentricity
$e=\bar{c}_{f}/\bar{a}$, which recovers to $0$ or approaches $1$ as the disk
becomes circular or needle-like. We choose the major axis of the disk to
denote its orientation
$\mathbf{e}_{s}=\sin\theta\mathbf{e}_{x}+\cos{\theta}\mathbf{e}_{y}$, which is
characterized by its angular deviation from $\mathbf{e}_{y}$.

Figure 1: Self-propulsion of an elliptic disk uniformly releasing chemical
solutes in a Newtonian solvent. $\mathbf{e}_{xy}$ denotes the laboratory
frame. The disk moves along an undulatory path with a translational velocity
of $\mathbf{U}$ and the color of the path is coded by the time $t$. The disk
orientation $\mathbf{e}_{s}$ coinciding with its major axis deviates from
$\mathbf{e}_{y}$ and $\mathbf{U}$ by angles $\theta$ and $\alpha$,
respectively. The inset shows the induced slip velocity $\mathbf{u}_{s}$ by
local solute gradients within a thin boundary layer of thickness $\lambda\ll
b$. $\mathbf{n}$ denotes the unit normal vector pointing away from the disk
surface. All the variables here are dimensionless.

We now describe the phoretic dynamics of the elliptic disk and the associated
governing equations. The disk’s surface $\Gamma_{d}$ uniformly emits or
absorbs solutes with a constant rate $\bar{A}$ (activity), hence,

$\left.\bar{D}\mathbf{n}\cdot\nabla\bar{c}\right|_{\Gamma_{d}}=-\bar{A}.$ (1)

Here, $\bar{c}$ is the solute concentration, $\bar{D}$ denotes the molecular
diffusivity of the solute, and $\mathbf{n}$ is the unit outward normal at the
surface. The solute interacts with the disk surface through a short-range
potential and here we focus on the classical thin-interaction-layer limit
$\bar{\lambda}\ll\bar{b}$, with $\lambda$ the thickness of the interaction
layer. Within this layer, the slip velocity along the disk surface induced by
the local tangential solute gradients reads

$\left.\bar{\mathbf{u}}_{s}\right|_{\Gamma_{d}}=\bar{M}\nabla_{s}\bar{c},$ (2)

where $\bar{\mathbf{u}}_{s}$ indicates the slip velocity,
$\nabla_{s}=\left(\mathbf{I}-\mathbf{n}\mathbf{n}\right)\cdot\nabla$ is the
surface gradient operator, and $\bar{M}$ denotes the constant phoretic
mobility. The signs of $\bar{A}$ and $\bar{M}$ depend on the physico-chemical
properties of the disk. Positive or negative $\bar{A}$ corresponds to the
solute emission or absorption at the disk surface, respectively. Note that the
spontaneous symmetry-breaking propulsion occurs only when $\bar{A}\bar{M}>0$
[52, 27], and the disk remains stable otherwise. Without any loss of
generality, this work will only consider positive $\bar{A}$ and $\bar{M}$.
Also, we assume that the inertia of both fluid and disk is negligible compared
to the viscous force because the Reynolds number $\mathrm{Re}$ is typically
small in the experiments [62, 47, 4, 26, 51, 24]. Hence, the fluid flow
surrounding the disk can be described approximately by the Stokes equation.

In the following, we introduce the dimensionless governing equations. We
choose $\bar{A}\bar{M}/\bar{D}$, $\bar{b}$, and $\bar{b}\bar{A}/\bar{D}$,
respectively, as the characteristic velocity $\bar{\mathcal{V}}$, length, and
concentration for the nondimensionalization. All variables below are
dimensionless unless otherwise specified. The dimensionless equations for the
velocity $\mathbf{u}$, pressure $p$, and concentration $c$ are

$\nabla\cdot\mathrm{\bm{\sigma}}=\mathbf{0},\;\;\;\nabla\cdot\mathbf{u}=0,$
(3) $\frac{\partial c}{\partial t}+\mathbf{u}\cdot\nabla
c=\frac{1}{\mathrm{Pe}}\Delta c,$ (4)

where the Péclet number $\mathrm{Pe}=\bar{A}\bar{M}\bar{b}/\bar{D}^{2}$
measures the ratio of flow advection to solute diffusion and
$\mathrm{\bm{\sigma}}=-p\mathbf{I}+\nabla\mathbf{u}+\left(\nabla\mathbf{u}\right)^{T}$
is the hydrodynamic stress tensor. At the disk surface, the constant flux
boundary condition (1) for the concentration reads

$\left.\mathbf{n}\cdot\nabla c\right|_{\Gamma_{d}}=-1.$ (5)

We consider a circular fluid domain with a finite radius $R$ and prescribe on
its exterior $\Gamma_{o}$

$\left.c\right|_{\Gamma_{o}}=0.$ (6)

It should be noted that equation (4) reduces to the Laplace equation when
$\mathrm{Pe}\rightarrow 0$, and the imposition of this Dirichlet boundary
condition at the outer boundary at a finite distance ensures the well-
posedness of our 2D model [27, 51]. Additionally, the domain radius $R$ is
known to affect the disk motion. Following Hu et al. [27] and Li [42], we set
$R=200$.

We perform numerical simulations in the frame co-moving with the disk center.
Hence, the boundary conditions for the velocity at the disk surface
$\Gamma_{d}$ and the outer boundary $\Gamma_{o}$ are

$\displaystyle\mathbf{u}|_{\Gamma_{d}}$
$\displaystyle=\mathbf{u}_{s}+\mathrm{\bm{\Omega}}\times\left(\mathbf{x}_{s}-\mathbf{x}_{c}\right),$
(7a) $\displaystyle\mathbf{u}|_{\Gamma_{o}}$ $\displaystyle=-\mathbf{U},$ (7b)

where $\mathbf{u}_{s}=\nabla_{s}c$ is the slip velocity at the disk surface,
$\mathbf{x}_{s}$ and $\mathbf{x}_{c}$ denote the coordinates of a general
point at the disk surface and the disk center, respectively. Here,
$\mathbf{U}$ and $\mathrm{\bm{\Omega}}$ denote the translational and
rotational velocities of the disk, respectively, which are determined by the
force-free and torque-free conditions [41]

$\displaystyle\mathbf{F}=\int\mathbf{n}\cdot\mathrm{\bm{\sigma}}\;\mathrm{d}\Gamma_{d}$
$\displaystyle=\mathbf{0},$ (8a)
$\displaystyle\mathbf{T}=\int\left(\mathbf{x}_{s}-\mathbf{x}_{c}\right)\times\left(\mathbf{n}\cdot\mathrm{\bm{\sigma}}\right)\mathrm{d}\Gamma_{d}$
$\displaystyle=\mathbf{0}.$ (8b)

Namely, the total force and torque exerted on the disk swimmer are zero.
Equations (3)-(8) form the complete set of governing equations for an elliptic
phoretic disk in the frame of co-moving with the disk.

###  II.2 Numerical method

We numerically solve the governing equations, using a finite-element-method
solver implemented in the commercial package COMSOL Multiphysics (I-Math,
Singapore). We adopt the moving mesh technique to tackle the deformation of
the fluid domain caused by disk rotation. Taylor-Hood and quadratic Lagrange
elements are employed to discretize the flow field,
$\left(\mathbf{u},p\right)$, and the concentration $c$, respectively. The
computational domain is discretized by approximately $85000-127000$ triangular
elements, and the mesh is locally refined near the disk. Our COMSOL
implementations are extensively validated against several published datasets,
as shown in Appendix A.

##  III Global linear stability analysis

Prior studies reveal that a chemically isotropic disk/droplet recovers a
stationary state, termed base state, at a sufficiently low $\mathrm{Pe}$. As
$\mathrm{Pe}$ grows beyond a critical value, an instability arises, leading
the disk/droplet to swim autonomously [27, 55, 42]. Our numerical results
indicate that the elliptic disk here exhibits analogous behavior, hence we
conduct a global LSA to examine the onset of global instability at a critical
Péclet number $\mathrm{Pe}_{c}^{(1)}$. The term ‘global instability’ here
refers to the instability of a base state within a domain with either no
homogeneous spatial direction or only one (in the case of a 3D domain). For
the 2D domain under consideration, the global instability recovers to the
local counterpart as the base state is homogeneous in one out of two spatial
directions [59].

We first decompose the space- ($\mathbf{x}$) and time-dependent state
variables $\left(c,\mathbf{u},p\right)$ into the sum of a base state and a
perturbation state as

$\displaystyle c\left(\mathbf{x},t\right)$
$\displaystyle=c_{b}\left(\mathbf{x}\right)+c^{\prime}\left(\mathbf{x},t\right),$
(9a) $\displaystyle\mathbf{u}\left(\mathbf{x},t\right)$
$\displaystyle=\mathbf{u}_{b}\left(\mathbf{x}\right)+\mathbf{u}^{\prime}\left(\mathbf{x},t\right),$
(9b) $\displaystyle p\left(\mathbf{x},t\right)$
$\displaystyle=p_{b}\left(\mathbf{x}\right)+p^{\prime}\left(\mathbf{x},t\right),$
(9c)

where the subscript $b$ denotes base-state fields and the primed variables are
infinitesimal perturbations. By substituting (9) into (3) and (4), and
retaining linear terms, we obtain

$\nabla\cdot\mathrm{\bm{\sigma}}^{\prime}=\mathbf{0},\;\;\nabla\cdot\mathbf{u}^{\prime}=0,$
(10) $\frac{\partial c^{\prime}}{\partial t}+\mathbf{u}_{b}\cdot\nabla
c^{\prime}+\mathbf{u}^{\prime}\cdot\nabla c_{b}=\frac{1}{\mathrm{Pe}}\Delta
c^{\prime}.$ (11)

The perturbations are assumed to vary exponentially in time with a complex
growth rate $\sigma=\sigma_{r}+i\sigma_{i}$, i.e.,

$\displaystyle c^{\prime}\left(\mathbf{x},t\right)$
$\displaystyle=\hat{c}\left(\mathbf{x}\right)e^{\sigma t},$ (12a)
$\displaystyle\mathbf{u}^{\prime}\left(\mathbf{x},t\right)$
$\displaystyle=\hat{\mathbf{u}}\left(\mathbf{x}\right)e^{\sigma t},$ (12b)
$\displaystyle p^{\prime}\left(\mathbf{x},t\right)$
$\displaystyle=\hat{p}\left(\mathbf{x}\right)e^{\sigma t}.$ (12c)

Consequently, equations (10) and (11) can be reformulated to

$\nabla\cdot\hat{\mathrm{\bm{\sigma}}}=\mathbf{0},\;\;\nabla\cdot\hat{\mathbf{u}}=0,$
(13)
$\sigma\hat{c}+\mathbf{u}_{b}\cdot\nabla\hat{c}+\hat{\mathbf{u}}\cdot\nabla
c_{b}=\frac{1}{\mathrm{Pe}}\Delta\hat{c}.$ (14)

By expanding the translational and rotational velocities similarly, we arrive
$\mathbf{U}=\mathbf{U}_{b}+\hat{\mathbf{U}}e^{\sigma t}$ and
$\mathrm{\bm{\Omega}}=\mathrm{\bm{\Omega}}_{b}+\hat{\mathrm{\bm{\Omega}}}e^{\sigma
t}$. Substituting them with (9) into (5) to (7) enables us to derive the
boundary conditions for $\hat{c}$ and $\hat{\mathbf{u}}$ at the disk surface
and outer boundary

$\left.\mathbf{n}\cdot\nabla\hat{c}\right|_{\Gamma_{p}}=0,\;\;\;\left.\hat{c}\right|_{\Gamma_{o}}=0,$
(15)
$\left.\hat{\mathbf{u}}\right|_{\Gamma_{p}}=\nabla_{s}\hat{c}+\hat{\mathrm{\bm{\Omega}}}\times(\mathbf{x}_{s}-\mathbf{x}_{c}),\;\;\;\left.\hat{\mathbf{u}}\right|_{\Gamma_{o}}=-\hat{\mathbf{U}}.$
(16)

Note that $\mathbf{U}_{b}$ and $\mathrm{\bm{\Omega}}_{b}$ disappear in (16)
corresponding to a stationary disk of the base state. The force-free and
torque-free conditions still hold as

$\displaystyle\hat{\mathbf{F}}$
$\displaystyle=\int\mathbf{n}\cdot\hat{\mathrm{\bm{\sigma}}}\;\mathrm{d}\Gamma_{p}=\mathbf{0},$
(17a) $\displaystyle\hat{\mathbf{T}}$
$\displaystyle=\int\left(\mathbf{x}_{s}-\mathbf{x}_{c}\right)\times\left(\mathbf{n}\cdot\hat{\mathrm{\bm{\sigma}}}\right)\mathrm{d}\Gamma_{p}=\mathbf{0}.$
(17b)

Equations (13)-(17) define an eigenvalue problem. The stability of the base
state is determined by the eigenvalue with the largest real part
$\sigma_{r}^{0}$, namely the leading eigenvalue, and the corresponding
perturbations $\left(\hat{\mathbf{u}},\hat{p},\hat{c}\right)^{T}$ are called
leading eigenmodes. The base state is stable when $\sigma_{r}^{0}<0$ but
unstable when $\sigma_{r}^{0}>0$. The Péclet number at which
$\sigma_{r}^{0}=0$ is precisely the critical Péclet number
$\mathrm{Pe}_{c}^{(1)}$ signifying the transition from a stationary state to
steady propulsion. Unless stated otherwise, the eigenvalues mentioned below
refer to the leading eigenvalues, and the superscript $0$ is omitted for
simplicity. We solve the eigenvalue problem with COMSOL using the eigenvalue
solver ARPACK. The validation of our approach is demonstrated in Appendix A.

##  IV Results

###  IV.1 Diverse behaviours of an elliptic phoretic disk

By increasing Péclet number $\mathrm{Pe}$, we demonstrate the diverse
$\mathrm{Pe}$-dependent behaviors of a disk with the eccentricity $e=0.87$.
For a sufficiently small $\mathrm{Pe}$ below a critical value
$\mathrm{Pe}_{c}^{(1)}$, the disk undergoes transient rotation before
recovering a stationary state. When $\mathrm{Pe}$ goes above
$\mathrm{Pe}_{c}^{(1)}$, e.g., $\mathrm{Pe}=0.4$, the stationary state
transits into the steady propulsion, in which the fore-aft symmetry in the
concentration profile is broken (see Figure 2f), and the resultant
concentration polarity induces a rectilinear motion with a constant swimming
speed, as depicted in Figure 2(a). Further increasing $\mathrm{Pe}$ to $0.62$,
the directed motion loses its stability, leading to a secondary instability
characterized by spontaneous chiral symmetry breaking. Consequently, the disk
repeatedly executes a circular path, accompanied by its own clockwise
rotation, as indicated in Figure 2(b). This regime featuring circular
trajectory and self-rotation is called the orbiting regime. Hu et al. [27] and
Li [42] observed analogous circular trajectories executed by circular swimmers
without self-rotation. At $\mathrm{Pe}=10$, the elliptic disk favors a wave-
like trajectory after a transient period of swimming straightforward (see
Figure 2c). From the viewpoint of the frame co-moving with the disk center,
the disk periodically swings like a pendulum (refer to Figure 2g). This
swimming behavior is analogous to that exhibited by an active prolate double-
core droplet, which moves along an undulatory trajectory [25]. Li [42] also
reported that a high-$\mathrm{Pe}$ active drop swims along a periodic zigzag
trajectory. As $\mathrm{Pe}$ grows to $25$, the elliptic disk surprisingly
recovers to steady propulsion, as illustrated in Figure 2(d). In
correspondence, chiral symmetry recovers. Unlike the ballistic motion after
rotation at low $\mathrm{Pe}$, the disk here travels straight from the onset
of instability without any rotation. The steady propulsion becomes unstable at
a higher $\mathrm{Pe}$, e.g., $\mathrm{Pe}=33$, accordingly, the disk swims
straightforward at an oscillating speed (see Figure 12 in Appendix B). Hu et
al. [28] and Kailasham & Khair [34] identified the similar motion of an
isotropic phoretic particle under the axisymmetric assumption. Figure 2(e)
indicates that the disk enters into a chaotic regime characterized by an
erratic trajectory at $\mathrm{Pe}=60$. In contrast to frequent intermittency
and random walk occurring in the chaotic regime for a circular disk, as
observed by Hu et al. [27], the change of velocity in both direction and
magnitude experienced by the elliptic disk is not drastic, yielding a less
chaotic trajectory.

Figure 2: An elliptic disk of eccentricity $e=0.87$ follows typical
trajectories depending on $\mathrm{Pe}$: (a) steady, (b) orbiting, (c)
periodic, (d) steady, and (e) chaotic. The color of trajectories is coded by
the time $t$. Green and red arrows denote the direction of the translational
velocity $\mathbf{U}$ and the disk orientation $\mathbf{e}_{s}$, respectively.
(f) Polarized concentration distribution with respect to $\mathbf{e}_{s}$ at
$\mathrm{Pe}=0.4$. (g) Periodic pendulum-like swinging of the disk at
$\mathrm{Pe}=10$ in the frame co-moving with its center. The star and triangle
symbols marked in (g) denote the two distinct moments when the disk reaches
the peak and trough of its trajectory, respectively, as shown in (c).

Having observed distinct swimming behaviours of a disk with a specific
eccentricity at varying $\mathrm{Pe}$, we then systematically examine how the
eccentricity $e$ affects these $\mathrm{Pe}$-dependent behaviours. The phase
diagram in Figure 3 shows that shifts in $\mathrm{Pe}$-dependent behaviours of
an elliptic disk with $e$ below $0.75$, resemble those of a circular
counterpart ($e=0$). The elliptic disk successively executes stationary,
steady, periodic, and chaotic motion with increasing $\mathrm{Pe}$. As $e$
exceeds $0.75$, the disk exhibits richer dynamics: the orbiting motion emerges
with broken chiral symmetry. In fact, we have also explored the scenarios at
$e>0.9$, e.g., $e=0.96$, and observed that the swimming dynamics of the
elliptic disk are almost dominated by orbiting and chaotic regimes (not shown
here). These two regimes have been revealed in the current phase diagram,
hence the characteristic locomotory modes can be well captured in the range of
$e$ considered.

We further discern from the phase diagram that the variation of $e$
categorizes certain identified swimming patterns into two types. In the
stationary or steady regime, the aforementioned observations indicate that the
disk experiences transient rotation before recovering a stationary state or
retaining a steady motion (see Figure 2). Nevertheless, these situations are
only applicable at a higher $e$. When $e<0.75$, the disk remains stationary or
swims straight without any rotation. Besides, two types of periodic motions
are depicted in the phase diagram: 1) a disk swings along a wavy trajectory
(light-colored squares); 2) a disk swims straight at an oscillating swimming
speed (dark-colored squares). The former and the latter are termed swinging
and straight periodic motions, respectively.

Figure 3: Phase diagram characterizing the behaviors of an elliptic phoretic
disk depending on its eccentricity $e$ and the Péclet number $\mathrm{Pe}$. It
shows five regimes: stationary, steady, orbiting, periodic, and chaotic. The
solid line denotes the LSA prediction.

###  IV.2 Spontaneous steady propulsion triggered by instability

Upon gaining a general understanding of the phase diagram, we then explore the
detailed swimming dynamics within each regime. We first focus on the
transition from the stationary to steady regime. The disk is stationary in the
base state. As $\mathrm{Pe}$ exceeds $\mathrm{Pe}_{c}^{(1)}$, an instability
arises, and the disk sets into steady propulsion along its major axis. Here,
we perform a global LSA, as introduced in Section III, to quantitatively
identify $\mathrm{Pe}_{c}^{(1)}$ and its dependence on the disk shape.
$\mathrm{Pe}_{c}^{(1)}$ predicted by the LSA monotonically decreases with $e$,
as depicted by the solid line in the above phase diagram, implying that the
ellipticity of the disk promotes its instability. Figure 4(a) depicts the
concentration field $\hat{c}$ of the global eigenmode at
$\mathrm{Pe}_{c}^{(1)}=0.43$ for $e=0.55$. We see that the symmetry of the
base concentration field $\hat{c}$ is broken in the direction of the major
axis at $\mathrm{Pe}_{c}^{(1)}$. The fore-aft asymmetric concentration
distribution induces a downward slip flow, as shown by the flow field
$\hat{\mathbf{u}}$ of the eigenmode in figure 4(b), driving the steady motion
of the disk along its major axis. In close proximity to
$\mathrm{Pe}_{c}^{(1)}$, the swimming speed $U_{mg}$ is proportional to
$\sqrt{\mathrm{Pe}-\mathrm{Pe}_{c}^{(1)}}$ (see Figure 4c), implying that the
instability occurs through a supercritical pitchfork bifurcation. This is
parallel to the observation of Hu et al. [27] and Li [42].

We further probe the physical mechanism underlying the instability by adopting
an approach resembling the energy budget analysis [3, 1]. By taking an inner
product denoted by $\left(\cdot,\cdot\right)$) of the equation (14) with
$\hat{c}$ in the $L^{2}$ space, we arrive

$\left(\sigma\hat{c},\hat{c}\right)+\left(\hat{\mathbf{u}}\cdot\nabla
c_{b},\hat{c}\right)=\frac{1}{\mathrm{Pe}}\left(\Delta\hat{c},\hat{c}\right).$
(18)

Using integration by parts, the divergence-free condition in (13), and the
boundary condition (15), we derive

$\displaystyle\left(\hat{\mathbf{u}}\cdot\nabla c_{b},\hat{c}\right)$
$\displaystyle=-\left(c_{b},\hat{\mathbf{u}}\cdot\nabla\hat{c}\right),$ (19a)
$\displaystyle\frac{1}{\mathrm{Pe}}\left(\Delta\hat{c},\hat{c}\right)$
$\displaystyle=-\frac{1}{\mathrm{Pe}}\parallel\nabla\hat{c}\parallel^{2},$
(19b)

with $\parallel\cdot\parallel$ denoting the $L^{2}$-norm. By substituting (19)
into (18), we obtain

$\sigma_{e}=\sigma_{c}+\sigma_{d},$ (20)

with

$\displaystyle\sigma_{c}$
$\displaystyle=\frac{\left(c_{b},\hat{\mathbf{u}}\cdot\nabla\hat{c}\right)}{\left(1,\hat{c}^{2}\right)},$
(21a) $\displaystyle\sigma_{d}$
$\displaystyle=-\frac{\parallel\nabla\hat{c}\parallel^{2}}{\left(\mathrm{Pe},\hat{c}^{2}\right)}.$
(21b)

Here, the growth rate $\sigma_{e}$ is introduced in (20) to be distinguished
from $\sigma$ obtained by the LSA. $\sigma_{c}$ and $\sigma_{d}$ represent the
contributions of advection and diffusion to $\sigma_{e}$, respectively. For a
steady motion, the imaginary part of the growth rate vanishes, thus, the
growth rate only has its real part, e.g., $\sigma=\sigma_{r}$. Figure 4(c)
shows that $\sigma_{e}$ and $\sigma$ lie on top of each other, giving us the
confidence to analyze the dominant physical ingredient that drives the
instability using (20). We naturally infer from (20) that $\sigma_{c}$ is
responsible for $\sigma_{e}$ turning positive by realizing that $\sigma_{d}$
is consistently negative, as confirmed by figure 4(c). Hence, as anticipated,
advection drives the instability and diffusion dampens the perturbation, the
balance between them dominates the phoretic dynamics of the elliptic disk: at
small $\mathrm{Pe}$, diffusion dominates and its homogenizing effect maintains
a fore-aft symmetric solute distribution. As $\mathrm{Pe}$ grows beyond
$\mathrm{Pe}_{c}^{(1)}$, advection suppresses diffusion and amplifies the
asymmetric solute disturbance. The slip flows triggered by the asymmetric
concentration distribution drive the disk to swim spontaneously.

Figure 4: Global eigenmodes at the critical Péclet number
$\mathrm{Pe}_{c}^{(1)}=0.43$ for an elliptic disk with $e=0.55$. The eigenmode
is characterized by (a) the perturbation concentration $\hat{c}$, and (b) the
perturbation velocity $\hat{\mathbf{u}}$ (red arrows) and its $y$ component
$\hat{v}$ (color map). (c) Dependence of time-averaged disk speed
$\left<U_{mg}\right>$ on $\mathrm{Pe}-\mathrm{Pe}_{c}^{(1)}$ at varying $e$.
In the vicinity of $\mathrm{Pe}_{c}^{(1)}$, $\left<U_{mg}\right>$ is
proportional to $\sqrt{\mathrm{Pe}-\mathrm{Pe}_{c}^{(1)}}$. (d)
$\mathrm{Pe}$-dependent growth rate $\sigma$ based on the LSA and that
$\sigma_{e}=\sigma_{c}+\sigma_{d}$ derived from the concentration perturbation
equation (14). The latter comprises the contributions $\sigma_{c}$ and
$\sigma_{d}$ from advection and diffusion, respectively.

Next, we analyze the steady propulsion of the autophoretic disk. We show in
figure 5(a) that the elliptic disk of $e=0.55$ swims as a puller attracting
the fluid from its front and rear at $\mathrm{Pe}=0.5$. In contrast, it
becomes a neutral-type swimmer at $\mathrm{Pe}=1.5$, as indicated in figure
5(b). In retrospect, Suda et al. [72] reported an analogous transition from a
puller swimming straight into a pusher-type swimmer executing unsteady
curvilinear motion when $\mathrm{Pe}$ increases. In that case, the droplet
motion is triggered by the concentration gradient caused by a point source of
surfactant at the surface. Besides, Li [42] observes that an active drop
transits from a steady pusher to a mixed pusher-puller propelling unsteadily
as $\mathrm{Pe}$ grows. For these droplet swimmers, the switching of their
disturbance flow occurs as they go from steady to unsteady motion.
Nevertheless, the disk swimmer here exhibits the puller-to-neutral transition
without a shift in its swimming patterns. We further examine the flow field
around a disk swimmer throughout the whole steady regime and summarize these
results in a phase diagram in figure 5(c). It is found that a circular ($e=0$)
or nearly-circular disk can only be a neutral swimmer, and the puller-to-
neutral transition emerges as $e\geq 0.20$. The critical Péclet number
$\mathrm{Pe}_{t}$ signifying this transition depends monotonically on $e$, as
indicated by the dashed line. When $e>0.65$, the transition disappears and the
disk solely swims as a puller. In fact, this peculiar transformation can be
understood by examining the solute distribution at the disk surface, as
discussed below. For visualization purpose, we normalize the concentration via
$\tilde{c}=\gamma\left(c_{max}-c\right)/\left(c_{max}-c_{min}\right)$, where
$c_{max}$ and $c_{min}$ denote the maximum and minimum concentrations at the
disk surface, respectively. The coefficient $\gamma$, e.g., $\gamma=0.5$,
restricts the normalized concentration $\tilde{c}$ in the range of
$[0,\gamma]$.

Figure 5: (a) A puller-type steady swimmer at $\mathrm{Pe}=0.5$ transitions to
(b) a neutral-type counterpart at $\mathrm{Pe}=1.5$, where the eccentricity
$e$ is $0.55$. The color maps indicate the distribution of $c$. The stagnation
point $\mathrm{P}$ coincides with the peak of $c$ at the disk surface. (c)
$\mathrm{Pe}-e$ phase diagram shows puller-type and neutral-type steady
swimmers demarcated by the dashed line. Solid lines separate regimes
identified in figure 3.  Figure 6: Normalized concentration at the surface of
a disk swimming steadily: (a)-(c) for an elliptic disk with $e=0.55$ at
$\mathrm{Pe}=0.3$, $0.5$, and $1.5$, respectively; (e)-(f) for a circular disk
at $\mathrm{Pe}=0.3$ and $0.5$, respectively. $c_{max}$ denotes the peak
concentration at the disk surface and $\beta$ is the polar angle with respect
to the disk orientation $\mathbf{e}_{s}$. The stagnation point $\mathrm{P}$
coincides with the location of $\tilde{c}_{max}$ at $\mathrm{Pe}=0.5$. (d):
transition from a bimodal to polarized concentration distribution at the
surface of an elliptic disk by increasing $\mathrm{Pe}$. (e): similar to (d)
but for a circular disk with a transition from an isotropic solute
distribution.

We demonstrate in figure 6(a) and (d) the bimodal distribution of normalized
concentration caused by the curvature variation of the disk surface in the
base state ($\mathrm{Pe}=0.3$). The maximum normalized concentration
$\tilde{c}_{max}$ is located at the left vertex of the minor axis
corresponding to $\beta=90^{\circ}$, where $\beta$ denotes the polar angle
with respect to $\mathbf{e}_{s}$. Here, we consider only the left half of the
disk ($\beta\in\left[0,180^{\circ}\right]$) for its symmetry. At
$\mathrm{Pe}=0.5$, the base state loses stability and solutes are advected
toward the disk’s rear ($\beta=180^{\circ}$). Correspondingly, the position of
$c_{max}$ migrates rearward from $\beta=90^{\circ}$ to the stagnation point
$\mathrm{P}$, as depicted in figure 6(b). Considering that the slip velocity
$\mathbf{u}_{s}=\nabla_{s}c$ is directed from low to high concentration at the
disk surface, the fluid is attracted from the front and rear of the elliptic
disk to the stagnation position $\mathrm{P}$, forming a puller-type swimmer
(see figure 5a). As $\mathrm{Pe}$ increases to $1.5$, the enhanced convective
transport of solutes causes the shift from a bimodal to a polarized
concentration profile with the peak concentration located at
$\beta=180^{\circ}$ (see figure 6c). Correspondingly, the slip flow induced by
the concentration gradient is driven from the front to the rear of the disk
(see figure 6b), generating a prototypical neutral-type swimmer. We comment
that as an elliptic disk swims steadily, a bimodal concentration profile at
its surface yields a puller-type swimmer, while a polarized one leads to a
neural-type counterpart. Different from the elliptic disk, a circular disk has
an isotropic base distribution of the solute ($\mathrm{Pe}=0.3$), as shown in
figure 6(e) and (g). For an unstable base state at $\mathrm{Pe}=1.5$, solutes
are advected rearward along the disk surface, leading to the direct transition
from an isotropic to a polarized concentration profile. Hence, the circular
disk behaves as a neutral swimmer solely.

###  IV.3 Chiral symmetry-breaking orbiting motion

Having observed the puller-neutral transition in the steady regime, we direct
our focus on the chiral symmetry-breaking orbiting motion of the elliptic
disk. Here, the disk approximately swims along a circular trajectory (see
figure 2b) with the rotational velocity $\Omega$ oscillating around a constant
value. We find that the translational velocity $\mathbf{U}$ (green arrow)
deviates from $\mathbf{e}_{s}$ (red arrow) by an angle $\alpha$. The evolution
of $\alpha$ resembles that of $\Omega$ oscillating slightly around
$15^{\circ}$ for $\mathrm{Pe}=4$, as depicted in figure 7(a). In fact, the
misalignment $\mathbf{U}$ and $\mathbf{e}_{s}$ can be rationalized by
examining the solute distribution $c(x,y)$ around the disk. Figure 7(b)
illustrates $c(x,y)$ and streamlines at the moment denoted by a hexagon in
Figure 7(a). The left-right asymmetry of the solute distribution about
$\mathbf{e}_{s}$ gives rise to the asymmetric slip velocity, as reflected by
the streamlines, hence causing the inconsistency in directions of $\mathbf{U}$
and $\mathbf{e}_{s}$. Besides, we notice that three stagnation points near the
rear of the disk in figure 7(b) gradually migrate toward the rear with
increasing $\mathrm{Pe}$ in the orbiting regime (not shown here). Accordingly,
the time-averaged $\left<\alpha\right>=\int_{0}^{T}\alpha\mathrm{d}t/T$ within
a time window $T$ decreases monotonically with $\mathrm{Pe}$, as shown in the
inset of figure 7(a). Intuitively, we infer $\mathbf{U}$ is aligned with
$\mathbf{e}_{s}$ ($\alpha=0$) as three stagnation points exactly coincide at
the rear of the disk, and indeed they do as $\mathrm{Pe}$ grows beyond $14$
where the disk executes a steady motion, e,g., $\mathrm{Pe}=25$ (streamlines
and concentration distribution are analogous to Figure 5b).

Figure 7: (a) Time evolution of the angle $\alpha$ between the translational
velocity $\mathbf{U}$ and the disk orientation $\mathbf{e}_{s}$ as the
elliptic disk with $e=0.87$ approximately follows a circular trajectory for
$\mathrm{Pe}=4$. The inset shows the monotonic decrease of the time-averaged
$\alpha$ with $\mathrm{Pe}$. (b) Solute distribution and streamlines at the
instant indicated by a hexagon in (a).

The orbiting regime is characterized by the continuous rotation of an elliptic
disk, which loses its chiral symmetry spontaneously. Conversely, a circular
disk maintains this symmetry irrespective of the $\mathrm{Pe}$ value. In
particular, for $\mathrm{Pe}\in[9,13]$, despite the asymmetric solute
distribution and streamlines resembling those in figure 7(b) occasionally, the
circular disk only executes a meandering motion [27] without rotation.
Ideally, a secondary stability analysis similar to the one performed above
will help understand the mechanism for the shape-induced chiral symmetry
breaking. This task is however technically challenging and will be pursued in
the future.

Notably, experiments on active droplets in a Hele-Shaw cell showed that the
droplet preferred to avoid the concentration trail emitted by itself at
earlier times, termed a self-avoiding walk phenomenon [26]. The elliptic disk
here does not avoid similarly but instead repeats its previous trajectory, as
shown in figure 2(b). The difference can be rationalized by a scaling
analysis. We estimate that the time $\bar{t}_{d}$ required for the
concentration to decay to zero is of the order of $\bar{b}^{2}/\bar{D}$. The
disk moves with the speed $U_{mg}\bar{\mathcal{V}}$, hence the traveling time
$\bar{t}_{p}$ it takes to execute a circular trajectory of perimeter
$\bar{L}_{p}$ is $\bar{L}_{p}/\left(U_{mg}\bar{\mathcal{V}}\right)$
approximately. Recalling the definitions of $\bar{\mathcal{V}}$ and
$\mathrm{Pe}$ in Section II, we rewrite $\bar{t}_{p}$ as
$\bar{L}_{p}\bar{b}/\left(U_{mg}\mathrm{Pe}\bar{D}\right)$. The ratio between
these two time scales is

$\displaystyle\frac{\bar{t}_{p}}{\bar{t}_{d}}=\frac{\bar{L}_{p}/\bar{b}}{U_{mg}\mathrm{Pe}}.$
(22)

At a specific $\mathrm{Pe}$, e.g., $\mathrm{Pe}=4$, in the orbiting regime,
$\bar{t}_{p}/\bar{t}_{d}\approx 87.2\gg 1$. This substantial ratio indicates
that the disk swimmer can barely sense and thus avoid its own rapidly decaying
chemical trace. In contrast, for a sufficiently larger $\mathrm{Pe}$, we may
observe the self-avoidance due to the significantly reduced travel time
$t_{p}$ [28, 26].

###  IV.4 Periodic swinging

The orbiting motion typically transitions to periodic swinging (refer to
figure 2c) as $\mathrm{Pe}$ increases, which will be analyzed below. Figure
8(a) illustrates that the time-evolving rotational velocity $\Omega$ of the
disk with $e=0.87$ at $\mathrm{Pe}=12.5$ is characterized by two phases.
First, $\Omega$ grows rapidly due to self-oscillation, and the dashed line
connecting local peaks $\Omega_{pk}$ of $\Omega$ indicates an exponential
growth of $\Omega$ in time. This trend is confirmed by the linear relationship
between $\mathrm{log}\;\Omega_{pk}$ and $t$ shown in the inset of figure 8(a).
Second, $\Omega$ saturates nonlinearly to a periodic state with a constant
amplitude $\Omega_{mg}$. The sinusoidal-like variation of $\Omega$ leads to a
wave-like trajectory. We plot $\Omega_{mg}^{2}$ as a function of $\mathrm{Pe}$
near the critical $\mathrm{Pe}\approx 13.5$ (star) in figure 8(b). The linear
dependence of $\Omega_{mg}^{2}$ on $\mathrm{Pe}$ implies that the steady
motion loses stability through a Hopf bifurcation. The critical $\mathrm{Pe}$
signifies the boundary between the periodic swinging and steady motion. We
further present in Figure 8(c) the phase portrait in the $\Omega-U$ plane,
where the unstable spiral at the origin grows continuously to an elliptic
stable limit circle, indicating the supercritical nature of the Hopf
bifurcation. Notably, when $e<0.87$, we observe a transition from steady to
periodic swinging with increasing $\mathrm{Pe}$, e.g., $e=0.81$ (see Figure
3). The Hopf bifurcation still holds for the onset of instability that
triggers this transition.

Figure 8: (a) Time evolution of the rotational velocity $\Omega$ of an
elliptic disk with $e=0.87$ at $\mathrm{Pe}=12.5$. The dashed line denotes
certain local peaks $\Omega_{pk}$ of $\Omega$, and the inset shows the linear
dependence of $\mathrm{log}\;\Omega_{pk}$ on $t$. (b) Linear variation of
$\Omega_{mg}^{2}$ in $\mathrm{Pe}$ near the critical $\mathrm{Pe}\approx 13.5$
(star), where $\Omega_{mg}$ denotes the constant amplitude of $\Omega$ at
$t>5400$, as depicted by the dotted line in (a). The periodic motion recovers
to steady propulsion as $\mathrm{Pe}$ grows beyond the critical value. (c)
Phase portrait in the $\Omega-U$ plane, with the color of the unstable spiral
coded by $t$.  Figure 9: (a) Mean square displacement (MSD) and velocity
autocorrelation function (VAF) for disks with different shapes at
$\mathrm{Pe}_{c}^{(2)}+10$, where $\mathrm{Pe}_{c}^{(2)}$ is the shape-
dependent critical Péclet number corresponding to the onset of chaos. (b) The
exponent $\lambda$ of the power-law scaling $\text{MSD}\propto\tau^{\lambda}$
versus $e$. (c) Monotonic dependence of $\lambda$ on $\mathrm{Pe}$ for the
disk of $e=0.2$. (d) Chaotic trajectories (color-coded by $t$) followed by an
elliptic disk of $e=0.2$ at $\mathrm{Pe}=8$, $21$ and $47$, respectively.

###  IV.5 Chaotic swimming dynamics

Finally, we analyze the chaotic motion of elliptic disks by examining their
mean square displacement (MSD) and the velocity autocorrelation function
(VAF). We calculate the $\text{MSD}\left(\tau\right)$ of a disk [50] that
depends on a time lag $\tau$ by

$\text{MSD}\left(\tau\right)=\frac{1}{T-\tau}\int_{0}^{T-\tau}\left[\mathbf{r}\left(t+\tau\right)-\mathbf{r}\left(t\right)\right]^{2},$
(23)

where $\mathbf{r}\left(t\right)$ denotes the time-dependent displacement of a
disk relative to its original position. Clearly, combining a large time period
$T$ and $\tau\ll T$ can improve the statistical confidence. We calculate VAF
using

$\text{VAF}\left(t\right)=\frac{1}{T}\int_{0}^{T}\frac{\mathbf{U}\left(t\right)\cdot\mathbf{U}\left(t+\Delta
t\right)}{\left|\mathbf{U}\left(t\right)\right|\left|\mathbf{U}\left(t+\Delta
t\right)\right|}\mathrm{d}t,$ (24)

where a sufficiently large $T$ is used to attain statistical invariance [9].

We examine how the eccentricity $e$ of a disk affects its chaotic motion.
First, we determine the shape-dependent $\mathrm{Pe}_{c}^{(2)}(e)$ when chaos
emerges. Then, we investigate the chaotic dynamics at
$\mathrm{Pe}(e)=\mathrm{Pe}_{c}^{(2)}(e)+10$ above that threshold by a fixed
offset, namely, $10$ here. The corresponding MSD and VAF are depicted in
figure 9(a). At early times, despite the diversity in the disk shape, they all
swim persistently resulting in the quadratically growing MSD in $\tau$. At
long times, the disk shape considerably affects its phoretic movement, leading
to a transition from a random-walking circular disk to a ballistically
swimming elliptic counterpart. For the former, we reproduce perfectly its
random walk behavior as reported by Hu et al. [27], Lin et al. [44], Hu et al.
[28], Hokmabad et al. [26], which features the linear scaling
$\text{MSD}\propto\tau$ and the decorrelation in velocity $\mathbf{U}$ due to
the rapidly changing swimming direction (see the inset). As the eccentricity
$e$ increases to $0.1$, $\text{MSD}\varpropto\tau^{1.5}$, reminiscent of the
self-avoidance walk identified experimentally [26]. We wonder whether this
growth scaling, $\propto\tau^{1.5}$, results from the self-avoidance walk. By
utilizing the scaling analysis mentioned above, we compare the decay time of
the chemical $\bar{t}_{d}$ and the traveling time $\bar{t}_{p}$ of the disk
and find that $\bar{t}_{p}/\bar{t}_{d}\approx 40.8\gg 1$ (see Equation 22),
precluding self-avoidance here. At $e=0.87$, the MSD’s quadratic growth over
lag time implies that the elliptic disk approximately executes the ballistic
motion, as confirmed by the trajectory shown in figure 2(e). Accordingly, the
VAF is almost constant in time, indicating a correlation in velocity
$\mathbf{U}$. It is worth noting that Hu et al. [28] and Morozov & Michelin
[54] also found the quadratic variations of MSD for the chaotic swimming of a
phoretic particle and an active droplet, respectively, under the axisymmetric
assumption.

By probing further in figure 9 the exponent $\lambda$ of power-law scaling
$\text{MSD}\propto\tau^{\lambda}$ at varying $e$, we notice that the normal
diffusion ($\lambda=1$) transitions to superdiffusion ($\lambda>1$) as $e$
increases. This observation demonstrates that the shape of a disk can
significantly affect its diffusion behavior. In addition, we explore the
effect of $\mathrm{Pe}$ on the chaotic motion of an elliptic disk, as depicted
in figure 9(c). In contrast to the eccentricity $e$, $\lambda$ monotonically
decreases with $\mathrm{Pe}$, implying the transition from superdiffusion to
normal diffusion. Accordingly, the elliptic disk undergoes a ballistic motion
near $\mathrm{Pe}_{c}^{(2)}$ and swims randomly at a larger $\mathrm{Pe}$,
e.g., $\mathrm{Pe}=47$ (see figure 9d). An analogous dependence of $\lambda$
on $\mathrm{Pe}$ for the chaotic motion of an axisymmetric spherical particle
was also reported by Kailasham & Khair [34];

##  V Conclusions and discussions

In this work, we numerically and theoretically investigate the swimming
dynamics of an elliptic disk uniformly releasing solutes in the creeping flow
regime. The disk with an eccentricity $e<0.75$ mimics a circular counterpart:
it shows stationary, steady, periodic, and chaotic behaviors dependent on
$\mathrm{Pe}$. When $e>0.75$, the disk attains an orbiting motion via an
instability spontaneously breaking the chiral symmetry.

By performing a global LSA, we theoretically predict the critical Péclet
number $\mathrm{Pe}_{c}^{(1)}$ above which a stationary disk becomes a steady
swimmer spontaneously triggered by instability. The ellipticity of the disk is
found to promote this instability. Besides the LSA, we calculate the growth
rate of perturbation via a method akin to energy budget analysis. The results
agree with the LSA predictions and provide new insight into the competing
roles of advection and diffusion.

We observe that the transition from a puller-type to a neutral-type steady
disk swimmer is induced by a stronger rearward advection of solute due to
growing $\mathrm{Pe}$. Accordingly, a bimodal concentration profile
corresponding to a puller-type swimmer becomes a polarized profile leading to
a neutral-type swimmer. The orbiting disk repeatedly swims along a circular
trajectory while simultaneously rotating.

Two distinct types of periodic motions are identified: swinging and straight
periodic motions. The former develops from a steady motion through a
supercritical Hopf bifurcation and features a wave-like trajectory. The latter
is characterized by a rectilinear motion with an oscillating swimming speed.
Finally, the effects of the disk shape and $\mathrm{Pe}$ on the chaotic motion
are examined. We uncover a shift from normal diffusion to superdiffusion with
growing eccentricity $e$: the former and latter correspond to a random-walking
circular disk and a ballistically swimming elliptic counterpart, respectively.
The influence of $\mathrm{Pe}$ on the disk’s diffusion behavior stands in
contrast to that of $e$.

It is important to note that both the elliptic camphor disk [36] and the
prolate composite droplet [25] observed in laboratory studies demonstrate
self-propulsion along their minor axes. In contrast, the elliptic phoretic
disk in this work swims steadily along its major axis. The differences in
these cases may be due to the failure of 2D simulations to capture the complex
3D physicochemical hydrodynamics present in experiments. Specifically, a
swimming camphor disk is driven by Marangoni flow at the air-liquid interface
and the resulting 3D subsurface flow, while the prolate droplet self-propels
within a Hele-Shaw cell. We are planning 3D studies to address these
discrepancies.

Figure 10: Validation of our numerical implementation against published works.
The numerical and published data are represented by markers and lines,
respectively. (a) Swimming speed of an isotropic autophoretic disk in a
circular domain of radius $R=200$ for varying $\mathrm{Pe}$, benchmarked by
[27]. Note that the time-averaged swimming speed
$\left<U_{mg}\right>=\int_{0}^{T}U_{mg}\;\mathrm{d}t$ within a time period $T$
recovers $U_{mg}$ for steady propulsion. (b) Swimming speed of an autophoretic
spherical particle versus $\mathrm{Pe}$ computed in an axisymmetric
configuration, in comparison to Michelin et al. [52]; the inset shows the
growth rate $\sigma$ of the unstable eigenmode versus $\mathrm{Pe}$. (c)
Swimming velocity of a spheroidal Janus particle as a function of $\zeta$,
validated against Popescu et al. [65]. $\zeta$ denotes the height of the
border dividing the particle into active and inert (marked in grey)
compartments, with $\zeta=0$ corresponding to the particle center. The
eccentricity $e$ of the spheroid recovers to zero for a spherical particle.

##  Appendix A Validation of numerical implementations

We show the validation of our COMSOL implementations for the numerical
simulation (Equations (3)-(8)) and stability analysis (Equations (13)-(17)).
First, we study the inertialess self-propulsion of an isotropic autophoretic
disk in a circular domain of radius $R=200$. This setup has been investigated
by Hu et al. [27] combining a stability analysis and simulations based on a
spectral method. The disk transits to a steady swimmer from its stationary
state at the first threshold $\mathrm{Pe}\approx 0.466$. Increasing
$\mathrm{Pe}$, the swimmer becomes unstable at the second critical condition
$\mathrm{Pe}\approx 4.65$, moving in a meandering manner. Our numerical data
shown in Figure 10(a) exactly recovers the reported $\mathrm{Pe}$-dependent
swimming speed and two critical $\mathrm{Pe}$ values. Besides the 2D case, we
examine the spontaneous motion of an axisymmetric isotropic particle in an
unbounded domain. In this case, instability occurs at $\mathrm{Pe}=4$ as
predicted by Michelin et al. [52]. Our data on swimming speed agrees well with
the numerical results therein, as depicted in Figure 10(b). Also, its inset
shows that the present LSA recovers the threshold $\mathrm{Pe}=4$ and the
eigenvalues depending on $\mathrm{Pe}$. As a side product, we probe the effect
of inertia $\mathrm{Re}$ on this scenario. The inertial effect systematically
enhances the swimming speed, becoming pronounced when $\mathrm{Re}>1$ (see
Figure 11a). The inertia-induced relative enhancement of the swimming speed
$\varepsilon$ in Figure 11(b) is observed to linearly scale with $\mathrm{Re}$
in this weak inertia regime, which deserves further theoretical underpinning.

Finally, we examine the ability of our implementation in handling non-
spherical phoretic swimmers. We focus on a spheroidal Janus colloid consisting
of active and inert compartments. Its propulsion depends on its eccentricity
$e$ and $\zeta$, where $\zeta$ denotes the height of the border separating the
two compartments. The obtained swimming speed agrees well with the reference
solutions over a wide range of $e$ and $\zeta$ [65], as shown in Figure 10(c).

Figure 11: (a) Swimming speed $U_{mg}$ of an isotropic phoretic particle
versus $\mathrm{Pe}$ when $\mathrm{Re}$ varies, which is calculated in an
axisymmetric setup. Our numerical data are compared with that of Michelin et
al. [52] at $\mathrm{Re}=0$. (b) Inertia enhances the swimming speed,
characterized by the linear relation between the relative enhancement
$\varepsilon=\left[U_{mg}(\mathrm{Re},\mathrm{Pe})-U_{mg}(\mathrm{Re}=0,\mathrm{Pe})\right]/U_{mg}(\mathrm{Re}=0,\mathrm{Pe})$
and $\mathrm{Re}$.

##  Appendix B Periodic motion of an elliptic phoretic disk

We identify two types of periodic motion: the first showcases a wave-like
trajectory caused by the periodic rotation of the elliptic disk (see Figure
2g); the second corresponds to unidirectional rectilinear motion with a time-
periodic swimming speed, as shown in Figure 12.

Figure 12: Straight periodic motion of an elliptic disk with $e=0.87$ at
$\mathrm{Pe}=33$. (a) Straight trajectory color-coded by time $t$. The
swimming direction $\mathbf{e}_{s}$ (red arrow) coincides with the
translational velocity $\mathbf{U}$ (green arrow). (b) Time evolution of the
swimming speed $U_{mg}$.

###### Acknowledgements.

We thank Qianhong Yang for helpful discussions and Wei-Fan Hu for generously
sharing the data utilized in the validation. L.Z. thanks Singapore Ministry of
Education Academic Research Fund Tier 2 (MOE-T2EP50221-0012 &
MOE-T2EP50122-0015) and Tier 1 (A-8000197-01-00) grants, and the Paris-NUS
joint research grant (ANR-18-IDEX-0001 & A-0009528-01-00). Some computation of
the work was performed on resources of the National Supercomputing Centre,
Singapore.

## References

  * Abubakar & Matar [2022] Abubakar, HA & Matar, OK 2022 Linear stability analysis of taylor bubble motion in downward flowing liquids in vertical tubes. J. Fluid Mech. 941, A2.
  * Ban et al. [2018] Ban, Takahiko, Sugiyama, Michiaki, Nagatsu, Yuichiro & Tokuyama, Hideaki 2018 Motion-based detection of lanthanides (III) using self-propelled droplets. J. Phys. Chem. B 122 (46), 10647–10651.
  * Bendiksen [1985] Bendiksen, Kjell H 1985 On the motion of long bubbles in vertical tubes. Int. J. Multiphas. Flow 11 (6), 797–812.
  * de Blois et al. [2019] de Blois, Charlotte, Reyssat, Mathilde, Michelin, Sébastien & Dauchot, Olivier 2019 Flow field around a confined active droplet. Phys. Rev. Fluids 4 (5), 054001\.
  * Boniface et al. [2021] Boniface, Dolachai, Cottin-Bizonne, Cécile, Detcheverry, François & Ybert, Christophe 2021 Role of marangoni forces in the velocity of symmetric interfacial swimmers. Phys. Rev. Fluids 6 (10), 104006.
  * Brady [2011] Brady, John F 2011 Particle motion driven by solute gradients with application to autonomous motion: continuum and colloidal perspectives. J. Fluid Mech. 667, 216–259.
  * Bunea & Glückstad [2019] Bunea, Ada-Ioana & Glückstad, Jesper 2019 Strategies for optical trapping in biological samples: Aiming at microrobotic surgeons. Laser Photonics Rev. 13 (4), 1800227\.
  * Campbell et al. [2019] Campbell, Andrew I, Ebbens, Stephen J, Illien, Pierre & Golestanian, Ramin 2019 Experimental observation of flow fields around active Janus spheres. Nat. Commun. 10 (1), 1–8.
  * Chen et al. [2021] Chen, Yibo, Chong, Kai Leong, Liu, Luoqin, Verzicco, Roberto & Lohse, Detlef 2021 Instabilities driven by diffusiophoretic flow on catalytic surfaces. J. Fluid Mech. 919.
  * Chen et al. [2019] Chen, Yufeng, Zhao, Huichan, Mao, Jie, Chirarattananon, Pakpong, Helbling, E Farrell, Hyun, Nak-seung Patrick, Clarke, David R & Wood, Robert J 2019 Controlled flight of a microrobot powered by soft artificial muscles. Nature 575 (7782), 324–329.
  * Datt et al. [2017] Datt, Charu, Natale, Giovanniantonio, Hatzikiriakos, Savvas G & Elfring, Gwynn J 2017 An active particle in a complex fluid. J. Fluid Mech. 823, 675–688.
  * Desai & Michelin [2021] Desai, Nikhil & Michelin, Sébastien 2021 Instability and self-propulsion of active droplets along a wall. Phys. Rev. Fluids 6 (11), 114103.
  * Dey et al. [2021] Dey, Ranabir, Buness, Carola M, Hokmabad, Babak Vajdi, Jin, Chenyu & Maass, Corinna C 2021 Oscillatory rheotaxis of active droplets in microchannels. arXiv preprint arXiv:2106.10222 .
  * Dey et al. [2022] Dey, Ranabir, Buness, Carola M, Hokmabad, Babak Vajdi, Jin, Chenyu & Maass, Corinna C 2022 Oscillatory rheotaxis of artificial swimmers in microchannels. Nat. Commun. 13 (1), 1–10.
  * Dreyfus et al. [2005] Dreyfus, Rémi, Baudry, Jean, Roper, Marcus L, Fermigier, Marc, Stone, Howard A & Bibette, Jérôme 2005 Microscopic artificial swimmers. Nature 437 (7060), 862–865.
  * Duan et al. [2015] Duan, Wentao, Wang, Wei, Das, Sambeeta, Yadav, Vinita, Mallouk, Thomas E & Sen, Ayusman 2015 Synthetic nano-and micromachines in analytical chemistry: sensing, migration, capture, delivery, and separation. Annu Rev Anal Chem 8, 311–333.
  * Dwivedi et al. [2021] Dwivedi, Prateek, Shrivastava, Atishay, Pillai, Dipin & Mangal, Rahul 2021 Rheotaxis of active droplets. Phys. Fluids 33 (8), 082108.
  * Ebrahimi et al. [2021] Ebrahimi, Nafiseh, Bi, Chenghao, Cappelleri, David J, Ciuti, Gastone, Conn, Andrew T, Faivre, Damien, Habibi, Neda, Hošovskỳ, Alexander, Iacovacci, Veronica, Khalil, Islam SM & others 2021 Magnetic actuation methods in bio/soft robotics. Adv. Funct. Mater. 31 (11), 2005137\.
  * Fu & Yan [2012] Fu, Jinglin & Yan, Hao 2012 Controlled drug release by a nanorobot. Nat. Biotechnol. 30 (5), 407–408.
  * Gao et al. [2010] Gao, Wei, Sattayasamitsathit, Sirilak, Manesh, Kalayil Manian, Weihs, Daniel & Wang, Joseph 2010 Magnetically powered flexible metal nanowire motors. J. Am. Chem. Soc. 132 (41), 14403–14405.
  * Ghosh & Fischer [2009] Ghosh, Ambarish & Fischer, Peer 2009 Controlled propulsion of artificial magnetic nanostructured propellers. Nano Lett. 9 (6), 2243–2245.
  * Golestanian et al. [2007] Golestanian, Ramin, Liverpool, TB & Ajdari, A 2007 Designing phoretic micro-and nano-swimmers. New J. Phys. 9 (5), 126.
  * Han et al. [2021] Han, E., Zhu, L., Shaevitz, J. W. & Stone, H. A. 2021 Low-Reynolds-number, biflagellated quincke swimmers with multiple forms of motion. Proc. Natl. Acad. Sci. USA 118 (29), e2022000118.
  * Hokmabad et al. [2022] Hokmabad, Babak Vajdi, Agudo-Canalejo, Jaime, Saha, Suropriya, Golestanian, Ramin & Maass, Corinna C 2022 Chemotactic self-caging in active emulsions. Proc. Natl. Acad. Sci. USA 119 (24), e2122269119.
  * Hokmabad et al. [2019] Hokmabad, Babak Vajdi, Baldwin, Kyle A, Krüger, Carsten, Bahr, Christian & Maass, Corinna C 2019 Topological stabilization and dynamics of self-propelling nematic shells. Phys. Rev. Lett. 123 (17), 178003.
  * Hokmabad et al. [2021] Hokmabad, Babak Vajdi, Dey, Ranabir, Jalaal, Maziyar, Mohanty, Devaditya, Almukambetova, Madina, Baldwin, Kyle A, Lohse, Detlef & Maass, Corinna C 2021 Emergence of bimodal motility in active droplets. Phys. Rev. X 11 (1), 011043\.
  * Hu et al. [2019] Hu, Wei-Fan, Lin, Te-Sheng, Rafai, Salima & Misbah, Chaouqi 2019 Chaotic swimming of phoretic particles. Phys. Rev. Lett. 123 (23), 238004.
  * Hu et al. [2022] Hu, Wei-Fan, Lin, Te-Sheng, Rafai, Salima & Misbah, Chaouqi 2022 Spontaneous locomotion of phoretic particles in three dimensions. Phys. Rev. Fluids 7 (3), 034003.
  * Iida et al. [2014] Iida, Keita, Kitahata, Hiroyuki & Nagayama, Masaharu 2014 Theoretical study on the translation and rotation of an elliptic camphor particle. Phys. D: Nonlinear Phenom. 272, 39–50.
  * Izri et al. [2014] Izri, Ziane, Van Der Linden, Marjolein N, Michelin, Sébastien & Dauchot, Olivier 2014 Self-propulsion of pure water droplets by spontaneous marangoni-stress-driven motion. Phys. Rev. Lett. 113 (24), 248302.
  * Jin et al. [2018] Jin, Chenyu, Hokmabad, Babak V, Baldwin, Kyle A & Maass, Corinna C 2018 Chemotactic droplet swimmers in complex geometries. J. Phys.: Condens. Matter 30 (5), 054003.
  * Jin et al. [2017] Jin, Chenyu, Krüger, Carsten & Maass, Corinna C 2017 Chemotaxis and autochemotaxis of self-propelling droplet swimmers. Proc. Natl. Acad. Sci. USA 114 (20), 5089–5094.
  * Joh & Fan [2021] Joh, Hyungmok & Fan, Donglei Emma 2021 Materials and schemes of multimodal reconfigurable micro/nanomachines and robots: Review and perspective. Adv. Mater. 33 (39), 2101965.
  * Kailasham & Khair [2022] Kailasham, R & Khair, Aditya S 2022 Dynamics of forced and unforced autophoretic particles. J. Fluid Mech. 948, A41.
  * Kasuo et al. [2019] Kasuo, Yui, Kitahata, Hiroyuki, Koyano, Yuki, Takinoue, Masahiro, Asakura, Kouichi & Banno, Taisuke 2019 Start of micrometer-sized oil droplet motion through generation of surfactants. Langmuir 35 (41), 13351–13355.
  * Kitahata et al. [2013] Kitahata, Hiroyuki, Iida, Keita & Nagayama, Masaharu 2013 Spontaneous motion of an elliptic camphor particle. Phys. Rev. E 87 (1), 010901.
  * Kohl et al. [2021] Kohl, Ryan, Corona, Eduardo, Cheruvu, Vani & Veerapaneni, Shravan 2021 Fast and accurate solvers for simulating janus particle suspensions in stokes flow. arXiv preprint arXiv:2104.14068 .
  * Koleoso et al. [2020] Koleoso, Mustaphis, Feng, Xue, Xue, Yuxiang, Li, Quan, Munshi, Tasnim & Chen, Xianfeng 2020 Micro/nanoscale magnetic robots for biomedical applications. Mater 8, 100085\.
  * Krüger et al. [2016] Krüger, Carsten, Klös, Gunnar, Bahr, Christian & Maass, Corinna C 2016 Curling liquid crystal microswimmers: A cascade of spontaneous symmetry breaking. Phys. Rev. Lett. 117 (4), 048003.
  * Lattuada & Hatton [2011] Lattuada, Marco & Hatton, T Alan 2011 Synthesis, properties and applications of Janus nanoparticles. Nano Today 6 (3), 286–308.
  * Lauga & Powers [2009] Lauga, Eric & Powers, Thomas R 2009 The hydrodynamics of swimming microorganisms. Rep. Prog. Phys. 72 (9), 096601.
  * Li [2022] Li, Gaojin 2022 Swimming dynamics of a self-propelled droplet. J. Fluid Mech. 934.
  * Li et al. [2022] Li, Meng, Pal, Aniket, Aghakhani, Amirreza, Pena-Francesch, Abdon & Sitti, Metin 2022 Soft actuators for real-world applications. Nat. Rev. Mater. 7 (3), 235–249.
  * Lin et al. [2020] Lin, Te-Sheng, Hu, Wei-Fan & Misbah, Chaouqi 2020 A direct poisson solver in spherical geometry with an application to diffusiophoretic problems. J. Comput. Phys. 409, 109362.
  * Lippera et al. [2020a] Lippera, Kevin, Benzaquen, Michael & Michelin, Sébastien 2020a Bouncing, chasing, or pausing: Asymmetric collisions of active droplets. Phys. Rev. Fluids 5 (3), 032201.
  * Lippera et al. [2020b] Lippera, Kevin, Morozov, Matvey, Benzaquen, Michael & Michelin, Sébastien 2020b Collisions and rebounds of chemically active droplets. J. Fluid Mech. 886.
  * Maass et al. [2016] Maass, Corinna C, Krüger, Carsten, Herminghaus, Stephan & Bahr, Christian 2016 Swimming droplets. Annu. Rev. Condens. Matter Phys. 7, 171–193.
  * Matsuda et al. [2016] Matsuda, Yui, Suematsu, Nobuhiko J, Kitahata, Hiroyuki, Ikura, Yumihiko S & Nakata, Satoshi 2016 Acceleration or deceleration of self-motion by the marangoni effect. Chem. Phys. Lett. 654, 92–96.
  * Meredith et al. [2020] Meredith, Caleb H, Moerman, Pepijn G, Groenewold, Jan, Chiu, Yu-Jen, Kegel, Willem K, van Blaaderen, Alfons & Zarzar, Lauren D 2020 Predator–prey interactions between droplets driven by non-reciprocal oil exchange. Nat. Chem. 12 (12), 1136–1142.
  * Michalet [2010] Michalet, Xavier 2010 Mean square displacement analysis of single-particle trajectories with localization error: Brownian motion in an isotropic medium. Phys. Rev. E 82 (4), 041914.
  * Michelin [2022] Michelin, Sébastien 2022 Self-propulsion of chemically active droplets. Annu. Rev. Fluid Mech. 55.
  * Michelin et al. [2013] Michelin, Sébastien, Lauga, Eric & Bartolo, Denis 2013 Spontaneous autophoretic motion of isotropic particles. Phys. Fluids 25 (6), 061701.
  * Moran & Posner [2017] Moran, Jeffrey L & Posner, Jonathan D 2017 Phoretic self-propulsion. Annu. Rev. Fluid Mech. 49, 511–540.
  * Morozov & Michelin [2019a] Morozov, Matvey & Michelin, Sébastien 2019a Nonlinear dynamics of a chemically-active drop: From steady to chaotic self-propulsion. J. Chem. Phys. 150 (4), 044110.
  * Morozov & Michelin [2019b] Morozov, Matvey & Michelin, Sébastien 2019b Self-propulsion near the onset of marangoni instability of deformable active droplets. J. Fluid Mech. 860, 711–738.
  * Nakata et al. [2006] Nakata, Satoshi, Kirisaka, Junko, Arima, Yoshie & Ishii, Toshio 2006 Self-motion of a camphanic acid disk on water with different types of surfactants. J. Phys. Chem. B 110 (42), 21131–21134.
  * Nasouri & Golestanian [2020a] Nasouri, Babak & Golestanian, Ramin 2020a Exact axisymmetric interaction of phoretically active Janus particles. J. Fluid Mech. 905.
  * Nasouri & Golestanian [2020b] Nasouri, Babak & Golestanian, Ramin 2020b Exact phoretic interaction of two chemically active particles. Phys. Rev. Lett. 124 (16), 168003.
  * Paredes et al. [2013] Paredes, Pedro, Hermanns, Miguel, Le Clainche, Soledad & Theofilis, Vassilis 2013 Order 104 speedup in global linear instability analysis using matrix formation. Comput Methods Appl Mech Eng 253, 287–304.
  * Paxton et al. [2004] Paxton, Walter F, Kistler, Kevin C, Olmeda, Christine C, Sen, Ayusman, St. Angelo, Sarah K, Cao, Yanyan, Mallouk, Thomas E, Lammert, Paul E & Crespi, Vincent H 2004 Catalytic nanomotors: autonomous movement of striped nanorods. J. Am. Chem. Soc. 126 (41), 13424–13431.
  * Paxton et al. [2006] Paxton, Walter F, Sundararajan, Shakuntala, Mallouk, Thomas E & Sen, Ayusman 2006 Chemical locomotion. Angew. Chem., Int. Ed. Engl. 45 (33), 5420–5429.
  * Peddireddy et al. [2012] Peddireddy, Karthik, Kumar, Pramoda, Thutupalli, Shashi, Herminghaus, Stephan & Bahr, Christian 2012 Solubilization of thermotropic liquid crystal compounds in aqueous surfactant solutions. Langmuir 28 (34), 12426–12431.
  * Peng & Schnitzer [2022] Peng, Gunnar G & Schnitzer, Ory 2022 Weakly nonlinear dynamics of a chemically active particle near the threshold for spontaneous motion. II. history-dependent motion. arXiv preprint arXiv:2211.09467 .
  * Picella & Michelin [2022] Picella, Francesco & Michelin, Sébastien 2022 Confined self-propulsion of an isotropic active colloid. J. Fluid Mech. 933.
  * Popescu et al. [2010] Popescu, Mihail Nicolae, Dietrich, S, Tasinkevych, M & Ralston, J 2010 Phoretic motion of spheroidal particles due to self-generated solute gradients. Eur. Phys. J. E 31 (4), 351–367.
  * Rao et al. [2015] Rao, K Jagajjanani, Li, Fei, Meng, Long, Zheng, Hairong, Cai, Feiyan & Wang, Wei 2015 A force to be reckoned with: a review of synthetic microswimmers powered by ultrasound. Small 11 (24), 2836–2846.
  * Schnitzer [2022] Schnitzer, Ory 2022 Weakly nonlinear dynamics of a chemically active particle near the threshold for spontaneous motion: adjoint method. Phys. Rev. Fluids .
  * Seemann et al. [2016] Seemann, Ralf, Fleury, Jean-Baptiste & Maass, Corinna C 2016 Self-propelled droplets. The European Physical Journal Special Topics 225 (11), 2227–2240.
  * Sharifi-Mood et al. [2016] Sharifi-Mood, Nima, Mozaffari, Ali & Córdova-Figueroa, Ubaldo M 2016 Pair interaction of catalytically active colloids: from assembly to escape. J. Fluid Mech. 798, 910–954.
  * Sitti et al. [2015] Sitti, Metin, Ceylan, Hakan, Hu, Wenqi, Giltinan, Joshua, Turan, Mehmet, Yim, Sehyuk & Diller, Eric 2015 Biomedical applications of untethered mobile milli/microrobots. P. IEEE 103 (2), 205–224.
  * Soto et al. [2021] Soto, Fernando, Karshalev, Emil, Zhang, Fangyu, Esteban Fernandez de Avila, Berta, Nourhani, Amir & Wang, Joseph 2021 Smart materials for microrobots. Chem. Rev. 122 (5), 5365–5403.
  * Suda et al. [2021] Suda, Saori, Suda, Tomoharu, Ohmura, Takuya & Ichikawa, Masatoshi 2021 Straight-to-curvilinear motion transition of a swimming droplet caused by the susceptibility to fluctuations. Phys. Rev. Lett. 127 (8), 088005.
  * Suematsu et al. [2010] Suematsu, Nobuhiko J, Ikura, Yumihiko, Nagayama, Masaharu, Kitahata, Hiroyuki, Kawagishi, Nao, Murakami, Mai & Nakata, Satoshi 2010 Mode-switching of the self-motion of a camphor boat depending on the diffusion distance of camphor molecules. J. Phys. Chem. C 114 (21), 9876–9882.
  * Suematsu et al. [2019] Suematsu, Nobuhiko J, Saikusa, Kazumi, Nagata, Toshiki & Izumi, Shunsuke 2019 Interfacial dynamics in the spontaneous motion of an aqueous droplet. Langmuir 35 (35), 11601–11607.
  * Suga et al. [2018] Suga, Mariko, Suda, Saori, Ichikawa, Masatoshi & Kimura, Yasuyuki 2018 Self-propelled motion switching in nematic liquid crystal droplets in aqueous surfactant solutions. Phys. Rev. E 97 (6), 062703.
  * Tang et al. [2020] Tang, Songsong, Zhang, Fangyu, Gong, Hua, Wei, Fanan, Zhuang, Jia, Karshalev, Emil, Esteban-Fernández de Ávila, Berta, Huang, Chuying, Zhou, Zhidong, Li, Zhengxing & others 2020 Enzyme-powered janus platelet cell robots for active and targeted drug delivery. Sci. Robot. 5 (43), eaba6137.
  * Thutupalli et al. [2018] Thutupalli, Shashi, Geyer, Delphine, Singh, Rajesh, Adhikari, Ronojoy & Stone, Howard A 2018 Flow-induced phase separation of active particles is controlled by boundary conditions. Proc. Natl. Acad. Sci. USA 115 (21), 5403–5408.
  * Thutupalli & Herminghaus [2013] Thutupalli, Shashi & Herminghaus, Stephan 2013 Tuning active emulsion dynamics via surfactants and topology. Eur. Phys. J. E 36 (8), 1–10.
  * Tomlinson [1862] Tomlinson, Charles, II. 1862 On the motions of camphor on the surface of water. Proc. Roy. Soc. London (11), 575–577.
  * Vyskocil et al. [2020] Vyskocil, Jan, Mayorga-Martinez, Carmen C, Jablonska, Eva, Novotny, Filip, Ruml, Tomas & Pumera, Martin 2020 Cancer cells microsurgery via asymmetric bent surface au/ag/ni microrobotic scalpels through a transversal rotating magnetic field. ACS Nano 14 (7), 8247–8256.
  * Wentworth et al. [2022] Wentworth, Ciera M, Castonguay, Alexander C, Moerman, Pepijn G, Meredith, Caleb H, Balaj, Rebecca V, Cheon, Seong Ik & Zarzar, Lauren D 2022 Chemically tuning attractive and repulsive interactions between solubilizing oil droplets. Angew. Chem., Int. Ed. Engl. 61 (32), e202204510.
  * Yang et al. [2023] Yang, Qianhong, Jiang, Maoqiang, Picano, Francesco, Jin, Chenyu & Lailai, Zhu 2023 Collective behaviors of isotropic phoretic microswimmers: from crystalline solids to active turbulence. arXiv preprint arXiv:2304.01593 .
  * Yariv & Kaynan [2017] Yariv, Ehud & Kaynan, Uri 2017 Phoretic drag reduction of chemically active homogeneous spheres under force fields and shear flows. Phys. Rev. Fluids 2 (1), 012201.