# Modeling of chemically active particles at an air-liquid interface

Shun Imamura1,2<EMAIL_ADDRESS>Toshihiro Kawakatsu1
1Department of Physics, Graduate School of Science, Tohoku University, Sendai
980-8578, Japan  2Mathematics for Advanced Materials-OIL, AIST-Tohoku
University, Sendai 980-8577, Japan

###### Abstract

The collective motion of chemically active particles at an air-liquid
interface is studied theoretically as a dynamic self-organization problem.
Based on a physical consideration, we propose a minimal model for self-
propelled particles by combining hydrodynamic interaction, capillary
interaction, driving force by Marangoni effect, and Marangoni flow. Our model
has successfully captured the features of chemically active particles, that
represent dynamic self-organized states such as crystalline, chain, liquid-
like and spreading states.

##  I Introduction

The collective behavior of chemically active particles in fluids is important
not only for the development of active matter physics[1, 2, 3, 4], which
focuses on the collective phenomena of self-propelled particles, but also for
a deeper understanding of the dynamic self-organization of systems that move
by chemical mechanisms such as bacteria[5, 6]. Chemically active particles
have been realized using colloids[7, 8, 9, 10], camphor particles[11, 12, 13,
14, 15], and droplets[16, 17, 18, 19].

Experimentally, the collective behavior of chemically active particles, such
as camphor particle system[20, 21, 22] and droplet system[23, 24, 25],
constrained on an air-liquid interface are particularly interesting as models
of two-dimensional wet active matter systems[2]. The wet system is defined as
those systems that are mainly dominated by hydrodynamic interaction. Here,
hydrodynamic interaction has a long-range nature and a lack of action-reaction
law, and can be an important factor in creating complex collective
behaviors[26]. For example, in an experiment by Soh et al., camphor particles
placed in a circular container show a dynamic self-organization where the
particles are arranged with a constant interval when the population of
particles is large[20]. Tanaka et al. observed complex dynamics in droplet
systems, where the behavior of the droplets changes with time, as if they
dance spontaneously[18, 25]. The physical mechanisms of these complex
collective dynamics are still unknown.

The motion of the chemically active particles at an air-liquid interface was
explained based on the following characteristic physical mechanisms[11, 20,
27, 28, 29, 30, 31, 32, 33, 34, 35, 36];

  * •

hydrodynamic interaction between particles (HI)[37, 38, 39],

  * •

lateral capillary force (capillary interaction) due to the deformation of the
air-liquid interface (CF)[33],

  * •

self-propelling force induced by the inhomogeneity in the surface tension due
to the gradient of the surrounding surfactant concentration field (SP)[11, 27,
15],

  * •

interaction between particles mediated by the concentration field of
surfactant (IC)[11],

  * •

Marangoni flow inside the liquid induced by the change in the surface tension
of the air-liquid interface due to surfactants (MF)[28, 31, 40, 29, 32, 35],

where the last 3 effects are originating from the concentration gradient of
the surfactant, and hereafter we will refer them as “Marangoni effects”.

In a colloidal particle system, Masoud and Shelley dealt with HI, IC, and MF
using direct numerical simulation[29], and Dominguez et al. dealt with HI, CF,
and MF analytically using the reaction diffusion equation[33]. For the camphor
particle system, Soh et al. dealt with HI, SP, IC, and MF numerically[20], and
Hirose et al. solved the diffusion equations for CF, SP, and IC numerically to
analyze the particle motion[36]. For self-propelled droplet systems, Yabunaka
and Yoshinaga studied the motion of two-particle systems by numerically and
analytically treating concentrations and fluid flow field, by incorporating
HI, SP, and IC[34]. Despite these successes in describing the collective
behaviors, the coupling between the concentration field and the fluid flow
field makes it difficult to identify the main factors of the physical
phenomena. It is also difficult to perform the direct numerical calculations
of many-particle systems. Therefore, it is important to incorporate each
physical element into a simple coarse-grained model in a way that the elements
can be turned on and off analytically.

In the present study, we propose a minimal model that can deal with
interacting many self-propelled particles while incorporating all elementary
physical mechanisms (hydrodynamic interaction, capillary interaction, and
Marangoni effect). The advantage of our model is that we need to solve neither
the concentration field nor the fluid flow fields, which reduces the
computational cost considerably.

The present article is organized as follows. First, we show the detail of our
proposed model in Sec. II. Then, the simulation method is explained in Sec.
III. In Sec. IV, we report the results of the simulations on single-particle,
two-particle and many-particle systems. Finally, in Sec. V, we draw our
conclusion.

##  II Model

Here, we construct a model of self-propelled particles by introducing the
hydrodynamic interaction, the capillary interaction, and the Marangoni
effects. For the Marangoni effect, we treat the self-propelling force due to
the gradient of the surfactant concentration field and the fluid flow due to
the difference of the interfacial tension of the air-liquid interface,
separately. Figure 1 shows schematic illustrations of individual physical
processes.

Figure 1:  Five physical mechanisms introduced in our model; long-range
hydrodynamic interaction (RPY type mobility tensor with mirror image method
applied), short-range capillary interaction, self-propelling force induced by
the Marangoni effect, short-range interaction due to the concentration field
of surfactant, and Marangoni flow that produces long-range repulsive force.

To describe the equations of motion for the particles, surfactant
concentration, and the liquid flow, we introduce a cartesian coordinate system
where the $z$-axis is set to the vertical direction to the air-liquid
interface and $x$ and $y$ axes in the interface (i.e. $z=0$). We also assume
that the air region and the liquid region correspond to $z>0$ and $z<0$,
respectively.

The equations of motion for the $i$-th particle are given by

$\displaystyle m\frac{d}{dt}\bm{v}_{i}^{p}(t)$
$\displaystyle=\bm{F}_{i}(t)+\bm{K}_{i}(t),$ (1)
$\displaystyle\frac{d}{dt}\bm{r}_{i}^{p}(t)$
$\displaystyle=\bm{v}_{i}^{p}(t),$ (2)

where $m$ is mass of a particle and $\bm{r}_{i}^{p}(t)$ and
$\bm{v}_{i}^{p}(t)$ represent the position and the velocity of $i$-th particle
at time $t$, respectively. $\bm{F}_{i}$ and $\bm{K}_{i}$ represent the driving
force and the viscous drag force, respectively. We model the Marangoni effects
based on an assumption that surfactant is emitted and diffuses from inside the
particle. To describe such a diffusion process, we assume a reaction-diffusion
equation for the surfactant emitted from an isolated particle located at the
origin of the system as follows[41, 42];

$\frac{\partial c}{\partial
t}+\bm{\nabla}\cdot\left(c\bm{v}^{f}\right)=D\bm{\nabla}^{2}c-\kappa\left(c-c_{\infty}\right)+A\delta(r),$
(3)

where, $c(\bm{r},t)$ represents the concentration field of the surfactant, and
$c_{\infty}$ is its value at the infinite distance from the particle in the
steady state. In Eq. (3), the first, second and third terms on the right-hand
side represent diffusion, consumption, and emission of the surfactant,
respectively, where $D,\kappa$ and $A$ are positive constants that represent
the diffusion coefficient, the consumption rate of surfactant that dissolves
into the bulk water from its surface per unit time, and the emission rate of
the surfactant. In this model of the surfactant diffusion, we approximate the
particle as a point particle for the sake of the analytical convenience[12].

By solving the reaction-diffusion equation, Eq. (3), in the steady state, we
can obtain the following concentration field,

$\displaystyle c(\bm{r})=\frac{A}{2\pi
D}\exp{\left[\frac{\bm{r}\cdot\bm{v}^{f}}{2D}\right]}K_{0}\left(2r\sqrt{\frac{\kappa}{4D}+\left(\frac{\bm{v}^{f}}{4D}\right)^{2}}\right),$
(4)

where $K_{0}(x)$ is the modified Bessel function of the second kind of order
0. For the derivation of Eq. (4), see Appendix A. The concentration field Eq.
(4) decays exponentially with the characteristic diffusion length
$\lambda=\sqrt{D/\kappa}$. We should note that, by assuming the source of the
surfactant as the delta function form, the solution of Eq. (3) corresponds to
the Green’s function, with which we can obtain the time evolution of the
surfactant concentration field emitted from a source with any shape, such as a
step function or a Gaussian function, by using the convolution integration.

We introduce the driving force $\bm{F}_{i}$ for $i$-th particle as

$\bm{F}_{i}(t)=\bm{F}_{i}^{\text{M}}(t)+\sum_{j(j\neq i)}\left(\bm{F}_{j\to
i}^{\text{C}}(t)+\bm{F}_{j\to i}^{\text{exc}}(t)\right),$ (5)

where $\bm{F}^{\text{M}}_{i}$ represents the force acting on the $i$-th
particle caused by the Marangoni effect, and $\bm{F}^{\text{C}}_{j\to i}(t)$
and $\bm{F}^{\text{exc}}_{j\to i}(t)$ are the lateral capillary force and the
excluded volume effect from $j$-th particle to $i$-th particle, respectively.
In the following, we will describe the detail of each force separately.

The Marangoni force $\bm{F}^{\text{M}}$ is divided into two parts as
$\bm{F}^{\text{M}}_{i}=\bm{F}^{\text{M-self}}_{i}+\bm{F}^{\text{M-dist}}_{i}$
(See Appendix B for details of calculation). $\bm{F}^{\text{M-self}}_{i}$ is
the self-propelling force given by

$\displaystyle\bm{F}^{\text{M-self}}_{i}=-\frac{\alpha
AR}{D}I_{1}\left(\frac{v^{f}_{i}R}{2D}\right)K_{0}\left(W_{i}R\right)\hat{\bm{v}}^{f}_{i},$
(6)

where $v^{f}_{i}$ is the fluid flow velocity at the position of $i$-th
particle, $R$ is the radius of the particle, $\alpha$ is the reduction rate
interfacial tension due to surfactant defined in Eq. (38),
$W_{i}=2\sqrt{\frac{\kappa}{4D}+\left(\frac{v^{f}_{i}}{4D}\right)^{2}}$, and
$I_{1}(x)$ is the modified Bessel function of the first kind of order 1. We
should note that, when the particle is moving with a velocity
$\bm{v}^{p}_{i}(t)$, the replacement of $\bm{v}^{f}_{i}$ by
$\bm{v}^{f}_{i}-\bm{v}^{p}_{i}$ should be made. Therefore, in Eq. (6), we
should regard $\bm{v}^{f}_{i}=\bm{v}^{f}(\bm{r}^{p}_{i}(t),t)$. On the other
hand, $\bm{F}^{\text{M-dist}}_{i}=\sum_{j\neq i}\bm{F}^{\text{M-dist}}_{j\to
i}$ is the interaction force due to the surfactant emitted by the neighboring
particles given by

$\displaystyle\bm{F}^{\text{M-dist}}_{j\to i}$ $\displaystyle=-\frac{\alpha
AR^{2}}{2D}\exp{\left[\frac{\bm{v}^{f}_{i}\cdot\bm{r}_{ij}}{2D}\right]}$
$\displaystyle\
\times\left[K_{0}(W_{j}r_{ij})\frac{\bm{v}^{f}_{i}}{2D}-W_{j}K_{1}(W_{j}r_{ij})\hat{\bm{r}}_{ij}\right].$
(7)

where $\bm{r}_{ij}(t)=\bm{r}^{p}_{i}(t)-\bm{r}^{p}_{j}(t)$,
$r_{ij}=|\bm{r}_{ij}|$ and $\hat{\bm{r}}_{ij}=\bm{r}_{ij}/r_{ij}$ is the unit
vector in the relative direction between $i$-th and $j$-th particles, and
$K_{1}(x)$ is the modified Bessel function of the second kind of order 1.

$\bm{F}^{\text{C}}_{j\to i}(t)$ in the second term of Eq. (5) represents the
driving force caused by the lateral capillary force[43] from $j$-th particle
to $i$-th particle given by

$\bm{F}^{\text{C}}_{j\to
i}(t)=-2\pi\gamma_{0}qQ_{i}Q_{j}K_{1}(qr_{ij})\hat{\bm{r}}_{ij}(t)\ \ \ (i\neq
j),$ (8)

where $\gamma_{0}$ represents the surface tension of the bare air-liquid
interface without surfactant, and $q^{-1}$ the capillary length and $Q_{i}$
the so-called “capillary charge” of the $i$-th particle, defined by
$Q_{i}=R^{c}_{i}\sin{\psi_{i}}$. $R^{c}$ represents the radius of contact line
between three phases, i.e. air-liquid-particle, and $\psi$ represents the
constant angle at the meniscus between the particle and the air-liquid
interface.

Finally, $\bm{F}^{\text{exc}}_{j\to i}(t)$ in the second term of Eq. (5)
represents the excluded volume effect given by the WCA potential[44] between
contacting particles.

We assume an overdamped dynamics where the driving force $\bm{F}_{i}(t)$ in
Eq. (5) is balanced by the viscous drag force $\bm{K}_{i}$ defined by

$\displaystyle\bm{K}_{i}(t)$
$\displaystyle=-\zeta\left(\bm{v}^{p}_{i}(t)-\bm{v}^{f}_{i}(t)\right),$ (9)
$\displaystyle\bm{v}^{f}_{i}(t)$ $\displaystyle=\bm{u}(\bm{r}^{p}_{i}(t))+\int
d\bm{r}^{\prime}\bm{T}(\bm{r}^{p}_{i}-\bm{r}^{\prime})\cdot\bm{F}^{f}(\bm{r}^{\prime},t),$
(10) $\displaystyle\bm{F}^{f}(\bm{r},t)$
$\displaystyle=\sum_{i}\zeta\left(\bm{v}^{p}_{i}(t)-\bm{v}^{f}_{i}(t)\right)\delta(\bm{r}^{p}_{i}(t)-\bm{r}),$
(11)

where $\zeta$ is the friction coefficient of the particle floating on the
liquid-air interface. We assume that the fluid velocity at the $i-$th particle
position $\bm{v}^{f}_{i}(t)$ is described by the Stokes equation
$\mu\bm{\nabla}^{2}\bm{v}^{f}(\bm{r},t)=\bm{\nabla}P(\bm{r},t)-\bm{F}^{f}(\bm{r},t)$,
where the inertia term is neglected and the incompressibility condition
$\bm{\nabla}\cdot\bm{v}^{f}(\bm{r},t)=0$ is imposed. Here, $\mu$,
$P(\bm{r},t)$ and $\bm{F}^{f}(\bm{r},t)$ represent the viscosity of the fluid,
the pressure field and the external force field, respectively.
$\bm{u}(\bm{r},t)$ denotes the Marangoni flow. (See Appendix B.3 for details
of calculation.) $\bm{T}(\bm{r})$ is the Green’s function for hydrodynamic
interaction under the boundary condition that the fluid velocity in the
$z$-direction vanishes at the air-liquid interface ($z=0$) and is given by
$\bm{T}(\bm{r}-\bm{r}^{\prime})=\bm{G}(\bm{r}-\bm{r}^{\prime})+\bm{G}(\bm{r}-\bm{r}^{\prime}+2h\bm{e}_{z})\cdot\bm{P}_{z},$
where $h$ is the $z$-coordinate of the center of mass of the particle,
$\bm{P}_{z}=\bm{1}-2\bm{e}_{z}\bm{e}_{z}$ represents the reflection operator
by a mirror image[39] and $\bm{e}_{z}$ is the unit vector in the
$z$-direction. $\bm{G}(\bm{r})$ expresses the RPY type mobility tensor[37, 38]

$\displaystyle\bm{G}(\bm{r})=\frac{1}{8\pi\mu
r}\left[(\bm{1}+\hat{\bm{r}}\hat{\bm{r}})+\frac{1}{3}\left(\frac{R}{r}\right)^{2}(\bm{1}-3\hat{\bm{r}}\hat{\bm{r}})\right],$
(12)

where the first term on the right-hand side is a stokeslet, and the second
term is a source doublet. The boundary condition at the air-fluid interface is
specified by the balance between the forces acting on the interface as
$\left.\frac{\partial v^{f}_{x}}{\partial z}\right|_{z=0}=\left.\frac{\partial
v^{f}_{y}}{\partial z}\right|_{z=0}=0,\ v^{f}_{z}|_{z=0}=0$.

##  III Simulation method

In the present study, we consider the case that the center of mass of each
particle is located on the air-liquid interface, i.e., $h\sim 0$, where the
effective velocity field at $xy$-plane is described by $\bm{T}(\bm{r})\sim
2\bm{G}(\bm{r})$. Due to such a simplification, the parameter $\zeta$, which
is given by $\zeta\sim 3\pi\mu R$, is the friction constant of the particle at
the air-liquid interface[45]. Although our model is a minimal model on the
coarse-grained scale, it is still too complicate to be solved analytically.
Thus, we introduce a further simplification that $R_{c}$ and $\psi$ are common
to all the particles, i.e., $R^{c}_{i}=R$ and $\psi_{i}=\psi$. Therefore,
$Q_{i}=Q=R^{c}\sin{\psi}\sim R\sin{\psi}$. In addition, we expand the velocity
field $\bm{v}^{f}_{i}(t)$ in Eqs. (10) and (11) up to second order in
$(\zeta\bm{T})$ to obtain analytically tractable model equations.

Let us consider the correspondence between our simulation and the existing
experiments[46, 15]. First, we rewrite the reaction-diffusion equation Eq. (3)
and the equation of motion Eq. (1) in non-dimensional forms by using the units
of length, time and energy, $L_{0},T_{0}$ and $E_{0}$, where these unit
quantities are estimated for the experimental situation[46, 15] as
$L_{0}\equiv R\sim 10^{-3}$ [m], $T_{0}\equiv\tau=m/\zeta$ [s], and
$E_{0}\equiv 2\pi\gamma_{0}R^{2}\sim 10^{-7}$ [J], respectively. Here, $\tau$
is the relaxation time for the particle motion and can be estimated as
$\tau\sim 1$ [s] for $m\sim 10^{-5}$ [kg].

The dimensionless parameters for the capillary interaction used in the
simulations are $\psi\sim 0.1$ [rad] and $\text{Bo}\equiv(qR)^{2}\sim 0.1$
(determines the capillary length $q^{-1}$), the latter being related to the
lateral capillary force [See Eq. (8)]. Here, Bo represents the Bond number
which is the ratio of the buoyancy force to the surface tension between the
particle and the liquid.

On the other hand, the dimensionless quantities in reaction-diffusion equation
for the surfactant [See Eq. (3)] are defined as follows;

$\displaystyle A^{*}=AT_{0},\ \kappa^{*}=\kappa T_{0},\
D^{*}=D\frac{T_{0}}{L_{0}^{2}},$ (13)

where $D^{*}$ corresponds to the dimensionless diffusion length
$\lambda=\sqrt{D/\kappa}$. Furthermore, we define the following dimensionless
quantities for the equation of motion for the particles [See Eq. (1)]

$\displaystyle\text{Pe}=\frac{vR}{D},\ \text{Ma}=\frac{\alpha AR}{\mu D^{2}},$
(14)

where Pe is the Péclet number and is defined as the nondimensional self-
propelling speed of the particles[47, 15], and Ma is the solute Marangoni
number, a dimensionless quantity defined by the ratio of the driving force of
the concentration gradient to the viscous friction force. Furthermore, in the
direct particle-particle interaction potential, we define
$\gamma_{0}^{*}=E_{0}T_{0}\sin^{2}{\psi}/(L_{0}^{2}\zeta)$, where
$\gamma_{0}^{*}$ corresponds to the ratio of the lateral capillary force to
the viscous friction force.

##  IV Results and Discussion

###  IV.1 Single-particle system

First, we discuss the behavior of single-particle systems. Figure 2 denotes
(a) the color display of the concentration field of surfactant around a
particle, (b) its cross-sections along $x$\- and $y$-axes, and (c) swimming
speed (Pe) as a function of Marangoni number. The data shown in (a) and (b)
are obtained using Eq. (4), and the curve shown in (c) is calculated using Eq.
(1). Experimentally, the Marangoni number is controlled by changing the
viscosity, for example by adding glycerin to water[12, 15].

Figure 2:  Self-driven motion for the single-particle system. (a) Color
display of the concentration field of surfactant around the particle, where
the particle is moving in the positive direction of the $y$-axis, (b) its
cross-sections along $x$\- and $y$-axes, and (c) swimming speed (Pe) as a
function of Marangoni number. The parameters are set as $D^{*}=0.15$ and
$\kappa^{*}=0.01$. The red circle represents the particle in (a), and the red
region represents the inside of the particle in (b). The data shown in (a) and
(b) and the dot in (c) correspond to the case with $\bm{v}^{f}=(0.0,-0.1)$.

The non-dimensional equation of motion for a single particle is represented as

$\displaystyle\frac{d}{dt}v=-v+F^{\text{M-self}}(v),$ (15)

where $v$ is the swimming speed of the particle. Let us consider the range of
stability of the steady state solution $v_{0}$ of Eq. (15). To show such a
condition, we expand $F^{\text{M-self}}(v)$ around $v=0$ as

$\displaystyle F^{\text{M-self}}(v)=C^{(1)}v+C^{(3)}v^{3}+\cdots,$ (16)

where

$\displaystyle C^{(n)}=\frac{1}{n!}\left.\frac{\partial^{n}}{\partial
v^{n}}F^{\text{M-self}}(v)\right|_{v=0}\ \ \ n=1,2,...\ .$ (17)

Here, $F^{\text{M-self}}(0)=0$ and $C^{(n)}=0$ for $n=2,4,6,...$, because
$F_{\text{M-self}}(v)$ is an odd function of $v$ due to the symmetry of the
system. Therefore, we obtain

$\displaystyle\frac{dv}{dt}$
$\displaystyle=\left(-1+C^{(1)}\right)v+C^{(3)}v^{3}+\cdots,$ (18)

with which we can clarify the condition for stably self-propelling motion as

$\displaystyle-1+C^{(1)}$ $\displaystyle>0,$ (19) $\displaystyle C^{(3)}$
$\displaystyle<0.$ (20)

Here,

$\displaystyle C^{(1)}$
$\displaystyle=\frac{\text{Ma}}{12\pi}K_{0}\left(\sqrt{\frac{\kappa^{*}}{D^{*}}}\right),$
(21) $\displaystyle C^{(3)}$
$\displaystyle=\frac{\text{Ma}}{384\pi(D^{*})^{2}}\left[K_{0}\left(\sqrt{\frac{\kappa^{*}}{D^{*}}}\right)-4\sqrt{\frac{D^{*}}{\kappa^{*}}}K_{1}\left(\sqrt{\frac{\kappa^{*}}{D^{*}}}\right)\right].$
(22)

We can confirm that the condition Eq. (20) is always satisfied because the
parameters $D^{*}$ and $\kappa^{*}$ are positive and $K_{0}(x)-4K_{1}(x)/x<0$
for $x>0$. When this condition is satisfied, the steady-state velocity
obtained with the expansion shown in Eq. (18) up to the third order of the
velocity $v$ is given by

$\displaystyle v_{0}=\sqrt{\frac{C^{(1)}-1}{-C^{(3)}}}.$ (23)

This result means that there is a threshold value for Ma where a single
particle starts to move spontaneously (See Fig. 2(c)). Similarly to this
result, the experimental data shows a positive correlation between Pe and Ma,
and the self-propelling speed has a threshold value for Ma[12]. Thus, our
model can reproduce the onset of the self-propelling motion induced by a
spontaneous symmetry breaking of the concentration field.

The self-driven velocity of each particle is given by Eq. (23), which means
that the threshold is determined by $C^{(1)}=1$. Substituting Eqs. (19) and
(21) into this condition, the expression of the threshold value of Ma in the
self-driven motion is obtained. In this expression of the threshold value,
$D^{*}$ and $\kappa^{*}$ play opposite roles, i.e. larger $D^{*}$ decreases
the threshold value of Ma while larger $\kappa^{*}$ increases the threshold
value of Ma. These tendencies can be understood as follows. When $D^{*}$ is
large, the concentration distribution of surfactant spreads over a wide area
in an asymmetric manner, resulting in an promoted mobility. On the other hand,
when $\kappa^{*}$ is large, the concentration distribution of surfactant
decays faster before it spreads over a distance, resulting in a decay of the
concentration gradient and a decrease in mobility. In our simulations, we
chose the values of the parameters $D^{*}$ and $\kappa^{*}$ based on the
validation shown in Appendix C.

In Fig. 2(c), the self-propelling motion of each particle occurs when Ma is
larger than a threshold value. The definition of Ma in Eq. (14) means that a
large Ma corresponds to a small viscosity $\mu$, a small diffusion constant
$D$, a large emission rate of surfactant $A$, or a large reduction rate of the
interfacial tension due to surfactant $\alpha$. All these properties tend to
enhance the asymmetric distribution of the surfactant around the particle,
which leads to a large self-driving force. When the self-propelling force
exceeds the other viscous drag force induced by the Marangoni flow, the
condition Eq. (19) is satisfied and the particle starts to move spontaneously.

###  IV.2 Two-particle system

Before discussing the dynamical behavior of many particle systems, we
investigate the elementary components of the interaction between two particles
as functions of their separation distance. Figure 3 shows the individual
components of the force acting between two particles. As was shown in Fig.
2(c), the self-driven motion does not occur in Fig. 3(a) ($\text{Ma}=1.0$) but
occurs in Fig. 3(b) ($\text{Ma}=50.0$). Figure 3 shows that there is a
characteristic inter-particle distance where the attraction and repulsion are
switched.

Figure 3:  Forces acting between two particles for
$D^{*}=0.15,\kappa^{*}=0.01,\gamma_{0}^{*}=0.1$, and (a) Ma $=1.0$ and (b) Ma
$=50.0$. “Marangoni flow” and “Capillary” denote the repulsive viscous force
due to the Marangoni flow and the attractive force caused by capillary
interaction, respectively.

In Fig. 3 (a), the self-propelling force is small due to the small Ma as was
shown in Fig. 2(c). For this reason, the main force acting between the two
particles is the capillary force caused by the depression of the water
surface. As the interparticle distance increases, the capillary force decays
rapidly, while the Marangoni flow decays slowly, leading to the dominance of
the repulsive interaction induced by the Marangoni flow. A competition between
these two tendencies defines a threshold distance, inside and outside of which
the interparticle interaction changes its nature from attractive to repulsive
ones.

In the case of large Ma Fig. 3(b), both self-propelling force and Marangoni
flow are large (As Ma is proportional to $\alpha$, a large Ma means that the
change in the interfacial tension caused by the surfactant is large.). In this
case, the behavior at short distances is dominated by the interplay between
the self-propelling force and the viscous drag force due to the Marangoni
flow. In the initial rest state, the direction of the initial particle motion
is determined by the self-propelling force induced by the high surfactant
concentration in the interparticle region, which reduces the surface tension
in that region. This inhomogeneity in the surface tension generates an outward
force due to the higher surface tension in the outside region, leading to an
initial repulsive force. Then, the Marangoni flow enhances this repulsive
interaction. On the other hand, at large interparticle distance, the capillary
interaction and the self-propelling force $F^{\text{M-self}}$ determine the
direction of the particle motion. In the initial rest state, even though the
capillary force is weak, it induces an attractive interaction, which leads to
the particle motion toward the attractive direction. This initial motion
induces the self-driven motion in the attractive direction because the self-
propelling force is insensitive to the separation between the particles.
Therefore, when the distance between the two particles is large, the self-
propelling force acts as the attractive force.

Figure 4 shows phase diagrams and trajectories of individual cases for the
two-particle system. These phase diagrams are made under the initial
conditions (a) $r_{12}(0)=3.0$ (inside of the characteristic inter-particle
distance) and (b) $r_{12}(0)=8.0$ (outside of the characteristic distance)
(See Fig. 3). The states of two-particle system are classified into three
motions; contacted, repulsive, and non-contacted motions, respectively. The
trajectories of each motion are shown in Fig. 4(c). Additionally, the boundary
line between the regions of repulsive and non-contacted motions is determined
by the threshold in the motion of the single particle system shown for Fig.
2(c).

Figure 4:  Phase diagrams and the typical trajectories of the two-particle
systems for the relative motion of two particles for $D^{*}=0.15$ and
$\kappa^{*}=0.01$. Phase diagrams are obtained for (a) $r_{12}(0)=3.0$ and (b)
$r_{12}(0)=8.0$, where $r_{12}(0)$ denotes the inter-particle distance at
$t=0$ and the phase boundaries are the guide to the eyes. Figure (c) shows
trajectories of the two particles starting from $r_{12}(0)=8.0$ for various
values of Ma and $\gamma_{0}^{*}$.

###  IV.3 Many-particle system

Let us discuss the behavior of collective motion of the particles by showing
simulation results, such as phase diagrams and order parameters. In Fig. 5, we
show 4 typical collective behaviors; i.e. crystalline, chain, spreading and
liquid-like states. (See the Supplemental Material [48] for movies of these
behaviors.) These collective behaviors correspond to experimentally observed
states. For example, the crystalline (Fig. 5(a)) and the chain states (Fig.
5(b)) reproduce the collective dynamics of self-propelled droplets[25], and
the spreading state (Fig. 5(c)) corresponds to the dynamic self-organization
of the camphor disks[20].

Figure 5:  Four typical collective motions; (a) crystalline, (b) small chains
and clusters, (c) spreading and (d) liquid-like states, respectively. The
directions of velocities of the particles are shown by arrows and colors
defined in the color legend in (a).

In order to identify the phase boundaries of each state shown in Fig. 5, we
introduce several order parameters to characterize the collective behaviors.
The crystalline state can be characterized by the 6-fold bond-orientational
order parameter[49] defined by

$\displaystyle\psi_{6}$
$\displaystyle=\frac{1}{N}\sum_{j=1}^{N}|\psi_{6}^{j}|,$ (24)
$\displaystyle\psi_{6}^{j}$
$\displaystyle=\frac{1}{Z_{j}}\sum_{k}^{Z_{j}}e^{i6\theta_{jk}},$ (25)

where $Z_{j}$ is the coordination number of $j$-th particle obtained from a
Voronoi construction for the particle configuration, and $\theta_{jk}$ is the
angle between a reference axis and the direction of the bond between $j$-th
particle and its $k$-th neighbor. $\psi_{6}=1$ means perfect hexagonal
ordering, whereas completely disordered structures give $\psi_{6}=0$.

As the above bond-orientational order parameter $\psi_{6}$ gives a large value
not only for crystalline structure but also a straight string-like structures,
$\psi_{6}$ cannot distinguish chain state from crystalline state. Thus, we
introduce an additional orientational order parameter defined by

$\displaystyle\phi$ $\displaystyle=\frac{1}{M}\sum_{i}\phi^{i},$ (26)
$\displaystyle\phi^{i}$
$\displaystyle=\frac{1}{{}_{{}^{\mathcal{N}_{i}}}C_{2}}\sum_{(j,k)\in\mathcal{S}_{i},j\neq
k}\left[\frac{2}{3}\left(\frac{1}{2}-\frac{\bm{r}_{ij}\cdot\bm{r}_{jk}}{r_{ij}r_{jk}}\right)\right],$
(27)

where $M$ represents the total number of particles that are in contact with
two or more other particles, $\mathcal{N}_{i}$ is the number of particles in
contact with the $i$-th particle, ${}_{{}^{\mathcal{N}_{i}}}C_{2}$ is the
binomial coefficient, and $\mathcal{S}_{i}$ denotes the search region for the
other particles that contact the $i$-th particle, i.e. a circle with a radius
of the order of the particle diameter centered at the center of mass of $i$-th
particle. The order parameter $\phi$ is defined for clusters composed of three
or more particles, and this parameter takes the value $\phi=1$ in the case of
the chain state, and $\phi=0$ in the case of the crystalline state. Here, the
order parameter $\phi^{i}$ is defined for $i$-th particle, where $\phi^{i}=1$
in the case that the relative positions of neighboring particles of the $i$-th
particle align in a straight line.

Figure 6 shows the phase diagram of parameter regions for each state
determined using the order parameters introduced above. Here, we take the
Marangoni number Ma and the magnitude of the capillary force $\gamma_{0}^{*}$
as independent parameters. In Fig. 6(a), the capillary interaction is dominant
in the crystalline state. On the other hand, when the capillary force is
small, the states are classified according to the Marangoni number into
spreading and liquid-like states. The boundary between the spreading state and
the liquid-like state locates at the threshold value of Ma for the self-driven
motion found in Fig. 2(c). This means that the spreading state is dominated by
the Marangoni flow only, and the liquid-like states are driven by the self-
driven motion. In Fig. 6(b), the chain state is realized in a limited domain
of Ma $\sim 100.0$ and $\gamma_{0}^{*}\sim 5.0$. This means that the chain
state is caused by a balance of all physical factors introduced in our model;
i.e. hydrodynamic interaction, capillary interaction, and Marangoni effect.

We show the phase diagram for many-particle system in Fig. 6(c). The
boundaries of this diagram are guide to the eyes that are determined based on
the results of the order parameter calculations shown in Figs. 6 (a) and (b).
The region of contacted, repulsive, and non-contacted motions for two-particle
system in Fig. 4 correspond to the crystalline, diffusion, and liquid states
in Fig. 6, respectively. The region of the chain state is a newly added phase
in the many-particle system.

In the region of large $\gamma_{0}^{*}$, the attractive capillary force due to
the depression of the water surface causes the particles to coagulate into a
crystalline state. When $\gamma_{0}^{*}$ is small, the attractive capillary
fore is small, and the Marangoni effect is dominant. When Ma is large, the
self-propelling force becomes large, and the particles start to move freely
and realizes the liquid state. In this case, the Marangoni flow acts as a
repulsive interparticle force, which prevents the particles to coalesce at
short distance. On the other hand, when Ma and $\gamma_{0}^{*}$ are small,
particles are repelling with each other to form the spreading state. The chain
state occurs when both Ma and $\gamma_{0}^{*}$ are large. Since
$\gamma_{0}^{*}$ is large, there is a large attraction between the particles
due to the depression of the water surface. This causes the particles to form
clusters. Since Ma is also large, however, the particles try to keep moving.
Such a competition results in a chain state where the two mechanisms are well-
balanced, i.e. the particles try to move spontaneously due to the Marangoni
effect while the capillary interaction keeps the clusters.

Figure 6:  Stable regions of individual states on the Ma-$\gamma_{0}^{*}$
plane for $N=30$ systems. The stable regions are determined with use of the
order parameters, (a) bond-orientational order parameter $\psi_{6}$ and (b)
orientational order parameter $\phi$. In this analysis, we averaged over 50
runs starting from different initial configurations. (c) From the order
parameters, we can draw phase diagram for many-particle systems, where the
phase boundaries are guide to the eyes.

##  V Conclusion

We have developed a minimal model to describe chemically active particles at
an air-liquid interface. In the modeling, the interactions between particles
are decomposed into hydrodynamic interactions, capillary interactions, driving
forces due to Marangoni effect, and Marangoni flow.

Analysis on the equation of motion for the single-particle system showed the
existence of a threshold for the self-driven motion when the viscosity of
liquid is changed. (i.e. Ma is increased.) This finding was confirmed by
numerical calculations, where the model parameters are chosen so that the
model corresponds to the experimental situation[12, 15]. These results show
that our simple model can capture the essential properties of the experimental
systems.

For the two-particle systems, we discussed the behavior of the inter-particle
interaction as a function of the separation distance, which shows a switching
from repulsive to attractive at a certain threshold distance. Such a
separation distance dependence of the inter-particle interaction is in good
agreement with the previous simulations reported by Soh et al.[20] for the
case without the capillary interaction. Based on these results, we can
construct the phase diagrams for these two-particle systems, with which we can
understand the elementary behavior of many-particle system.

Decomposing the inter-particle interaction into contributions from different
physical elements listed in Figure 1, we find that the concentration field and
the repulsion by Marangoni flow give the most dominant contribution.

Simulations on many-particle systems reproduced the collective behaviors found
in the existing experiments, from which we can identify the major contribution
for each state, i.e. the spreading state is caused by Marangoni flow as a
major physical factor, the crystalline state by capillary interaction, the
liquid-like state by self-propelling force, and the chain state by all
physical elements included in our model. Our model, which includes the effect
of the hydrodynamic interactions, can reproduce both crystalline and chain
states[25] and the self-organization of particles with regular intervals[20].
Furthermore we could identify the stable region for each state expressed by
two parameters, i.e. Marangoni number Ma and the amplitude of capillary
interaction $\gamma_{0}^{*}$.

With a further extension and an improvement of our model by introducing time-
dependent model parameters, we will try to explain time-dependent non-steady
collective behaviors. By using our model, one can explore the behavior of
many-particle system, with which we hope to observe typical dynamic phenomena
such as motility induced phase separation in chemically active particles.
Applying our model to phenomena with much larger length and time scales, such
as scaling behaviors, will be another interesting extension of the current
study, where a further coarse-graining and multiscale treatment would be
necessary.

###### Acknowledgements.

The present work is partially supported by the Grant-in-Aid for Scientific
Research from The Ministry of Education, Culture, Sports, Science and
Technology of Japan (Grant No. 19H01858).

##  Appendix A The derivation of the steady state concentration field of
surfactant

We consider the following reaction-diffusion equation for surfactant density
$\Gamma(\bm{r},t)$

$\displaystyle\frac{\partial\Gamma}{\partial
t}+\bm{\nabla}\cdot\left(\Gamma\bm{v}_{f}\right)=D\bm{\nabla}^{2}\Gamma-\kappa\left(\Gamma-
c_{\infty}\right)+A\delta(\bm{r}),$ (28)

where $\delta(\bm{r})$ is Dirac’s delta function. We solve this equation under
the quasi-static approximation where the velocity field $\bm{v}_{f}$ is
assumed to be in a steady state. Then, the concentration field of the
surfactant at time $t$ generated at a previous time $t^{\prime}$ can be
expressed in terms of the Green’s function. This Green’s function
$\Gamma_{t^{\prime}}(\bm{r},t-t^{\prime})$ is obtained by solving

$\displaystyle\frac{\partial}{\partial
t}\Gamma_{t^{\prime}}+\bm{v}_{f}\cdot\bm{\nabla}\Gamma_{t^{\prime}}=D\bm{\nabla}^{2}\Gamma_{t^{\prime}}-\kappa\Gamma_{t^{\prime}}+A\delta(\bm{r})\delta(t-t^{\prime}),$
(29)

where the concentration field for steady state of Eq. (28) is given by

$\displaystyle\Gamma(\bm{r})=\int_{-\infty}^{t}dt^{\prime}\Gamma_{t^{\prime}}(\bm{r},t-t^{\prime}).$
(30)

We define the Fourier transform for position $\bm{r}$ as

$\displaystyle\widetilde{G}_{k}(\bm{k},t)$
$\displaystyle=\frac{1}{2\pi}\int_{-\infty}^{\infty}\int_{-\infty}^{\infty}G(\bm{r},t)e^{-i\bm{k}\cdot\bm{r}}dxdy,$
(31) $\displaystyle G(\bm{r},t)$
$\displaystyle=\frac{1}{2\pi}\int_{-\infty}^{\infty}\int_{-\infty}^{\infty}\widetilde{G}_{k}(\bm{k},t)e^{i\bm{k}\cdot\bm{r}}dk_{x}dk_{y},$
(32)

and for time $t$ as

$\displaystyle\widetilde{G}_{\omega}(\bm{r},\omega)$
$\displaystyle=\frac{1}{\sqrt{2\pi}}\int_{-\infty}^{\infty}G(\bm{r},t)e^{-i\omega
t}dt,$ (33) $\displaystyle G(\bm{r},t)$
$\displaystyle=\frac{1}{\sqrt{2\pi}}\int_{-\infty}^{\infty}\widetilde{G}_{\omega}(\bm{r},\omega)e^{i\omega
t}dt.$ (34)

Using these Fourier transforms, we can obtain the Green’s function as

$\displaystyle\Gamma_{t^{\prime}}(\bm{r},s)$ $\displaystyle=\frac{A}{4\pi
Ds}\exp{\left[-\left(\kappa
s+\frac{\left(\bm{r}-\bm{v}_{f}s\right)^{2}}{4Ds}\right)\right]}.$ (35)

Therefore, the following expression for the Green’s function is obtained from
Eqs. (30) and (35);

$\displaystyle\Gamma(\bm{r})$ $\displaystyle=\frac{A}{2\pi
D}\exp{\left[\frac{\bm{r}\cdot\bm{v}_{f}}{2D}\right]}K_{0}\left(2r\sqrt{\frac{\kappa}{4D}+\left(\frac{\bm{v}_{f}}{4D}\right)^{2}}\right),$
(36)

where $K_{0}(x)$ is the $0$-th order modified Bessel function of the second
kind. Using the above Green’s function, we can express the steady-state
concentration field for the source with any shape $f(\bm{r})$ as

$\displaystyle c(\bm{r})$ $\displaystyle=\int
d\bm{r}^{\prime}f(\bm{r}^{\prime})\Gamma(\bm{r}-\bm{r}^{\prime}).$ (37)

##  Appendix B Derivation of the driving force caused by Marangoni effect

The driving force due to Marangoni effect $\bm{F}^{\text{M}}_{i}$ is divided
into two parts as
$\bm{F}^{\text{M}}_{i}=\bm{F}^{\text{M-self}}_{i}+\bm{F}^{\text{M-dist}}_{i}$,
where $\bm{F}^{\text{M-self}}_{i}$ is the self-propelling force and
$\bm{F}^{\text{M-dist}}_{i}$ is the interaction force due to the surfactant
distribution emitted by the neighboring particles. In addition, we introduce
driving force by Marangoni flow caused by the change of the surface tension of
air-liquid interface[40].

###  B.1 Self-propelling force

Let us consider self-propelling force caused by the Marangoni effects. The
local surface tension of the air-liquid interface is assumed to be linearly
dependent on the surfactant concentration as

$\gamma(\bm{r})=\gamma_{0}-\alpha c(\bm{r}),$ (38)

where this surface tension is defined at the interface between the liquid and
the air. By integrating the force due to this surface tension along the three-
phase contact line (hereafter denoted as $C$) around the particle, we can
obtain the self-propelling force as

$\bm{F}^{\text{M-self}}=\oint_{C}(\gamma_{0}-\alpha
c(\bm{r}))\widehat{\bm{r}}\ dl=-\alpha\oint_{C}c(\bm{r})\widehat{\bm{r}}\ dl,$
(39)

where $\hat{\bm{r}}$ denotes unit normal vector from the center of mass of the
particle to the point on the contact line. This integration leads to the
expression of the driving force as

$\displaystyle\bm{F}^{\text{M-self}}_{i}=-\frac{\alpha
AR}{D}I_{1}\left(\frac{v^{f}_{i}R}{2D}\right)K_{0}\left(W_{i}R\right)\hat{\bm{v}}^{f}_{i},$
(40)

where $W_{i}=2\sqrt{\frac{\kappa}{4D}+\left(\frac{v^{f}_{i}}{4D}\right)^{2}}$,
$R$ indicates the radius of the particle and $I_{1}(x)$ is the 1-st order
modified Bessel function of the first kind. Note that, when the $i$-th
particle is moving with a velocity $\bm{v}^{p}_{i}$, we should replace
$\bm{v}^{f}_{i}$ by $\bm{v}^{f}_{i}-\bm{v}^{p}_{i}$, where
$\bm{v}^{f}_{i}=\bm{v}^{f}(\bm{r}^{p}_{i}(t),t)$.

###  B.2 The influence of the surrounding particles

Let us consider the inter-particle force caused by the concentration field of
the surfactant emitted by the other surrounding particles. By using the steady
state solution of the concentration field, we can easily describe the inter-
particle forces as

$\displaystyle\bm{F}^{\text{M-dist}}_{i}$ $\displaystyle=\sum_{j\neq
i}\bm{F}^{\text{M-dist}}_{j\to i},$ (41)
$\displaystyle\bm{F}^{\text{M-dist}}_{j\to i}$
$\displaystyle=-\alpha\oint_{C_{i}}c(\bm{r}+\bm{r}_{ij})\widehat{\bm{r}}\ dl$
(42)
$\displaystyle\sim-\alpha\oint_{C_{i}}R\bm{\nabla}c(\bm{r}_{ij})\cdot\hat{\bm{r}}\hat{\bm{r}}dl,$
(43)

where $C_{i}$ is the contact line on the surface of the $i$-th particle. We
keep the leading order term in the expansion of the solution with respect to
the ratio between the particle radius and the inter-particle distance
$R/r_{ij}$. Based on this leading order approximation, the inter-particle
force is calculated as follows

$\displaystyle\bm{F}^{\text{M-dist}}_{j\to i}$ $\displaystyle=-\frac{\alpha
AR^{2}}{2D}\exp{\left[\frac{\bm{v}^{f}_{i}\cdot\bm{r}_{ij}}{2D}\right]}$
$\displaystyle\
\times\left[K_{0}(W_{j}r_{ij})\frac{\bm{v}^{f}_{i}}{2D}-W_{j}K_{1}(W_{j}r_{ij})\hat{\bm{r}}_{ij}\right],$
(44)

where $K_{0}(x)$ and $K_{1}(x)$ are the $0$-th and $1$-st order modified
Bessel functions of the second kind, respectively.

###  B.3 Repulsive interaction caused by the Marangoni flow

The Marangoni flow is induced by the unbalance of the Marangoni stress on the
air-liquid interface, and is expressed as follows[40];

$\displaystyle\bm{u}(\bm{r})=\int
d\bm{r}^{\prime}2\bm{G}_{0}(\bm{r}-\bm{r}^{\prime})\cdot\left[-\alpha\bm{\nabla^{\prime}}_{\parallel}\Phi(\bm{r}^{\prime})\right],$
(45)

where $\nabla_{\parallel}\equiv(\partial_{x},\partial_{y})$,
$2\bm{G}_{0}(\bm{r})$ denotes Oseen tensor on the air-liquid interface, and
$\Phi(\bm{r})=\sum_{i}c(\bm{r}-\bm{r}^{p}_{i})$ represents the total
concentration field at the location $\bm{r}$. We can derive the Marangoni flow
by solving the reaction-diffusion equation of surfactant Eq. (3) within the
Stokes approximation as

$\displaystyle u(\bm{r})=\sum_{i}\frac{A\alpha}{8\mu
D}f\left(\frac{|\bm{r}-\bm{r}_{i}^{p}|}{\lambda}\right),$ (46)

where $f(0)=0$ and

$\displaystyle f(x)=\frac{L_{1}(x)+L_{-1}(x)}{2}-I_{1}(x)+\frac{1}{\pi}$ (47)

for $x\neq 0$, and $L_{\nu}(x)$ is the modified Struve functions of order
$\nu$. We can confirm that the Marangoni flow at long distances decays as
$1/r^{2}$ because $\bm{u}(\bm{r})=\int
d\bm{r}^{\prime}\bm{G}_{0}(\bm{r}-\bm{r}^{\prime})\cdot\bm{\nabla^{\prime}}_{\parallel}\delta(\bm{r}^{\prime})=-\bm{\nabla}_{\parallel}\bm{G}_{0}(\bm{r})\propto
1/r^{2}$.

##  Appendix C Validation of the model

For a comparison between the present study and the previous studies[20], we
consider a case where the capillary interaction is neglected, i.e.
$\gamma_{0}^{*}=0.0$. (For the effect of capillary interaction, see Fig. 3.)
Figure 7 shows the dependences of individual forces on the particle distance
for the case with $D^{*}=0.15$. The self-propelling force $F^{\text{M-self}}$
is almost zero for all distances, while the inter-particle force
$F^{\text{M-dist}}$ mediated by the surfactant concentration field is large
and attractive for short distances but is screened beyond the diffusion length
$\lambda=\sqrt{D^{*}/\kappa^{*}}\sim 3.87$. The repulsive force by the
Marangoni flow gives the largest contribution to the total force, leading to
the long-range repulsive interaction between particles.

Figure 7:  Distance dependences of the individual forces for the case with
$D^{*}=0.15$, $\text{Ma}=1.0$, $\kappa^{*}=0.01$, and $\gamma_{0}^{*}=0.0$.
The black, pink, green, red, and blue lines represent total force acting on
the particle, self-propelling force $F^{\text{M-self}}$, the inter-particle
interaction force mediated by the surfactant concentration field
$F^{\text{M-dist}}$, viscous friction force by the hydrodynamic interaction,
and repulsive force by Marangoni flow, respectively. The dots show the data
from the previous study by Soh et al.[20] and the horizontal axis is in a
logarithmic scale.

## References

  * Vicsek and Zafeiris [2012] T. Vicsek and A. Zafeiris, Phys. Rep. 517, 71 (2012).
  * Marchetti  _et al._ [2013] M. C. Marchetti, J.-F. Joanny, S. Ramaswamy, T. B. Liverpool, J. Prost, M. Rao, and R. A. Simha, Rev. Mod. Phys. 85, 1143 (2013).
  * Gompper  _et al._ [2020] G. Gompper, R. G. Winkler, T. Speck, A. Solon, C. Nardini, F. Peruani, H. Löwen, R. Golestanian, U. B. Kaupp, L. Alvarez, _et al._ , J. Phys. Condens. Matter 32, 193001 (2020).
  * Popescu [2020] M. N. Popescu, Langmuir 36, 6861 (2020).
  * Lauga and Powers [2009] E. Lauga and T. R. Powers, Rep. Prog. Phys. 72, 096601 (2009).
  * Koch and Subramanian [2011] D. L. Koch and G. Subramanian, Ann. Rev. Fluid Mech. 43, 637 (2011).
  * Paxton  _et al._ [2004] W. F. Paxton, K. C. Kistler, C. C. Olmeda, A. Sen, S. K. St. Angelo, Y. Cao, T. E. Mallouk, P. E. Lammert, and V. H. Crespi, J. Am. Chem. Soc. 126, 13424 (2004).
  * Ebbens and Howse [2010] S. J. Ebbens and J. R. Howse, Soft Matter 6, 726 (2010).
  * Bechinger  _et al._ [2016] C. Bechinger, R. Di Leonardo, H. Löwen, C. Reichhardt, G. Volpe, and G. Volpe, Rev. Mod. Phys. 88, 045006 (2016).
  * Zöttl and Stark [2016] A. Zöttl and H. Stark, J. Phys. Condens. Matter 28, 253001 (2016).
  * Kohira  _et al._ [2001] M. I. Kohira, Y. Hayashima, M. Nagayama, and S. Nakata, Langmuir 17, 7124 (2001).
  * Nagayama  _et al._ [2004] M. Nagayama, S. Nakata, Y. Doi, and Y. Hayashima, Physica D 194, 151 (2004).
  * Kitahata  _et al._ [2004] H. Kitahata, S.-i. Hiromatsu, Y. Doi, S. Nakata, and M. R. Islam, Phys. Chem. Chem. Phys. 6, 2409 (2004).
  * Nakata  _et al._ [2015a] S. Nakata, M. Nagayama, H. Kitahata, N. J. Suematsu, and T. Hasegawa, Phys. Chem. Chem. Phys. 17, 10326 (2015a).
  * Boniface  _et al._ [2019] D. Boniface, C. Cottin-Bizonne, R. Kervil, C. Ybert, and F. Detcheverry, Phys. Rev. E 99, 062605 (2019).
  * Hanczyc  _et al._ [2007] M. M. Hanczyc, T. Toyota, T. Ikegami, N. Packard, and T. Sugawara, J. Am. Chem. Soc. 129, 9386 (2007).
  * Toyota  _et al._ [2009] T. Toyota, N. Maru, M. M. Hanczyc, T. Ikegami, and T. Sugawara, J. Am. Chem. Soc. 131, 5012 (2009).
  * Tanaka  _et al._ [2015] S. Tanaka, Y. Sogabe, and S. Nakata, Phys. Rev. E 91, 032406 (2015).
  * Maass  _et al._ [2016] C. C. Maass, C. Krüger, S. Herminghaus, and C. Bahr, Annu. Rev. Condens. Matter Phys. 7, 171 (2016).
  * Soh  _et al._ [2008] S. Soh, K. J. Bishop, and B. A. Grzybowski, J. Phys. Chem. B 112, 10848 (2008).
  * Ikura  _et al._ [2013] Y. S. Ikura, E. Heisler, A. Awazu, H. Nishimori, and S. Nakata, Phys. Rev. E 88, 012911 (2013).
  * Nishimori  _et al._ [2017] H. Nishimori, N. J. Suematsu, and S. Nakata, J. Phys. Soc. Jpn. 86, 101012 (2017).
  * Nagai  _et al._ [2005] K. Nagai, Y. Sumino, H. Kitahata, and K. Yoshikawa, Phys. Rev. E 71, 065301(R) (2005).
  * Chen  _et al._ [2009] Y.-J. Chen, Y. Nagamine, and K. Yoshikawa, Phys. Rev. E 80, 016303 (2009).
  * Tanaka  _et al._ [2017] S. Tanaka, S. Nakata, and T. Kano, J. Phys. Soc. Jpn. 86, 101004 (2017).
  * Kano  _et al._ [2017] T. Kano, K. Osuka, T. Kawakatsu, and A. Ishiguro, J. Phys. Soc. Jpn. 86, 124004 (2017).
  * Yabunaka  _et al._ [2012] S. Yabunaka, T. Ohta, and N. Yoshinaga, J. Chem. Phys. 136, 074904 (2012).
  * Ikura  _et al._ [2012] Y. S. Ikura, R. Tenno, H. Kitahata, N. J. Suematsu, and S. Nakata, J. Phys. Chem. B 116, 992 (2012).
  * Masoud and Shelley [2014] H. Masoud and M. J. Shelley, Phys. Rev. Lett. 112, 128304 (2014).
  * Nishi  _et al._ [2015] K. Nishi, K. Wakai, T. Ueda, M. Yoshii, Y. S. Ikura, H. Nishimori, S. Nakata, and M. Nagayama, Phys. Rev. E 92, 022910 (2015).
  * Nakata  _et al._ [2015b] S. Nakata, R. Tenno, A. Deguchi, H. Yamamoto, Y. Hiraga, and S. Izumi, Colloids Surf. A Physicochem. Eng. Asp. 466, 40 (2015b).
  * Matsuda  _et al._ [2016] Y. Matsuda, N. J. Suematsu, H. Kitahata, Y. S. Ikura, and S. Nakata, Chem. Phys. Lett. 654, 92 (2016).
  * Domínguez  _et al._ [2016a] A. Domínguez, P. Malgaretti, M. Popescu, and S. Dietrich, Soft Matter 12, 8398 (2016a).
  * Yabunaka and Yoshinaga [2016] S. Yabunaka and N. Yoshinaga, J. Fluid Mech. 806, 205 (2016).
  * Kitahata and Yoshinaga [2018] H. Kitahata and N. Yoshinaga, J. Chem. Phys. 148, 134906 (2018).
  * Hirose  _et al._ [2020] Y. Hirose, Y. Yasugahira, M. Okamoto, Y. Koyano, H. Kitahata, M. Nagayama, and Y. Sumino, J. Phys. Soc. Jpn. 89, 074004 (2020).
  * Rotne and Prager [1969] J. Rotne and S. Prager, J. Chem. Phys. 50, 4831 (1969).
  * Yamakawa [1970] H. Yamakawa, J. Chem. Phys. 53, 436 (1970).
  * Perkins and Jones [1991] G. Perkins and R. Jones, Physica A 171, 575 (1991).
  * Domínguez  _et al._ [2016b] A. Domínguez, P. Malgaretti, M. N. Popescu, and S. Dietrich, Phys. Rev. Lett. 116, 078301 (2016b).
  * Nakata  _et al._ [2000] S. Nakata, Y. Hayashima, and H. Komoto, Phys. Chem. Chem. Phys. 2, 2395 (2000).
  * Hayashima  _et al._ [2001] Y. Hayashima, M. Nagayama, and S. Nakata, J. Phys. Chem. B 105, 5353 (2001).
  * Kralchevsky and Nagayama [1994] P. A. Kralchevsky and K. Nagayama, Langmuir 10, 23 (1994).
  * Weeks  _et al._ [1971] J. D. Weeks, D. Chandler, and H. C. Andersen, J. Chem. Phys. 54, 5237 (1971).
  * Dörr  _et al._ [2016] A. Dörr, S. Hardt, H. Masoud, and H. A. Stone, J. Fluid Mech. 790, 607 (2016).
  * Suematsu  _et al._ [2014] N. J. Suematsu, T. Sasaki, S. Nakata, and H. Kitahata, Langmuir 30, 8101 (2014).
  * Michelin  _et al._ [2013] S. Michelin, E. Lauga, and D. Bartolo, Phys. Fluids 25, 061701 (2013).
  * [48] See Supplemental Material at [URL will be inserted by publisher] for a movie showing the collective motion in cases (a) through (d) for FIG. 5. .
  * Bialké  _et al._ [2015] J. Bialké, T. Speck, and H. Löwen, J. Non-Cryst. Solids 407, 367 (2015).