# Collective Motion of Quincke Rollers with Fully Resolved Hydrodynamics

Shun Imamura1 Kohei Sawaki1 John J. Molina1 Matthew S. Turner1,2 Ryoichi
Yamamoto1<EMAIL_ADDRESS>1Department of Chemical Engineering,
Kyoto University, Kyoto 615-8510, Japan
2Department of Physics, University of Warwick, Coventry CV4 7AL, UK

###### Abstract

A Quincke roller is a unique active particle that can run and tumble freely on
a flat plate due to the torque generated by a uniform DC electric field
applied perpendicular to the plate[1]. A system involving many such particles
exhibits a variety of collective dynamics, such as the disordered gas, polar
liquid, and active crystal states[1, 2, 3, 4, 5, 6, 7, 8]. We performed direct
numerical simulations of a three-dimensional system containing many self-
rotating particles to explicitly resolve the hydrodynamic interactions among
rotating particles. The collective motion depends on the magnitude of the
dipole moments induced on the dielectric particles, the area fraction of
particles, and the strength of interparticle attraction. We find that the
highly ordered polar liquid state is destabilized by the hydrodynamic
interaction between rotating particles at high densities: the near-field
lubrication interaction becomes dominant over far-field effects as the
interparticle separation becomes shorter.

When a uniform DC electric field is applied to a dielectric particle in an
electrolyte fluid, spontaneous rotation of the particle can be observed due to
an electrohydrodynamic instability known as the Quincke effect [9] (see Fig.
1). A Quincke roller is a self-propelled particle whose propulsion is driven
by the hydrodynamic traction force exerted on a particle rotating due to the
Quincke effect near a substrate [1]. The propulsion velocity of the roller is
proportional to the angular velocity, which can be controlled via the
magnitude of the external electric field [10, 11].

A variety of distinct regimes have been observed in experimental studies,
including gas, polar liquid, active solid and cluster states [1, 3, 4, 5, 6,
7, 8], as well as an order–disorder transition between the polar liquid and
gas states [1], mobility enhancement due to paired-up states [3], and
solidification at high density [5]. Numerical studies have also been
conducted, mostly using two-dimensional particle-based models. These have had
partial success in reproducing some experimental observations [1, 6]. The far-
field hydrodynamic interaction (HI) between Quincke rollers has been
considered using Stokesian dynamics [12, 13, 14]. A continuum model has also
been used to investigate the linear stability of the liquid–solid flocking
transition [1, 5]. In previous numerical studies the near-field HI has been
greatly simplified, even for systems at high densities where a more faithful
treatment of the lubrication effect between rotating particles becomes
essential.

Figure 1:  A schematic illustration of a Quincke roller. When a DC electric
field $E_{0}$ greater than a threshold $E_{Q}$ is applied to a dielectric
particle located near the lower of two plates in an electrolyte fluid, the
particle spontaneously rotates and is propelled to move on the plate due to an
electrohydrodynamic instability called the Quincke effect. Here, $\bm{P}$ is
the induced electric dipole moment, $\bm{N}$ is the torque exerted due to the
Quincke effect, $v_{0}$ $(\propto E_{0})$ is the propulsion velocity, and $H$
is the distance between the lower and upper electrode plates.

In this paper, we report direct numerical simulations (DNS) of a three-
dimensional system containing many self-rotating particles in which we
explicitly resolve the HIs between them. We first present the details of our
numerical implementation method. Then we test the validity of the present DNS
method by comparing our results with theoretical predictions for a single
particle system. The role and importance of the near-field hydrodynamics is
then examined for two- and three-particle systems. Finally, we present a
comprehensive analysis of our DNS results for many-particle systems. Here the
collective motion depends on the magnitude of the induced dipole moments, the
area fraction of particles, and the strength of the interparticle attraction.

Many computational simulation methods that include HIs have been proposed for
systems of microswimmers, including colloidal rollers [15]. One example is the
smoothed profile (SP) method. This has been applied to study the structure and
dynamics in a variety of particle dispersions, such as colloidal systems and
self-propelled microswimmers [16]. The SP method can describe both near-field
and far-field hydrodynamic effects reasonably accurately and is applicable to
dense systems [17, 18]. Here, we use dimensionless simulation units throughout
the length given by the grid spacing, the time and the energy by the fluid
density and the fluid viscosity, and the electric dipole by the dielectric
constant in vacuum (see Methods for more details).

First, as an evaluation of the SP method, we confirm that the speed of the
single roller is in good agreement with lubrication theory [19, 20], and that
the flow field around the roller agrees with previous simulations that used a
rotlet representation corrected by the Blake tensor [21, 13, 14, 22] (see
Supplementary Information). The advantage of the SP method is that it can
accurately handle HIs between rollers. We first focus on how the cluster
states of two and three rollers are affected by these lubrication effects
before later studying many-body systems.

Figure 2:  Behavior of a cluster consisting of two rollers. (a) The viscous
dissipation rate per volume,
$w=[(\bm{\nabla}\bm{u})+(\bm{\nabla}\bm{u})^{T}]^{2}/2$, for the face-to-back
($\theta=0$) and back-to-back ($\theta=-\pi$) states for radius $a=8$
particles, and the blue circles represent the boundary line of each particle.
The viscous dissipation in the space between the two rollers is higher than it
is elsewhere, and the dissipation in the face-to-back state is larger than
that in the back-to-back state. (b) The counter-rotating (face-to-face or
back-to-back) states dissipate less heat and might therefore be expected to be
more stable than the co-rotating (face-to-back) state, based on the principle
of minimum energy dissipation. Here, the face-to-face (back-to-back) state
represents the case in which the two rollers are being propelled toward (away
from) each other. Results shown for different roller sizes $a=2$, $4$, $6$,
and $8$. (c) Stable steady state for different dipole strengths draws a
hysteresis curve, where $a=2$. Arrows and colors in this curve indicate the
direction of the loop. The HI stabilizes the back-to–back motion $\theta=-\pi$
at $\mu_{\parallel}$ is small, and the electrostatic interaction stabilizes
the face-to-back motion $\theta=0$ at $\mu_{\parallel}$ is large. Here, the
strength of attraction sets $\epsilon=10$, and all units are dimensionless,
see Methods for more details.

The stable cluster states for two-body systems can be classified into a co-
rotating (face-to-back) state, in which the rollers are propelled in the same
direction, $\theta=0$, and a counter-rotating (face-to-face $\theta=\pi$ or
back-to-back $\theta=-\pi$) state, in which they rotate in opposite
directions, where $\theta$ is the relative angle between the axes of Quincke
rotation of the rollers (see Supplementary Movie 1). In order to understand
the stabilizing structure of clusters by HI, we consider the viscous
dissipation of the system. The viscous dissipation of the fluid around the
clusters in the co-rotating and counter-rotating states is shown in Fig. 2
(a). The viscous dissipation rate is defined as $W=\int w(\bm{r})d\bm{r}$
where the density
$w(\bm{r})=\left[(\bm{\nabla}\bm{u}(\bm{r}))+(\bm{\nabla}\bm{u}(\bm{r}))^{T}\right]^{2}/2$.
According to the principle of minimum energy dissipation [23, 24, 25], a
nonequilibrium system prefers to be in a state with lower energy dissipation.
In this system, the viscous dissipation in a counter-rotating state is lower
than that in the co-rotating state, see Fig. 2 (b). In addition, from the
discussion regarding linear stability [26] and the forces acting between the
rollers, we find that the back-to-back state is more stable than the face-to-
back and face-to-face states due to the HI (see Supplementary Information). On
the other hand, the face-to-back state is stabilized by electrostatic
interactions at the larger dipole strength. This transition between face-to-
back and back-to-back exhibits hysteresis-like behavior surprisingly, see Fig.
2 (c).

Figure 3:  The rotational order parameter $\Psi$ in a three-body system. When
the dipole strength $\mu_{\parallel}$ is small, the HI propels the particles
outward from each other, stabilizing the immobile state of the cluster. On the
other hand, when $\mu_{\parallel}$ is large, the dipole interaction stabilizes
the rotating cluster state.

The behavior of the three-roller cluster can be explained by extending the
analysis of the two-roller case. When the dipole interaction is weak, the HI
stabilizes the immobile state, but when the dipole interaction is strong, the
rotating cluster state is stabilized (see Fig. 3 and Supplementary Movie 2).
This is a result of the dipole interaction stabilizing the rotation axis of
each particle either inward or outward relative to the geometric center of the
cluster. To evaluate this quantitatively, we introduce a rotational order
parameter for each cluster,
$\Psi=\braket{}{\bm{n}_{i}\cdot\widehat{{\bm{r}}^{\prime}}_{i}}{}$, where
$\braket{\cdots}$ denotes an average over the particles and time, the vector
symbol $\widehat{\cdot}$ represents a unit vector, the unit vector
$\bm{n}_{i}$ refers to the axis of Quincke rotation for the $i$-th roller and
${\bm{r}}^{\prime}_{i}=\bm{r}_{i}-\bm{r}_{G}$, with $\bm{r}_{G}$ being the
center of mass of the cluster (see Supplementary Information). A continuous
transition is observed for the three-body cluster, between the immobile
disordered clusters and the rotating clusters.

Figure 4:  Phase diagram of a many-roller system. The blue-outlined region
represents polar liquid states, the gray-outlined region represents gas
states, the green-outlined region represents vortex states, the purple-
outlined region represents chain states, the orange-outlined region represents
rotating cluster states, the olive-outlined region represents disordered
cluster states, and the azure-outlined region represents polar cluster states.
See Supplementary Movies 3 and 4 for representative examples of all of these
phases.

In our simulations for a many-roller systems, we use steady-state solutions
for the dipole and angular velocity, which depend on the electric field
$E\equiv E_{0}/E_{Q}$[27, 1] (see Methods). The collective structures realized
by the many-roller system can be classified by a set of order parameters [28]
as follows: polar liquid states (blue-outlined region), gas states (gray),
vortex states (green), chain states (purple), rotating cluster states
(orange), disordered cluster states (olive), and polar cluster states (azure)
(see Fig. 4, Supplementary Information and Supplementary Movies 3 and 4). When
the applied electric field $E$ is small, each cluster is a rotating cluster
state in the lower-density region, but as the density increases, the clusters
tend to become more anisotropic, colliding and coalescing repeatedly to form
disordered clusters. At high enough density, all rollers form a single polar
cluster. On the other hand, when the applied electric field $E$ is large, as
the area fraction increases, we can find a transition from an ordered polar
state to a disordered gas or vortex state. When the repulsive force exerted by
$\mu_{\perp}$ and the attractive force are of the same order, the
dipole–dipole interaction given by $\mu_{\parallel}$ dominates, and a chain
state is formed.

Our new three-dimensional direct numerical simulation provides important
insights into the relationship between the interactions and the collective
dynamics. First, the differences between our observed phase-diagram and
previous simulation and experimental results can mainly be traced to the
breaking of polar liquid states at higher density[1, 4, 6]. As the applied
electric $E_{0}$ field increases, the repulsion due to $\mu_{\perp}$ exceeds
the attraction due to the electro-osmotic flow, causing the rollers to move
apart. Then, as the area fraction $\phi_{0}$ increases, the average distance
between the rollers changes, balancing the HIs and dipole interactions and
causing a transition from a disordered gas state to an ordered polar liquid
state, similar to the results by Bricard et al. [1]. As the area fraction
$\phi_{0}$ increases further, the interactions become unbalanced again, and
transitions to the disordered gas and vortex states from polar liquid states
are observed at $\phi_{0}\geq 0.2$ (see Supplementary Information). Our
results suggest that the disordered dynamics of the active solid state
observed by Geyer et al. [5] are caused by the polar state disappearing at
high density due to the lubrication effects. However, unlike in the
experiments, the solidification from polar liquid states has not been observed
in our simulations, which we believe is due to implementation differences,
e.g. different boundary conditions, system size and number of particles.
Second, we expect our simulations to reach the chain state whenever the dipole
$\mu_{\parallel}$ becomes dominant, in where isotropic forces between the
electrostatic forces by $\mu_{perp}$ and electro-osmotic attractive forces are
balanced at $E\sim 2.0$. Such a chain state has been observed actually under
different conditions with the proceeding researches[1, 6], including larger
colloid size and different concentration of conducting liquid
experimentally[29], and the active Brownian particles with dipole–dipole
interactions numerically[30]. Finally, they provide new insights into the
physical mechanisms responsible for various cluster motions observed in
previous experiments [4, 6, 8]. When the applied electric field $E_{0}$ is
small, the attraction due to the electro-osmotic flow is dominant and the
cluster motion is reproduced. A rotating cluster created by stabilization
through $\mu_{\parallel}$, and clusters with anisotropic shapes (e.g.,
elliptical), which tend to form at high densities, as well as select
translational motion in a polar cluster state. A disordered cluster arises due
to the lubrication effect and the shape change caused by frequent coalescence
and separation between clusters.

In conclusion, we take advantage of DNSs to include HI and investigate the
flocking transition mechanism of Quincke rollers and their dynamics in high-
density regions. The DNS of the fluid flow is performed using the SP method.
Our simulations reveal that the transition between cluster states is caused by
the balance between the disordering of the orientation due to the lubrication
effects and the ordering of the orientation due to the horizontal component of
the electric dipole, and the deformation of the cluster shape. Additionally,
our simulations show an intermediate transition to disordered states from
polar liquid states at high densities, due to an imbalance between
hydrodynamic and dipole interactions. Our proposed DNS method can be easily
extended to the collective dynamics of colloidal rollers in magnetic fields
[31, 14]. Our findings show that the HIs differ between the near and far
fields, which indicates the possibility of new improvements in particle-based
models.

### Acknowledgements

This work was supported by Grants-in-Aid for Scientific Research (JSPS
KAKENHI) under grant nos. JP 20H00129, 20H05619, and 20K03786. It was also
supported by Professional development Consortium for Computational Materials
Scientists (PCoMS). R.Y. acknowledges helpful discussions with Profs. Hajime
Tanaka and Akira Furukawa.

### Author contributions

All authors designed the research and wrote the paper. S.I., K.S. and R.Y.
performed the simulations and analysed the data. J.J.M. assisted with
numerical methods.

### Competing interests

The authors declare no competing interests.

### Additional information

Supplementary Information is available for this paper at [url will be inserted
by publisher].

### Methods

#### Smoothed Profile (SP) method

The presence of the $i$-th particle is described by the SP function
$\phi_{i}$, which takes values of $\phi_{i}=1$ in the solid domain and
$\phi_{i}=0$ in the fluid domain. These domains are smoothly connected by an
interface of finite width $\xi$ at the particle surface [17]. In the SP
method, the velocity field $\bm{u}$ is defined at all positions $\bm{x}$ and
times $t$ so as to interpolate between the host fluid, moving with velocity
$\bm{u}_{f}$, and the solid particles moving with velocity $\bm{u}_{p}$,

$\bm{u}(\bm{x},t)=(1-\phi)\bm{u}_{f}(\bm{x},t)+\phi\bm{u}_{p}(\bm{x},t).$ (1)

In the equation above, the particle field $\phi$ is to be treated as a sum
over the contribution of all particles. The two terms in (1) are defined as

$(1-\phi)\bm{u}_{f}(\bm{x},t)=\left(1-\sum^{N}_{i=1}\phi_{i}(\bm{x},t)\right)\bm{u}_{f}(\bm{x},t)$
(2)

and

$\phi\bm{u}_{p}(\bm{x},t)=\sum^{N}_{i=1}\phi_{i}(\bm{x},t)(\bm{V}_{i}+\bm{\Omega}_{i}\times\bm{r}_{i})$
(3)

where $\bm{r}_{i}=\bm{x}-\bm{R}_{i}$ is a radial vector originating at the
center of the $i$-th particle, with $\bm{R}_{i}$ that particle’s center,
$\bm{V}_{i}$ its velocity, and $\bm{\Omega}_{i}$ its angular velocity. The
total particle phase-field $\phi(\bm{x},t)$ is defined as the superposition of
the SP functions of all particles. The total velocity field is governed by a
modified Navier–Stokes equation, together with the incompressibility
condition:

$\displaystyle\rho_{f}(\partial_{t}+\bm{u}\cdot\bm{\nabla})\bm{u}$
$\displaystyle=-\bm{\nabla}p+\eta_{f}\bm{\nabla}^{2}\bm{u}+\rho_{f}\phi\bm{f}_{p},$
(4) $\displaystyle\bm{\nabla}\cdot\bm{u}$ $\displaystyle=0,$ (5)

where $\rho_{f}$ and $\eta_{f}$ are the density and viscosity of the fluid,
respectively, and $p$ is the pressure. The rigidity of the particle domains is
maintained by the force density field $\phi\bm{f}_{p}$, which is computed to
ensure momentum conservation between the fluid domain and the rigid-body
domain of the particles [17].

The motion of the Quincke rollers is determined by the Newton–Euler equations
of motion:

$\displaystyle\frac{d\bm{R}}{dt}$ $\displaystyle=\bm{V},$ (6) $\displaystyle
M\frac{d\bm{V}}{dt}$ $\displaystyle=\bm{F}^{H}+\bm{F}^{P}+\bm{F}^{E},$ (7)
$\displaystyle\dot{\bm{Q}}$
$\displaystyle=\mathrm{skew}(\bm{\Omega})\cdot\bm{Q}$ (8)
$\displaystyle\bm{I}\cdot\frac{d\bm{\Omega}}{dt}$
$\displaystyle=\bm{N}^{H}+\bm{N}^{P}+\bm{N}^{C},$ (9)

where $M=4\pi a^{3}\rho_{p}/3$ is the mass of a particle, $\bm{Q}$ the
particle’s rotation matrix ($\mathrm{skew}(\bm{\Omega})$ the skew-symmetric
angular velocity matrix), $\bm{I}=2a^{2}M\mathsf{I}/5$ the moment of inertia
($\mathsf{I}$ is the unit tensor), $a$ the particle radius and $\rho_{p}$ the
particle density. The hydrodynamic forces and torques due to the fluid-
particle interactions are denoted as $\bm{F}^{H}$ and $\bm{N}^{H}$,
respectively [17]. The forces and torques due to the inter-particle
interactions, $\bm{F}^{P}$ and $\bm{N}^{P}$, are derived from an interaction
potential $U=U^{S}+U^{A}+U^{D}$, which is assumed to be pairwise additive[32]
(see Supplementary Information). The corresponding potential between a pair of
particles $i$ and $j$, with dipole moments $\bm{P}_{i}$ and $\bm{P}_{j}$ and
separation $\bm{r}=\bm{R}_{j}-\bm{R}_{i}$ is

$\displaystyle U_{ij}^{S}(r)$
$\displaystyle=4\epsilon^{LJ}\left[\left(\frac{\sigma}{r}\right)^{36}-\left(\frac{\sigma}{r}\right)^{18}\right]+\epsilon^{LJ},$
(10) $\displaystyle U_{ij}^{A}(r)$ $\displaystyle=-\epsilon^{A}\exp{(-\kappa
r)}/r^{2},$ (11) $\displaystyle U_{ij}^{D}(\bm{r},\bm{P}_{i},\bm{P}_{j})$
$\displaystyle=\frac{1}{4\pi\epsilon_{0}}\left[\frac{\bm{P}_{i}\cdot\bm{P}_{j}}{r^{3}}-3\frac{(\bm{P}_{i}\cdot\bm{r})(\bm{P}_{j}\cdot\bm{r})}{r^{5}}\right].$
(12)

Here the $U_{ij}^{S}$ potential represents the excluded volume effect, and is
described by a truncated $18-36$ Lennard–Jones potential, with a cut-off at
$r_{c}=2^{1/18}\sigma$, with $\sigma=2a$ the particle diameter. The
$U_{ij}^{A}$ potential is used to describe the electro-osmotic attraction,
where the screening length is given by $\kappa=1/(3\sigma)$ [6]. The strengths
of the Lennard-Jones repulsion and electro-osmotic attraction are given by
$\epsilon^{LJ}$ and $\epsilon^{A}$, respectively. The $U_{ij}^{D}$ represents
the electrostatic dipole–dipole interactions and includes the interactions of
the image dipoles [1]. To mimic experimental conditions, we set the applied DC
electric field to be perpendicular to the walls, i.e.,
$\bm{E}_{0}=E_{0}\bm{e}_{z}$ where $\bm{e}_{z}$ is an unit vector of the $Z$
axis, and we set also the angular velocity vector of a single roller due to
its Quincke rotation $\bm{\Omega}_{0}=\Omega_{0}\bm{n}$, where the initial
orientation is randomly chosen to be parallel to the walls ($XY$-plane), i.e.,
at the time the electric field is applied, where the magnitude $\Omega_{0}$
depends on the scaled electric field $E_{0}$ (see Eq. (14)). The electric
dipole is defined as $\bm{P}=P_{\parallel}\bm{e}_{p}+P_{\perp}\bm{e}_{z}$,
where $P_{\parallel}$ and $P_{\perp}$ are the components parallel and
perpendicular to the wall, respectively, and $\bm{e}_{p}$ is a unit vector in
the direction of $\bm{E}_{0}\times\bm{\Omega}_{0}$. These electrostatic
interactions are calculated via the Ewald method [33, 34, 35]. Here,
$P_{\parallel}$ gives the dipole–dipole interaction on the plane, and
$P_{\perp}$ gives an isotropic repulsion [1]. $\bm{F}^{E}$ represents the
(external) forces exerted by gravity and the wall, where the wall–particle
interaction is again represented by a truncated Lennard–Jones potential,
similar to the particle-particle steric interactions, that is chosen to be
strong enough to exclude the particles from the wall. Finally, to implement
the Quincke roller dynamics, we introduce a constraint torque $\bm{N}^{C}$,
derived from a potential $U^{C}$[36], to fix the rotation denoted by
$\bm{\Omega}_{0}$ to lie in the plane (see Supplementary Information).

#### Parameters and Units

##### The electric dipole moments and the rotation speed by Quincke rotation
for the steady-state solution:

Let us consider the correspondence between the theoretical and simulation
parameters. In Quincke rotation, in accordance with the time evolution of an
electric dipole[27, 1], the steady-state solution for the dipole moment
induced in a dielectric sphere can be expressed as follows:

$\displaystyle
P_{\perp}=P_{\perp}^{\epsilon}E-\frac{P_{\perp}^{\sigma}}{E},P_{\parallel}=P_{\perp}^{\sigma}\sqrt{1-\frac{1}{E^{2}}},$
(13)

where $E\equiv E_{0}/E_{Q}$, $P_{\perp}^{\epsilon}\equiv
4\pi\epsilon_{0}a^{3}\chi^{\infty}E_{Q}$ is the static contribution from the
dielectric polarization, and the $z$ component
$P_{\perp}^{\sigma}\equiv\epsilon_{0}/(\tau\epsilon_{l}\mu_{r}E_{Q})$ and the
$xy$ component $P_{\parallel}$ are the dynamic contributions from the Quincke
rotation for $E>1$. $E_{Q}$ is the threshold value for Quincke rotation,
$\tau$ is the Maxwell–Wagner time, $\mu_{r}=(8\pi\eta_{f}a^{3})^{-1}$ is the
mobility of the rotational motion, and
$\chi^{\infty}\equiv\frac{\epsilon_{p}-\epsilon_{l}}{\epsilon_{p}+2\epsilon_{l}}$,
with $\epsilon_{p}$ and $\epsilon_{l}$ the dielectric permittivities of the
sphere and the liquid medium, respectively, and $\epsilon_{0}$ the dielectric
permittivity of vacuum. Additionally, the angular velocity of a roller in the
steady state is described as

$\displaystyle\Omega_{0}=\frac{1}{\tau}\sqrt{E^{2}-1}.$ (14)

Note that the accuracy of approximation by this steady-state solutions is
guaranteed by the fact that the unit of time $t_{0}$ is sufficiently longer
than the Maxwell–Wagner time $\tau$, e.g. $t_{0}\gg\tau$ (see the paragraph of
units about a definition of the unit $t_{0}$).

##### Units:

Within the SP method, the system is discretized on a rectangular grid with
spacing $\Delta$, under periodic boundary conditions in all directions. This
allows us to employ a pseudo-spectral method, and significantly reduce the
computational cost. The units of all variables are based on the fundamental
simulation parameters $\Delta$, $\eta_{f}$ and $\rho_{f}$, the grid spacing,
fluid viscosity, and fluid density respectively. The unit of length is
$l_{0}=\Delta$, time $t_{0}=\rho_{f}\Delta^{2}/\eta_{f}$, and that of energy
is $e_{0}=\eta_{f}^{2}\Delta/\rho_{f}$. The electric dipole is measured in
units of $p_{0}=\sqrt{4\pi\epsilon_{0}e_{0}l_{0}^{3}}$. Here, the viscosity of
the liquid is $\eta_{f}\sim 4.3\times 10^{-3}$ Pa $\cdot$ s, the density ratio
between the colloids and the liquid is $\rho_{p}/\rho_{f}=1.19/0.77\sim
1.5$[11], the dielectric permittivities $\epsilon_{p}\sim
2.6\epsilon_{0}$[10], $\epsilon_{l}\sim 2.0\epsilon_{0}$[37], and the ratio
$\chi^{\infty}\sim 0.08$[1]. Based on these physical quantities, we can
estimate the parameter units to $l_{0}=\Delta=1$ $\mu$m,
$t_{0}=\frac{\rho_{f}\Delta^{2}}{\eta_{f}}\sim 0.1$ $\mu$s,
$e_{0}=\frac{\eta_{f}^{2}\Delta}{\rho_{f}}\sim 10^{-14}$ J, and
$p_{0}=\sqrt{4\pi\epsilon_{0}e_{0}l_{0}^{3}}\sim 10^{-21}$ C$\cdot$m. Because
the Maxwell-Wagner time $\tau\sim 1$ ms[1], the accuracy of the steady-state
solutions Eqs. (13) and (14) is not sufficiently guaranteed in this unit
$t_{0}$. However, to simplify the problem and to build our fundamental
understanding in the collective system, we proceed with the simulation under
the assumption that the contribution from the steady-state is dominant.

##### Simulation parameters:

We set the simulation parameters from these theoretical values, and we define
$\mu_{\parallel}$ and $\mu_{\perp}$ as the dimensionless dipole components
corresponding to $P_{\parallel}$ and $P_{\perp}$, respectively, where in what
follows all units are dimensionless. The system size is
$(L_{x},L_{y},L_{z})=(128,128,32)$ and we introduce a flat wall parallel to
the $XY$ plane, with a total depth of $d=5$, i.e., $H$ is fixed for all
simulations. The particle Reynolds number is defined as
$Re\equiv\rho_{f}\sigma^{2}\Omega_{0}/\eta_{f}$, and $Re\sim 10^{-5}$ in
experiments[1, 11], where we assume that $Re\sim 1$ to reduce the
computational cost. Because the excluded volume effect is sufficiently larger
than the electro-osmotic attraction, we set
$\epsilon\equiv\epsilon^{LJ}=\epsilon^{A}$ for simplicity, and we takes the
interfacial width of the particle domain $\xi=2$ with high reliability[17].
For single, two, and three rollers simulation, we set $Re=1$. To investigate
the stable structures of cluster states in the simulation for two and three
rollers, we assume that the strength of repulsive force $\mu_{\perp}=0$ and
the attractive force $\epsilon=10.0$ for two rollers and $\epsilon=30.0$ for
three rollers to keep cluster motion, where the radius $a=2,4,6,$ and $8$ for
the two rollers and $a=2$ for the three rollers. Additionally, for many-
rollers simulation, we set $a=2$, the strength of the attractive force
$\epsilon=10.0$, and the area fraction is varied in the range
$\phi_{0}=0.01-0.8$, based on the steady-state solution Eqs. (13) and (14),
each input parameter is determined from the electric field $E$ as follows: the
angular velocity is $\Omega_{0}\sim 0.025\times\sqrt{E^{2}-1}$, the dipole
moment is $\mu_{\parallel}\sim 2.6\times\sqrt{1-\frac{1}{E^{2}}}$, and the
ratio between the dipole components is $|\mu_{\perp}|/|\mu_{\parallel}|\sim
E$.

Supplemental Information for the Manuscript

“Collective Motions of Quincke Rollers with Fully Resolved Hydrodynamics”

Shun Imamura, Kohei Sawaki, John J. Molina, Matthew S. Turner, and Ryoichi
Yamamoto

Department of Chemical Engineering, Kyoto University, Kyoto 615-8510, Japan

##  I Simulation Details

###  I.1 Derivation of forces and torques from the potential

The potential $U$ depends not only on the relative position of the particles
$\bm{r}$ but also on their relative orientations in space, as given by the
rotation axis for Quincke rotation and the dipole moments of the particles.
The interparticle forces and torques can be obtained from the potential as
follows [32]:

$\displaystyle\bm{F}^{P}$ $\displaystyle=-\frac{\partial U}{\partial\bm{r}},$
(15) $\displaystyle\bm{N}^{P}$
$\displaystyle=-\sum_{m}\widehat{\bm{a}}_{m}\times\frac{\partial
U}{\partial\widehat{\bm{a}}_{m}},$ (16)

where the $\widehat{\bm{a}}_{i}$ are the (body-fixed) ortho-normal principal
axis vectors ($m=x^{\prime},y^{\prime},z^{\prime}$).

###  I.2 A torque to constrain the Quincke rotation

Figure 5:  Schematic diagram of the torque used to constrain the axis of
rotation. The angle between the applied electric field $\bm{E}_{0}$ and the
angular velocity by the Quincke rotation $\bm{\Omega}_{0}$ is constrained to
be perpendicular by a harmonic potential.

We introduce an external torque inspired by a Quincke roller. It has been
analytically proven that the angular velocity vector generated by Quincke
rotation is perpendicular to the external electric field $\bm{E}_{0}$[10]. In
our model, to simplify the treatment, we impose a constant angular velocity
vector $\bm{\Omega}_{0}$ on a rigid sphere and a constraint torque
$\bm{N}^{C}$ defined by a harmonic potential is introduced to maintain the
rotation axis of the sphere to a plane, see Fig. 5. Specifically, the
constraint potential and torque are defined from the directions of
$\bm{E}_{0}$ and $\bm{\Omega}_{0}$ as follows:

$\displaystyle U^{C}$
$\displaystyle=\frac{K}{2}\left(\widehat{\bm{E}}_{0}\cdot\widehat{\bm{\Omega}}_{0}\right)^{2},$
(17) $\displaystyle\bm{N}^{C}$
$\displaystyle=K\left(\widehat{\bm{E}}_{0}\cdot\widehat{\bm{\Omega}}_{0}\right)\left(\widehat{\bm{E}}_{0}\times\widehat{\bm{\Omega}}_{0}\right).$
(18)

The vector symbol $\widehat{\cdot}$ represents a unit vector. $K>0$ is the
amplitude of the constraint torque, whose value must be sufficiently large
relative to the hydrodynamic torque $\bm{N}^{H}$ because the axis of rotation
should be fixed. The constraint torque is obtained from the potential as
$\displaystyle{\widehat{\bm{\psi}}\cdot\bm{N}^{C}=-\frac{\partial
U^{C}}{\partial\psi}=-\frac{\partial
U^{C}}{\partial\widehat{\bm{\Omega}}_{c}}\cdot\frac{\partial\widehat{\bm{\Omega}}_{0}}{\partial\psi}}$
for angle $\psi$ about any axis represented by a unit vector
$\widehat{\bm{\psi}}$, where we use the rotation formula
$\displaystyle{\frac{\partial\widehat{\bm{\Omega}}_{0}}{\partial\psi}=\widehat{\bm{\Omega}}_{0}\times\widehat{\bm{\psi}}}$
[36].

###  I.3 Simulation procedure

The simulation procedure is as follows:

  1. 1.

We obtain $\phi^{n}$ and $\bm{u}^{n}$ from the configuration
$\\{\bm{R}_{i}^{n}\\}$ and the velocities $\\{\bm{V}_{i}^{n}\\}$ and
$\\{\bm{\Omega}_{i}^{n}\\}$, where the superscript $n$ denotes the time step,
and then we calculate $\bm{u}^{*}$ using $\bm{u}^{n}$.

  2. 2.

We update the configuration $\\{\bm{R}_{i}^{n},\bm{Q}_{i}^{n}\\}$ to
$\\{\bm{R}_{i}^{n+1},\bm{Q}_{i}^{n+1}\\}$ using the velocities
$\\{\bm{V}_{i}^{n},\bm{\Omega}_{i}^{n}\\}$ and obtain the updated
$\phi^{n+1}$.

  3. 3.

We calculate the hydrodynamic force $\bm{F}_{i}^{H}$ and torque
$\bm{N}_{i}^{H}$ from the updated profile function $\phi^{n+1}$ and then
further calculate the forces $\bm{F}_{i}^{P},\bm{F}_{i}^{E}$ and torques
$\bm{N}_{i}^{P},\bm{N}_{i}^{C}$ by inter-particle interactions and external
field.

  4. 4.

From the forces and torques, we obtain the updated velocities
$\\{\bm{V}_{i}^{n+1}\\}$ and $\\{\bm{\Omega}_{i}^{n+1}\\}$; then, we return to
step 1.

The equations of motion are solved using the second-order Adams–Bashforth
method. The time step is taken to be $\Delta
t=\rho_{f}/(\eta_{f}k_{max}^{2})$, where $k_{max}=2\pi/\Delta$ is the maximum
wavenumber.

###  I.4 Numerical validations

Figure 6:  (a) Normalized roller velocity $v_{0}/(a\Omega_{0})$ as a function
of the distance between the surface of the particle and the wall, $\delta$,
for comparison with lubrication theory (far field: Goldman et al., 1967; near
field: Cichocki et al., 1998). (b) Flow field around a single roller.

We present a numerical validation for the propulsion of a colloidal roller.
The propulsion speed $v_{0}/(\Omega_{0}a)$ as a function of the distance
between the surface of the particle and the wall is derived from lubrication
theory[19, 20] as follows:

$\displaystyle f^{tr}$
$\displaystyle=\begin{cases}\frac{1}{8}\left(\frac{a}{h}\right)^{4}\left(1-\frac{3}{8}\frac{a}{h}\right)&\left(\frac{a}{h}\ll
1\right)\\\
\frac{2}{15}\ln{\left(\frac{\delta}{h}\right)}+0.25726+\frac{86}{375}\left(\frac{\delta}{a}\right)\ln{\left(\frac{\delta}{h}\right)}&\left(\frac{\delta}{a}\ll
1\right)\end{cases},$ (19) $\displaystyle f^{tt}$
$\displaystyle=\begin{cases}-\left[1-\frac{9}{16}\left(\frac{a}{h}\right)+\frac{1}{8}\left(\frac{a}{h}\right)^{3}-\frac{45}{256}\left(\frac{a}{h}\right)^{4}-\frac{1}{16}\left(\frac{a}{h}\right)^{5}\right]^{-1}&\left(\frac{a}{h}\ll
1\right)\\\
-\frac{8}{15}\ln{\left(\frac{\delta}{h}\right)}+0.95429-\frac{64}{375}\left(\frac{\delta}{a}\right)\ln{\left(\frac{\delta}{h}\right)}&\left(\frac{\delta}{a}\ll
1\right)\end{cases},$ (20)

where $h=a+\delta$. We define $F^{t}$ as the drag force acting on the particle
for translational motion and $F^{r}$ as the translational hydrodynamic force
for a rotating particle; these forces are defined as $F^{tt}=6\pi\mu
af^{tt}(\delta)v_{0}$ and $F^{tr}=6\pi\mu a^{2}f^{tr}(\delta)\Omega_{0}$.
Since the forces acting on the particle satisfy the force-free condition,
$F^{tt}+F^{tr}=0$, we obtain the relation
$v_{0}(\delta)/(\Omega_{0}a)=-f^{tr}(\delta)/f^{tt}(\delta)$. Our simulations
using the SP method can reproduce the results of lubrication theory, and the
accuracy increases as the particle radius increases; see Fig. 6(a). This
relation has been confirmed in experiments on Quincke rollers by Pradillo et
al.[11] and is a sufficiently reliable result.

The flow field around a single roller appears as shown in Fig. 6 (b). This
flow field reproduces the results of the multiblob method with a rotlet
corrected by the Blake tensor[21, 14, 22]. Here, we perform computations for a
system that accounts for the effects between the upper and lower wall
surfaces, which is slightly different from the system treated in the multiblob
method. In the multiblob method, the mobility matrix is computed via the
superposition of two-body interactions based on the Stokesian dynamics
method[13], whereas in our SP method, the hydrodynamic equations are treated
precisely by direct numerical calculation.

##  II Stabilization by the lubrication effect for a cluster of two rollers

###  II.1 Linear stability analysis for reorientation

Figure 7:  (a) Streamlines and vorticity around a single roller. (b) The
torque around a single roller caused by the lubrication effect, where the
vorticity is zero on the dotted line. (c) The reorientation of the lubrication
effect for two rollers clustered by isotropic attraction, where the green
arrows show the directions of rotation by hydrodynamic interaction when
perturbed. The direction of perturbation is the same (left side) or opposite
(right side) from each roller. The co-rotating state is unstable, while the
counter-rotating (face-to-face and back-to-back) states have configurations
that are stable against perturbations.

In this section, we discuss how a cluster of two rollers, formed by isotropic
attraction, is stabilized by the lubrication effect of the HI. According to
Faxén’s law, due to the HI between the rollers, the rollers are subjected to a
torque proportional to the vorticity
$\Omega_{f}=\frac{1}{2}\left(\bm{\nabla}\times\bm{u}\right)_{Z}$ caused by the
fluid flow[26]. The streamlines and the vorticity in the flow field on the
plane across the center of gravity of the roller are shown in Fig. 7 (a). The
vorticity is large when the relative distance from the center of gravity of
the single roller is approximately equal to the roller diameter, indicating
that the influence of the lubrication effect is large. The torque exerted by
the flow field around the single roller based on the vorticity is shown in the
schematic diagram in Fig. 7 (b). As shown in Fig. 7 (c), the fixed point lies
on the line of zero vorticity, and there are three possible reorientations:
face-to-back (co-rotating), face-to-face (counter-rotating), and back-to-back
(counter-rotating) states. When evaluating the stability of these states
against any perturbations, it can be seen that the counter-rotating states are
more stable than the co-rotating state. This stability is consistent with the
discussion by the minimum principle of viscous dissipation rate in main text.
However, we cannot know which of the face-to-face and back-to-back states is
more stable from this analysis only.

###  II.2 Forces acting on two rollers

Figure 8:  Forces due to the lubrication effect between two rollers. Blue, red
and green lines show the forces acting on the face-to-back, back-to-back and
face-to-face states, respectively. The insets show the log-log plot, where the
dotted lines have the slope $r^{-3}$, and the yellow arrows show the
directions of the hydrodynamic forces in the schematic illustrations. (a) The
relative force due to the HI, where positive and negative values represent
repulsive and attractive forces, respectively. We can see that the back-to-
back state which has a attractive forces differ from others is preferred, the
other states have repulsive forces, and the relative force acting on the face-
to-back (co-rotating) state is smaller than the forces acting on the other
states. (b) The force acting on the center of gravity due to the HI. The force
acts on same direction in the face-to-back state, accelerating the movement of
the two rollers.

To determine the stability and features of the states for two rollers under
the lubrication effect, we consider the hydrodynamic forces acting on a
cluster formed by two rollers. The forces $F^{H}_{1}$ and $F^{H}_{2}$ are
defined as the forces acting on rollers 1 (blue) and 2 (red), respectively,
and the hydrodynamic forces are divided into two components, namely, the
relative force $(F^{H}_{1}-F^{H}_{2})/2$ and the force acting on the center of
gravity $(F^{H}_{1}+F^{H}_{2})/2$; see Fig. 8. The relative force acts as an
attractive force in the back-to-back state, and this state is preferred under
the lubrication effect. In contrast, the relative hydrodynamic force acts as a
repulsive force in the face-to-back or face-to-face states, and this states is
not preferred under the lubrication effect, as shown in Fig. 8(a). Only in the
face-to-back (co-rotating) state the force act in the direction of propulsion,
causing the two rollers to speed up; see Fig. 8(b). This acceleration
phenomenon has been observed experimentally, and it is known that the
acceleration effect increases with higher density and higher rotation
speed[3]. Thus, the lubrication effect also improves the mobility.

##  III Analysis of many-roller system

###  III.1 Order parameters

Figure 9:  Fingerprint of each state created by the order parameters. The
polar order parameter (POP), the bond-orientational order parameter (BOP), the
rotational order parameter (ROP) and the chain order parameter (COP) are used
to classify the states. The letters ‘H’ and ‘L’ denote high and low values of
the order parameters, respectively, and the threshold values are defined as
$0.6$ (POP), $0.40$ (BOP), $0.25$ (ROP) and $0.6$ (COP).

To characterize each state, we define a polar order parameter (POP), a bond-
orientational order parameter (BOP), a rotational order parameter (ROP) and a
chain order parameter (COP)[28], and we show the corresponding fingerprints of
each state (see Fig. 9). Note that the POP is defined as the average of the
propulsion directions of all particles, $\left|\braket{\bm{n}_{i}}\right|$,
where $\bm{n}_{i}$ shows the unit vector of rotation axis by Quincke rotation
of the $i$-th roller and $\braket{\cdots}$ denotes an average over the
particles and time, while the BOP and COP reflect whether the particles are
arranged in a hexagonal configuration or a straight line. The BOP is defined
as $\braket{}{\psi_{i}}{}$,

$\displaystyle\psi_{i}=\sum_{j=1}^{Z_{i}}\exp{(i6\theta_{ij})}/Z_{i},$ (21)

where $Z_{i}$ is the coordination number of the $i$-th particle obtained from
a Voronoi construction for the particle configuration and $\theta_{ij}$ is the
angle between a reference axis and the bond between the $i$-th particle and
its $j$-th neighbor. $\psi=1$ indicates perfect hexagonal ordering, whereas a
completely disordered structure corresponds to $\psi=0$. The ROP is defined as
$\braket{}{\braket{\hat{\bm{\Omega}}_{0,i}\cdot\hat{{\bm{r}}^{\prime}}_{i}}_{c}}{}$,
where $\braket{\cdots}_{c}$ denotes taking the average for the particles that
make up each cluster in each snapshot and the operation $\braket{\cdots}$
denotes taking the average over the values for each cluster.
${\bm{r}}^{\prime}_{i}=\bm{r}_{i}-\bm{r}_{G}$, with $\bm{r}_{G}$ being the
center of gravity for each cluster, where a cluster is defined as a set of
rollers that are in contact with each other. The COP is defined as
$\braket{c_{i}}$,

$\displaystyle
c_{i}=\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}_{ik}}{r_{ij}r_{ik}}\right)\right],$
(22)

where $M$ represents the 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, and $\mathcal{S}_{i}$ denotes the region where other
particles are in contact with the $i$-th particle, i.e., a circle with a
radius on the order of the particle diameter centered on the $i$-th
particle[28].

###  III.2 Two-body correlation function for the angular velocity vector

Figure 10:  The two-body angular velocity correlations $C_{\omega\omega}(r)$
for each state; (a) polar, disordered, and rotating cluster states. (b) gas,
polar liquid, vortex, and chain states.

Furthermore, to provide an understanding of the characteristic features of
each behavior, we introduce the two-body angular velocity correlations as
follows:

$\displaystyle
C_{\omega\omega}(r)=\frac{\Braket{\bm{n}_{i}\cdot\bm{n}_{j}\delta(r-r_{ij})}}{2\pi
r\Delta r\phi_{0}(r)},$ (23)

where $\phi_{0}(r)$ is the local area fraction of the rollers and $2\pi
r\Delta r\phi_{0}(r)$ is the number of rollers in a ring-shaped parcel of
width $\Delta r$ and radius $r$. This function takes values of
$C_{\omega\omega}\to 1$ in the parallel state and $C_{\omega\omega}\to-1$ in
the antiparallel state. Fig. 10(a) shows functions related to cluster motion,
which are characterized by larger fluctuations than those observed in Fig.
10(b) for motion under high applied electric fields. In particular, the
rotating cluster function is characterized by a negative peak around the
particle diameter and oscillating between positive and negative values, as in
Fig. 10(a). The function for chain state forms peaks at intervals roughly
corresponding to the diameter scale $2a$, as in Fig. 10(b), while the function
in the vortex state shows a maximum directional order at the diameter and then
decays as a longer distance between rollers. In the vortex state, the
directional order is maximal at the diameter and then decays with distance. It
can be seen that the disorder of the disordered cluster and gas states and the
order of the polar clusters and liquids are manifested as the correlations in
the long-distance.

###  III.3 The breaking of the polar state for denser region

Figure 11:  The POP depends on the density of the rollers at $E=3.0$. The
average distance occupied each roller is defined as $l/a\equiv
1/\sqrt{\phi_{0}}$, and the shaded area represents the region inside the
roller, $l/a\leq 1$, where $\sigma$ is the diameter of the particle. The
dotted line represents the experimental findings of Bricard et al., 2013 at
$E=1.39$, and it is seen that our simulation results reproduce the experiments
in the lower-density region. On the other hand, in the higher-density region,
the ordering of the flocking is rapidly broken due to the lubrication effect
of the HI as the average radius of the region occupied by a single particle
becomes smaller than $l/a\sim 2$.

Here, we observe that the POP depends on the density of the rollers, and our
results reproduce the experiments[1] in the lower-density region. We consider
the characteristic distance between particles at which the polar state breaks
down at the higher-density region, i.e., the characteristic distance at which
the dipole–dipole interaction and the HI lubrication effect switch in
dominance. The average area occupied by each roller is $\bar{S}\equiv
L^{2}/N=\pi a^{2}/\phi_{0}$, and we assume that the shape of this region is
circular and define the characteristic distance as
$l\equiv\sqrt{\bar{S}/\pi}=a/\sqrt{\phi_{0}}$. As $l/a$ decays with increasing
density, the POP rapidly becomes disordered below $l/a\sim 2$, as shown in
Fig. 11. This characteristic distance, which is equivalent to the diameter
scale $2a$, corresponds to the range where the lubrication effect has a large
influence on the interparticle interaction and is considered to represent the
distance at which the dipole–dipole interaction and lubrication effects switch
in dominance in a many-roller system. As the result, the Lubrication effect of
HI is thought to promote further disordering at the higher-density region for
a many-roller system.

##  IV Descriptions of Supplementary Animations

We provide animations capturing the time evolution of each state, as obtained
from simulations.

Supplementary Movie 1

The motions of a two-roller can be classified into two types: co-rotating
($\theta=0$) and counter-rotating ($\theta=-\pi,\pi$) states. In the co-
rotating state, the cluster moves in one direction at a faster speed than that
of a single roller. On the other hand, in a counter-rotating state, the
cluster is immobile because each roller moves in the opposite direction of the
other. In these animations, the color bar represents the normalized amplitude
of the velocity field, and the arrows represent the direction of the velocity
field at each point.

Supplementary Movie 2

When the electric dipole moment strength is small, the axes of rotation are
disordered by the lubrication effect. As a result, the cluster is immobile. On
the other hand, when the dipole strength is large, the cluster starts to
rotate due to the stabilization of the rotational axes by the dipole–dipole
interaction. The rollers are colored in accordance with the direction of the
rotational axis, and the arrows also show the direction.

Supplementary Movie 3

This diagram was obtained from simulations run with parameters corresponding
to the theory in the steady state. The colors of the frames represent the
types of states; blue indicates polar liquid states, gray indicates gas
states, green indicates vortex states, purple indicates chain states, orange
indicates rotating cluster states, olive indicates disordered cluster states,
and azure indicates polar cluster states.

Supplementary Movie 4

A gas state, in which each roller is moving randomly. The gas states are
thought to be caused by an imbalance between HI and dipole interactions
between the rollers. The long-range HI is dominant at low densities, while the
short-range dipole interaction is dominant at high densities. Parameters:
$E=3.0$ and $\phi_{0}=0.01,0.2$.

A polar liquid state, in which many rollers move in one direction and their
distances from each other remain almost constant during the motion. Isotropic
repulsion causes the rollers to separate from each other because of the
dominant interaction of $\mu_{\perp}$. The rollers are not scattered, but are
regularly arranged due to the periodic boundary condition. At the right
density, the distance between the rollers is just enough to balance the
interaction between the horizontal components of HI and dipole, which is
thought to produce directional order. Parameters: $E=3.0$ and $\phi_{0}=0.08$.

A vortex state, in which all the rollers generate enough vortex to cover the
entire system. The average distance between the rollers is almost of the order
of the diameter, and not only the electric dipole interaction but also the
lubrication effect is strongly acting. As a result, the torque due to the
lubrication effect gradually bends the propulsive direction of the neighboring
particles, which is thought to create a vortex. Parameters: $E=3.0$ and
$\phi_{0}=0.8$.

A chain state, in which the front and back of the rollers are attached to each
other. Because the isotropic attraction due to the electrohydrodynamic flow is
balanced by the isotropic repulsion by $\mu_{\perp}$, $\mu_{\parallel}$ is
dominant, creating a chain. Since the parallel arrangement of rollers like a
chain is one of the fixed points in the lubrication effect, the chain state is
maintained as long as there is no perturbation due to motion more intense.
Parameters: $E=2.0$ and $\phi_{0}=0.08$.

A rotating cluster state, in which each cluster rotates clockwise or
counterclockwise. The rotational motion of the cluster is caused by the
stabilization of each roller’s rotational axis by the dipole–dipole
interaction. Parameters: $E=1.25$ and $\phi_{0}=0.08$.

A disorder cluster state, in which there is no order regarding the direction
of propulsion and rotation for the entire cluster. There are two possible
mechanisms for the generation of these disordered clusters. First, when the
dipole is small and the HI is dominant, the cluster becomes disordered due to
the lubrication effect. On the other hand, as higher density, the proportion
of clusters that have anistropic shape when the dipole is large increases, and
the rotational order within the clusters due to the dipole is decreased. As a
result, each cluster makes random translational and rotational motions, which
repeatedly coalesce and separate, creating a disordered state as a whole.
Parameters: $E=1.25$ and $\phi_{0}=0.5$.

A polar cluster state, in which all the rollers are moving in one direction.
As higher density sufficiently, all the rollers come together to form a large
cluster that spans the boundaries of the system. Because of the anisotropic
shape of this cluster, the rotational order is small and the motion is almost
translational. Parameters: $E=1.25$ and $\phi_{0}=0.8$.

## References

  * [1] Bricard, A., Caussin, J.-B., Desreumaux, N., Dauchot, O. & Bartolo, D.  Emergence of macroscopic directed motion in populations of motile colloids.  _Nature_ 503, 95–98 (2013).
  * [2] Bricard, A. _et al._ Emergent vortices in populations of colloidal rollers.  _Nature communications_ 6, 1–8 (2015).
  * [3] Lu, S. Q., Zhang, B. Y., Zhang, Z. C., Shi, Y. & Zhang, T. H.  Pair aligning improved motility of quincke rollers.  _Soft Matter_ 14, 5092–5097 (2018).
  * [4] Karani, H., Pradillo, G. E. & Vlahovska, P. M.  Tuning the random walk of active colloids: From individual run-and-tumble to dynamic clustering.  _Physical review letters_ 123, 208002 (2019).
  * [5] Geyer, D., Martin, D., Tailleur, J. & Bartolo, D.  Freezing a flock: Motility-induced phase separation in polar active liquids.  _Physical Review X_ 9, 031043 (2019).
  * [6] Mauleon-Amieva, A. _et al._ Competing active and passive interactions drive amoebalike crystallites and ordered bands in active colloids.  _Physical Review E_ 102, 032609 (2020).
  * [7] Zhang, B., Karani, H., Vlahovska, P. M. & Snezhko, A.  Persistence length regulates emergent dynamics in active roller ensembles.  _Soft Matter_ 17, 4818–4825 (2021).
  * [8] Liu, Z. T. _et al._ Activity waves and freestanding vortices in populations of subcritical quincke rollers.  _Proceedings of the National Academy of Sciences_ 118 (2021).
  * [9] Quincke, G.  Ueber rotationen im constanten electrischen felde.  _Annalen der Physik_ 295, 417–486 (1896).
  * [10] Pannacci, N., Lobry, L. & Lemaire, E.  How insulating particles increase the conductivity of a suspension.  _Physical review letters_ 99, 094503 (2007).
  * [11] Pradillo, G. E., Karani, H. & Vlahovska, P. M.  Quincke rotor dynamics in confinement: rolling and hovering.  _Soft Matter_ 15, 6564–6570 (2019).
  * [12] Swan, J. W. & Brady, J. F.  Simulation of hydrodynamically interacting particles near a no-slip boundary.  _Physics of Fluids_ 19, 113306 (2007).
  * [13] Balboa Usabiaga, F. _et al._ Hydrodynamics of suspensions of passive and active rigid particles: a rigid multiblob approach.  _Communications in Applied Mathematics and Computational Science_ 11, 217–296 (2017).
  * [14] Driscoll, M. _et al._ Unstable fronts and motile structures formed by microrollers.  _Nature Physics_ 13, 375–379 (2017).
  * [15] Shaebani, M. R., Wysocki, A., Winkler, R. G., Gompper, G. & Rieger, H.  Computational models for active matter.  _Nature Reviews Physics_ 2, 181–199 (2020).
  * [16] Yamamoto, R., Molina, J. J. & Nakayama, Y.  Smoothed profile method for direct numerical simulations of hydrodynamically interacting particles.  _Soft Matter_ 17, 4226–4253 (2021).
  * [17] Nakayama, Y. & Yamamoto, R.  Simulation method to resolve hydrodynamic interactions in colloidal dispersions.  _Physical Review E_ 71, 036707 (2005).
  * [18] Nakayama, Y., Kim, K. & Yamamoto, R.  Simulating (electro) hydrodynamic effects in colloidal dispersions: smoothed profile method.  _The European Physical Journal E_ 26, 361–368 (2008).
  * [19] Goldman, A. J., Cox, R. G. & Brenner, H.  Slow viscous motion of a sphere parallel to a plane wall―i motion through a quiescent fluid.  _Chemical engineering science_ 22, 637–651 (1967).
  * [20] Cichocki, B. & Jones, R.  Image representation of a spherical particle near a hard wall.  _Physica A: Statistical Mechanics and its Applications_ 258, 273–302 (1998).
  * [21] Blake, J. & Chwang, A.  Fundamental singularities of viscous flow.  _Journal of Engineering Mathematics_ 8, 23–29 (1974).
  * [22] Delmotte, B.  Hydrodynamically bound states of a pair of microrollers: A dynamical system insight.  _Physical Review Fluids_ 4, 044302 (2019).
  * [23] Strutt, J. W.  Some general theorems relating to vibrations.  _Proceedings of the London Mathematical Society_ 1, 357–368 (1871).
  * [24] Onsager, L.  Reciprocal relations in irreversible processes. i.  _Physical review_ 37, 405 (1931).
  * [25] Onsager, L.  Reciprocal relations in irreversible processes. ii.  _Physical review_ 38, 2265 (1931).
  * [26] Lauga, E.  _The fluid dynamics of cell motility_ , vol. 62 (Cambridge University Press, 2020).
  * [27] Das, D. & Saintillan, D.  Electrohydrodynamic interaction of spherical particles under quincke rotation.  _Physical Review E_ 87, 043014 (2013).
  * [28] Imamura, S. & Kawakatsu, T.  Modeling of chemically active particles at an air–liquid interface.  _The European Physical Journal E_ 44, 1–12 (2021).
  * [29] López Llobet, L. & Fernández-Nieves, A.  The important role of electrohydrodynamics in a system of quincke rollers.  _Bachelor theses, Universitat de Barcelona_ (2020).  URL http://hdl.handle.net/2445/171367.
  * [30] Liao, G.-J., Hall, C. K. & Klapp, S. H.  Dynamical self-assembly of dipolar active brownian particles in two dimensions.  _Soft Matter_ 16, 2208–2223 (2020).
  * [31] Kaiser, A., Snezhko, A. & Aranson, I. S.  Flocking ferromagnetic colloids.  _Science advances_ 3, e1601469 (2017).
  * [32] Allen, M. P. & Germano, G.  Expressions for forces and torques in molecular simulations using rigid bodies.  _Molecular Physics_ 104, 3225–3235 (2006).
  * [33] Aguado, A. & Madden, P. A.  Ewald summation of electrostatic multipole interactions up to the quadrupolar level.  _The Journal of chemical physics_ 119, 7471–7483 (2003).
  * [34] Allen, M. P. & Tildesley, D. J.  _Computer simulation of liquids_ (Oxford university press, 2017).
  * [35] Dwelle, K. A. & Willard, A. P.  Constant potential, electrochemically active boundary conditions for electrochemical simulation.  _The Journal of Physical Chemistry C_ 123, 24095–24103 (2019).
  * [36] Goldstein, H., Poole, C. & Safko, J.  _Classical mechanics_ (American Association of Physics Teachers, 2002).
  * [37] Schmidt, J. _et al._ Conductivity in nonpolar media: Experimental and numerical studies on sodium aot–hexadecane, lecithin–hexadecane and aluminum (iii)-3, 5-diisopropyl salicylate–hexadecane systems.  _Journal of colloid and interface science_ 386, 240–251 (2012).