Active colloid with externally induced periodic bipolar motility
and its cooperative motion
Airi N. Kato1,2 ,∗ Kazumasa A. Takeuchi1 ,† and Masaki Sano3‡
1

arXiv:2203.12882v3 [cond-mat.soft] 14 Jul 2022

Department of Physics, The University of Tokyo,
7-3-1 Hongo, Bunkyo-ku,Tokyo 113-0033, Japan.
2
Wenzhou Institute, University of Chinese Academy of Sciences, Wenzhou, Zhejiang 325001, China and3
Institute of Natural Sciences, School of Physics and Astronomy,
Shanghai Jiao Tong University, Shanghai 200240, China.
(Dated: July 15, 2022)
Active matter physics has been developed with various types of self-propelled particles, including
those with polar and bipolar motility and beyond. However, the bipolar motions experimentallyrealized so far have been either random along the axis or periodic at intrinsic frequencies. Here we
report another kind of bipolar active particles, whose periodic bipolar self-propulsion is set externallyat a controllable frequency. We used Quincke rollers—dielectric particles suspended in a conducting
liquid driven by an electric field—under an AC electric field instead of the usually used DC field.
Reciprocating motion of a single particle at the external frequency was observed experimentallyand characterized theoretically as stable periodic motion. Experimentally, we observed not only
the reciprocating motion but also non-trivial active Brownian particle (ABP)-like persistent motion
in a long time scale. This resulted in a Lorentzian spectrum around zero frequency, which is notaccounted for by a simple extension of the conventional model of Quincke rollers to the AC field. It
was found that ABP-like motion can be reproduced by considering the top-bottom asymmetry in
the experimental system. Moreover, we found a rotational diffusion coefficient much larger than thethermal one, as also reported in previous experiments, which may have resulted from roughness ofthe electrode surface. We also found self-organized formation of small clusters, such as doublets andtriplets, and characterized cooperative motion of particles therein. The AC Quincke rollers reported
here may serve as a model experimental system of bipolar active matter, which appears to deserve
further investigations.

1.

INTRODUCTION

Active matter is a class of intrinsically non-equilibrium systems, usually consisting of a collection of self-propelled
particles (SPPs) that move according to the individual polarity or axis of motility [1–4]. Characteristic featuresof active matter that do not appear in equilibrium systems, such as collective motions (e.g., ordering [1, 4, 5] anddynamic clustering [2, 6–9]) and giant number fluctuations [1, 4, 10–12], have been a target of intensive studies for
decades. These include experimental studies using systems from microscopic to macroscopic scales, and from biologicalto artificial systems. Examples of artificial experimental systems are shaken rods and disks [12, 13], Janus particles[6–9, 14–20], Quincke rollers [21–24], floating droplets [25, 26], and so on (see [2] for a review).
One of the main advantages of using artificial self-propelled particles is controllability. Typically, self-propulsion
speed can be controlled externally within some range; for example, speed of Janus particles driven by induced-charge
electrophoresis [16, 18–20] and that of Quincke rollers [21, 22] can be controlled by the strength of the electric field.
In such systems driven by an external field, active and passive states can also be switched by turning the field on
and off [23, 27, 28] or otherwise [19], leading to motions reminiscent of run-and-tumbling of E.coli. [29]. In thesecases, the motion is unidirectional, i.e., particles are self-propelled in the direction of the polarity while the field is
on. In contrast, some systems such as shaken rods [12] exhibit bipolar motion by stochastic reversal of self-propelling
direction. Together with biological examples such as Myxococcus xanthus [30] and neural progenitor cells [31], thesesystems have been studied as examples of active nematics [3]. There are also examples of bipolar active particles withperiodic reversal of self-propelling direction, such as self-propelled Belousov-Zhabotinsky droplets [32], a surface-wave
driven droplet [26], and Quincke rollers under a strong DC electric field [24]. However, in all these examples, theperiodic reversal takes place at a frequency intrinsic to each system, which is therefore uncontrollable.
Here we report that Quincke rollers, driven by an AC electric field instead of the usually used DC field, constitute anexperimental system of bipolar active particles with periodic reversal of self-propelling direction at an externally set
frequency (Fig. 1). We call our system “AC Quincke rollers” afterwards to distinguish from the original DC Quincke∗ airi@ucas.ac.cn† kat@kaztake.org‡ sano.masaki@sjtu.edu.cn

2

FIG. 1. Schematics of a Quincke roller under a DC field (a) and an AC field (b). (a) Under a DC field, a constant polarization
is formed and inclined in steady state, leading to a constant self-propelling velocity. (b) Under an AC field, the polarizationand the velocity become time-dependent. We found periodic bipolar motility at the frequency of the external field in this case.

rollers. A DC Quincke roller consists of a dielectric sphere immersed in a conducting liquid under a constant DC
electric field E0 . The particle rolls on the electrode by the Quincke effect [33] as shown in Fig. 1a. The Quincke effectoccurs under the condition that the relaxation time of ions in liquid, τl := l /σl , is shorter than that on the surfaceof the particle, τp := p /σp . Here, l,p is the dielectric constant of the liquid and the particle, respectively, and σl,pis the electric conductivity of the liquid and the particle, respectively. When the electric field E0 is applied, charges
accumulate on the surface of the particle, which makes an effective polarization P antiparallel to the field (Fig. 1a).
This configuration is unstable because a small perturbation of the polarization perpendicular to the field is amplified
by the electric torque P × E0 . Therefore, if this torque surpasses the viscous restitution torque, the particle starts
to rotate spontaneously. As a result, the particle rolls on the substrate electrode (Fig. 1a), leading to unidirectional
self-propelling motion in the two-dimensional plane. The direction of motion is selected by spontaneous symmetry
breaking. Note that the particle itself has no polarity without the external electric field. It is the polarization P that
determines the direction of motion (polarity), which exists only when the external field is applied. The fact that thepolarity is acquired only a posteriori under the external field is the key to realize the externally induced periodic
bipolar motion realized by our AC Quincke rollers.
In this paper, we focus on the behavior of the AC Quincke rollers in a dilute suspension. First we describe theexperimental setup (Sec. 2) and show experimental results on the motion of a single particle (Sec. 3). We foundthat a single particle exhibits periodic reciprocating motion with active Brownian particle (ABP)-like persistent net
motion. Next, the theoretical description of a Quincke roller under a DC field was generalized to that under an AC
field (Sec. 4). It was found that the reciprocating motion at the frequency of the external field corresponds to alimit cycle of the equations of motion. This supports our observation that AC Quincke rollers are a realization ofcontrollable SPPs with stable reciprocating motions at the frequency of the external field. The ABP-like motion was
reproduced by incorporating the asymmetry caused by the existence of the lower electrode. Then we focus on motionof two- and three-particle clusters (doublet and triplet) which were spontaneously formed and maintained for long
periods of time (at least 100–1000 periods) (Sec. 5). Characteristic correlations in particle motions are revealed andinterpreted in terms of particle interactions. We believe that the AC Quincke rollers reported here may provide a
new perspective on active matter systems with periodic bipolar self-propulsion, through nontrivial consequences thatsuch time dependence brings.

2.

EXPERIMENTAL SETUP AND METHODS

Quincke rollers under an AC electric field were experimentally realized as follows. Following Bricard et al. [21, 22] on
the DC Quincker rollers, we used polystyrene spherical colloids (Thermo Scientific GS0500, radius 2.5 µm) dispersed in
0.15 mol/L of AOT/hexadecane solution. As shown in Fig. 2, the suspension was sandwiched between two transparent
electrodes (glass plate coated with indium-tin oxide (ITO) film, Mitsuru Optical Co. Ltd.). The electrodes were spaced
by 30 µm thick double-sided polyethylene terephthalate tapes (NITTO, No.5603). The AC electric field was applied
by a function generator (Agilent, 33220A) through an amplifier (NF, TA-120). [34] The particles were observedby bright-field microscopy using an inverted microscope (Olympus, IX70) with an x40 objective lens (Olympus,

3

FIG. 2. Experimental setup. A suspension of polystyrene spheres including surfactant (AOT) is sandwiched between two
transparent electrodes. The particles roll on the lower substrate under an electric field, moving thereby on the quasi-twodimensional plane.

LUCPLFLN40X, NA 0.60) unless otherwise stipulated, and captured by a high-speed camera (Photron, Fastcam miniAX) at 3000 fps for 300 Hz and 2500 fps for 250 Hz. The image acquisition was synchronized by the trigger signal
from the function generator. Tracking was conducted using Particle Tracker 2D/3D [35], a plugin of ImageJ. When an
AC electric field was applied, we observed approximately periodic motion after a short transient time. We therefore
conducted all measurements after waiting long enough, specifically at least 5 seconds after the start of the voltage
application. In this paper, unless otherwise stipulated, we fixed the root-mean-square amplitude and the frequencyof the AC applied voltage at 150 V and 300 Hz or 250 Hz, respectively.
3.

EXPERIMENTAL RESULTS ON SINGLE AC QUINCKE ROLLERS

This section reports experimental results on single particle motion. Here we used isolated particles, or more
specifically, particles separated more than 100 µm from others. Typical trajectories of isolated particles at 300 Hz are
shown in Fig. 3a. In short time scales, the observed particles exhibited reciprocating motion at the frequency of theexternal field, with some extent of randomness in the direction and amplitude of the reciprocating motion (Fig. 3a,
inset). By contrast, in time scales much longer than the reciprocation period, we observed persistent motion whose
direction changed gradually (see Videos S1 and S2 in ESI†). Using each component of the single particle velocityv = (vx , vy ), we computed the energy spectral density
Svv (f ) := |v̂(f )|2 ,v̂(f ) :=

N −1

1 X
v(tk )e−2πif tk ,
N

(1)

k=0

with v(t) = vx (t) or vy (t), where tk is the time at the kth video frame and N is the total frame number. The results
are shown in Fig. 3b. To reduce statistical fluctuations in a low-frequency region, we computed the power spectra in
the inset of Fig. 3b by segmenting the single whole time series into 30 segments and then averaging them. This showsthat the velocities had the periodicity at the external frequency, 300 Hz. The spectra also showed another signal in a
low-frequency region (. 10 Hz) as shown in the inset of Fig. 3b, which corresponds to the slowly changing, persistentmotion we observed. This low-frequency mode was also observed for 200Hz and 250Hz (Fig. S2, ESI†). In fact, thespectral density in the low-frequency region is Lorentzian, which is expected from a simple one-particle model with asinusoidal and a nonzero constant self-propelling velocity v(t) = v0 + v1 cos ωt (see Supplemental Notes, ESI†), andalso seen for the ABP [17, 36]. More precisely, the spectral density S(f ) for the model with v(t) = v0 + v1 cos ωt is


2Dθ
Dθ
v12
Dθ
S(f ) = v02 2
+
+
,
(2)
Dθ + (2πf )2
2 Dθ2 + (ω + 2πf )2Dθ2 + (ω − 2πf )2
where Dθ is the rotational diffusion coefficient. This coefficient was estimated at Dθ = 31.4 ± 0.6 s−1 (Dθ−1 ∼ 0.03 s)
by the Lorentzian (the first term of Eq. (2)) fitting to the average of the four spectra in the inset of Fig. 3b for< 40 Hz; here the fitting was done with the peak value of the Lorentzian fixed at S(0). The velocities v0 and v1
were estimated at v0 = 7.3 ± 0.4 µm/s and v1 = 24 ± 8 µm/s by the peaks at 0 and 300 Hz. Note that this effective
rotational diffusion coefficient Dθ is four orders of magnitude larger than the thermal rotational diffusion coefficient4

FIG. 3. Single-particle behavior (300 Hz, 150 V). (a) Two independent trajectories (over 1091 periods = 3.639 s). The subpanels
show the positions of the particles recorded every half period for the first 10 periods. The initial positions are depicted by theblack dots. The two particles were captured simultaneously, with the distance varying in the range of 177–198 µm. (b) Energy
spectral densities. The particles 1 and 2 in the legend correspond to the red and blue trajectories in (a), respectively. The inset
shows an enlargement of a low-frequency region. To reduce statistical fluctuations at a low-frequency region, the single wholetime series was divided into 30 segments (each segment contains 36 periods), then the power spectra were averaged. The yellow
2Dθ
solid line shows the result of the fitting to the average of the four data by the Lorentzian spectrum v02 D2 +(2πfin Eq. (2). The)2
θ

data for f < 40 Hz was used for fitting. (c) The MSDs of three isolated particles (red, blue, and green dots) and their average(black solid line). The dashed lines are guides for eyes with the indicated slopes.

Dθ = kB T /(8πηa3 ) ∼ 10−3 s−1 , similarly to the case of Bricard et al. [22]. They reported Dθ−1 ∼ 0.31 s in the DC case.
In contrast to these Quincke experiments, such a considerable deviation between the measured rotational diffusion
coefficient of a spherical particle and the thermal value was not observed in experiments of thermally activated Janus
particles[17], which did not use an electrode. Therefore, this athermal rotational diffusion may be originated fromproperties of the electrode, such as its surface roughness. To inspect this possibility, we measured surface roughness of
the ITO film by using an atomic force microscope (AFM) and discussed possible effects of quenched noise originatedfrom the surface roughness in the Supplemental Notes (ESI†). The origin of the DC component of the velocity will
be discussed in Sec. 4.
To characterize the dynamics of the AC Quincke rollers, we also computed their mean square displacements (MSDs):

MSD(τ ) =

N −m

X
1
(x(tk+m ) − x(tk ))2 ,
N −m

(3)

k=1

where τ = (m/3000) s is the lag time, and x(tk ) is the particle position at time tk . The results for three isolatedparticles are shown in Fig. 3c. The MSDs oscillate with the period T = (1/300) s of the applied field and growalgebraically over a duration longer than 20T with a superdiffusive exponent 1.70 ± 0.07 as shown by the black dashedline in Fig. 3c. This estimate was obtained by using the data points taken at every 10 frames (1 period) in the range
T ≤ τ ≤ 20T , and the 95% confidence interval was given. We consider that this superdiffusion is associated with the
low-frequency mode that manifested itself in the long-time persistent motion. After this superdiffusive behavior, theMSD exponent crossovers to the diffusive one. This crossover from the nearly ballistic behavior to the diffusive one isanother aspect similar to the ABP, in addition to the aforementioned Lorentzian around 0 Hz in the spectral densities
characterized by the effective rotational diffusion.

5

FIG. 4. Single-particle behavior in the theoretical model (300 Hz, 150 V). (a) Time evolution of the non-dimensional polarizat= 82.1 µm/s. (b) Sketch of the reciprocatingtion P̃kσ (red) and P̃zσ (blue), and the particle velocity v (green). We set µraµτMWrolling mechanism of an AC Quincke roller. Note that the sign of P̃kσ is kept the same. (c) A trajectory of (P̃kσ , P̃zσ ) from atransient to the limit cycle. The red line shows the periodic cycle over one period. The color gradient shows the time change.
(d) The autospectra of the polarization (top) and the velocity (bottom) from the data in (a) in the periodic state. (e) The
autospectra of the polarization (top) and the velocity (bottom) in the biased model (α = 1). (f) The positions of a particle
x(t) with no bias (α = 0) and with a bias (α = 1).

4.

THEORY ON SINGLE AC QUINCKE ROLLERS

To gain insights on the single particle dynamics observed in the previous section, here we provide a theoreticaldescription of single AC Quincke rollers, by extending the theory for DC Quincke rollers formulated by Bricard etal. [21]. In the DC case, the locomotion of a single Quincke roller of radius a under an electric field E(t) = E0 is
described by its angular momentum Ω (t), which obeys the equation of motion
I

Ω
dΩ
l
= (P × E0 ) − µ−1
r Ω
dt
0

(4)

with the rotational inertia I, the dielectric constant of vacuum 0 , and the rotational mobility µr expected to beaffected by lubrication [21, 37–40].
The polarization P(t) evolves by the following electrodynamic equation [21, 41, 42]:12π0 a3dP+
P=−
E0 + Ω × (P − 4π0 a3 χ∞ E0 ),
dt
τMW
τMW
 −

(5)

 +2

with χ∞ = pp+2ll and the Maxwell-Wagner time τMW = p2σl l . Since the inertia term is negligible (µr I ∼ 10−8 s τMW in our system), we have Ω(t) = µr 0l (P × E0 ). At this point, it is convenient [21] to recognize that P(t) consistsof the following two contributions: P(t) = P (t)+Pσ (t), where P (t) = 4π0 a3 χ∞ E0 is the dielectric contribution due
to the permittivity discontinuity at the interface, and the rest Pσ (t) is due to the surface charge transport. The chargetransport contribution Pσ (t) can be further decomposed into the z-component Pzσ (t) and the component parallel tothe surface, Pσk (t). With this and Ω(t) = µr 0l (P × E0 ), the time evolution equations for Pzσ (t) and Pσk (t) read


dPzσ
1
4π0 a31
l+Pzσ = −
χ∞ +E0 + µr E0 (Pkσ )2 ,
dt
τMW
τMW
2
0
dPσk
1
l
+
Pσ = −µr E0 Pzσ Pσk .
dt
τMW k
0

(6)
(7)

For a strong enough field E0 , this set of the equations has a stable fixed point corresponding to the motion at a6constant translation velocity. It is given byv=−

laµt E0 Pσk ,0(8)with the cross mobility coefficient µt that relates the electric torque to the translational velocity [21].
Now we generalize this formulation to deal with an oscillating electric field, E(t) = E0 cos ωt. Replacing E0
in Eqs. (4) and (5) with E0 cos ωt and using P(t) = P (t) + Pσ (t) with P (t) = 4π0 a3 χ∞ E0 cos ωt, under theassumption that dielectric polarization takes place fast enough compared to T = 2π/ω (typical dielectric polarization
response is considered to be instantaneous, much faster than our external frequency: 300 Hz), we obtain the following
time evolution equations:

 ∞
dPzσ
1
χ + 1/2σ3∞+
P = −4π0 a E0
cos ωt − χ ω sin ωt
dt
τMW z
τMW
l
(9)
+ µr E0 (Pkσ )2 cos ωt,
0
dPσk
1
l
+
Pσk = −µr E0 Pzσ Pσk cos ωt.
(10)
dt
τMW
0
The rolling velocity of the particle isv(t) = −laµt E0 Pσk (t) cos ωt.
0

(11)

Now, we introduce the following non-dimensionalized quantities: t̃ := t/τMW , ω̃ := ωτMW , P̃ := ε0 EP0 a3 , and A :=µr εl E02 a3 τMW . The rolling velocity, along the axis of motion isv(t) = −aµt
AP̃kσ cos ω̃ t̃.
µr τMW(12)

Equations (9) and (10) are then rewritten as follows:dP̃zσ+ P̃zσ = −4π (χ∞ + 1/2) cos ω̃ t̃ − χ∞ ω̃ sin ω̃ t̃
dt̃
+ A(P̃kσ )2 cos ω̃ t̃,

dP̃σk
dt̃

+ P̃σk = −AP̃zσ P̃σk cos ω̃ t̃.

(13)
(14)

ω
= 300 Hz and τMW = 1 ms. Note that the direction of P̃σk (t)We numerically solved Eqs. (13) and (14) with 2π
in Eq. (14) is determined by the spontaneous symmetry breaking and changes by the effective rotational diffusion.
Figure 4a displays the time evolution of P̃kσ (t), P̃zσ (t), and v(t), from the initial condition (P̃kσ (0), P̃zσ (0)) = (0.1, 0.1).

The polarization (P̃kσ , P̃zσ ) reaches the periodic state, in which P̃zσ (t) oscillates with the period identical to that of

the external field, but P̃kσ (t) does so with the half period (fundamental frequency 2ω) without changing its sign. As
a result, the velocity v(t) oscillates at the period of the external field, corresponding to the reciprocating motion as
observed experimentally. This rolling mechanism of an AC Quincke roller is depicted in Fig. 4b. Physically, when theexternal field is applied, as in the DC case, the component P̃kσ is initially developed by instability due to ionic flow
around the particle and reaches a characteristic strength. In the AC case, the field is reversed every half period andthe ionic flow is also reversed accordingly in the z-direction, but the charge distribution remains biased in such a waythat the direction (sign) of P̃kσ is conserved. Figure 4c shows the trajectory in (P̃kσ , P̃zσ ) space, which goes to the redclosed loop (limit cycle) after a few cycles. Trajectories from different initial conditions went to the same limit cycle
as long as the sign of P̃kσ (0) was positive. If P̃kσ (0) < 0, the limit cycle to realize was the one obtained as a mirror
image (with respect to the P̃z axis) of the limit cycle shown in Fig. 4c. This is because, if (P̃kσ (t), P̃zσ (t)) is a solution

to Eqs. (13) and (14), (−P̃kσ (t), P̃zσ (t)) is also a solution. In our model, this relaxation to the limit cycle and theresulting reciprocating motion were observed in a broad range of the frequency, from the low frequency limit f → 0
up to f ∼ 1/τMW . For higher frequencies, the limit cycle does not result in reciprocating motion because P̃kσ (t) = 0(see Fig. S1, ESI†).

7
To compare the theory and the experimental result further, the autospectral densities of P̃kσ (t), P̃z (t) and v(t) in the

periodic state of the model are shown in Fig. 4d. In the model, P̃kσ (t) and P̃zσ (t) have only even and odd harmonics
respectively. Then the velocity v(t) has only odd harmonics according to Eq. (11), as confirmed in Fig. 4d. Since thethird harmonic peak of Svv is about 1/100 of the peak at the fundamental frequency 300 Hz (Fig. 4d), it is reasonable
that we did not see the odd higher harmonics in the experiment (Fig. 3b). In contrast, unlike the experiment, Svv (f )in the model does not show a low-frequency peak, implying that the reciprocating motion in this model does not
involve the long-time persistent motion that was observed experimentally (Sec. 3). We consider that this is because
the model does not capture the asymmetry of the surface charging rate at the upside and downside of the particle
near the electrode correctly. Zhang et al. pointed out this asymmetry and took it into account in their model, which
successfully explained the spontaneous oscillation of the Quincke roller under a strong DC field [24].
Here we introduce a biased model that incorporates the top-bottom asymmetry. For simplicity, we consider theeffect of the asymmetry by adding a constant biased term to the polarization in z-direction: Pz (t) = Pzσ (t)+Pz (t)+Pzα ,where Pzα = 4π0 E0 α with a constant dimensionless value α. This modifies Eqs. (13) and (14) as follows:dP̃zσ+ P̃zσ = −4π (χ∞ + 1/2) cos ω̃ t̃ + χ∞ α − χ∞ ω̃ sin ω̃ t̃
dt̃
+ A(P̃kσ )2 cos ω̃ t̃,

dP̃σk
dt̃

+ P̃σk = −AP̃σk cos ω̃ t̃(P̃zσ + 4πχ∞ α).

(15)
(16)

Note that these equations are not equivalent to the case with a nonzero DC electric field E(t) = E0 cos ωt + E1 evenformally. There is no DC component in the r.h.s. of Eqs. (15) and (16), but the DC component appears due tononlinearity of the equations. Another note is that the steady solution of these equations in the DC limit (ω → 0)
just results in the shift of the stable fix point from that of Eqs. (9) and (10). The modified equations (Eqs. (15) and(16)) indicate that both P̃kσ and P̃zσ can be neither odd nor even functions. It means that P̃kσ has a nonzero odd
component which generates even components in the velocity. The autospectra of the biased model with α = 1 areshown in Fig. 4e. The polarization and velocity have both even and odd peaks. Figure 4f shows the positions of aparticle x(t) for the models with no bias (α = 0) showing no net motion and a bias (α = 0) showing a net motion.
Because there is a nonzero DC component of the velocity in the biased model, the reciprocating motion accompanies
net motion. Although this model does not involve the rotational diffusion, the existence of DC component of velocity
leads to the Lorentzian by considering the athermal rotational coefficient which we discussed in Sec. 3. Therefore, the
ABP-like long-time persistent motion we observed experimentally (Sec. 3) is accounted for by the asymmetry due tothe electrode.
5.

MULTI-PARTICLE BEHAVIOR

Besides single-particle motion, we also observed formation of clusters of multiple particles, which were then maintained over 100–1000 periods or longer in our dilute system (see Video S3, ESI†). To confirm the cluster formation,
we first check spatial distribution of particles in a wider field of view (1031.2 µm square region) using a x10 objective
lens (Olympus, UPlanFL10X, NA0.30) and at 250 Hz. We measure the radial distribution function
g(r) =

X1
δ(|xij | − r),
Np (Np − 1)

(17)

i6=j

where Np is the number of the particles in the field of view, and xij := xi − xj with xi being the position of the ith
particle. Then we evaluate the correlation integral C(R) [43],
C(R) =

Z R

g(r)dr(18)

0

and compare with that of randomly placed particles, generated numerically by the uniform distribution. As shown in
Fig. 5a, in the uniformly random case, C(R) increases with a power law RD with D = 2 in the two-dimensional case[43]. By contrast, the experimental result differs significantly. The number of pairs with short interparticle distanceis significantly larger, indicating the formation of clusters. Based on this result, we consider that a pair of particles
forms a cluster if the interparticle distance is shorter than rth = 35 µm (black dotted line in Fig. 5a, chosen in sucha way that C(R) for R < rth is markedly different from that of the uniformly random case). Using this threshold, we8FIG. 5. Dynamic clustering and bound states of doublets and triplets. (a) Correlation integral C(R) for the experiment (250 Hz,blue) and for the numerically generated, randomly distributed particles (red). We used low-magnification (x10 objective lens)time-lapse images, taken at an interval of 100T to ensure statistical independence between images. The time-averaged numberof particles was N p ' 52. Error bars are the standard errors from 10 slices. In simulations, the same number of particles wererandomly placed in a box of the same size (1031.2 µm square), and average is taken over 10000 realizations. Error bars are thestandard errors from those realizations. The vertical dotted line indicates R = rth = 35 µm. (b) Cluster size distribution for
the experiment (250 Hz, blue) and for the numerically generated random particles (red). Same data as (a) are used both forthe experiment and for the random case. (c) The MSDs of isolated particles, and the MSDs of the centroids of the particles
forming doublets and triplets: MSDc (τ ). The averages over a few samples are shown. The number of samples was three forisolated particles, two for doublets (300 Hz), one for triplet (300 Hz), three for doublets (250 Hz), and one for triplet (250 Hz).
(d) The MSRD of the particle pairs forming doublets and triplets (solid line): MSRD(τ ). For the triplets, averages of MSRD(τ )over all pairs i, j = 1, 2, 3 are shown. Sample average was also taken, as in (c). The dashed gray line is MSRD(τ ) for a particle
pair that initially formed a doublet, but collapsed during the observation time. The dashed black lines show MSRD(τ ) between
a single particle and one of the particles forming a doublet, which join during the observation time to form a triplet (seeVideo S12, ESI†). From the same set of the three particles, MSRD(τ ) between the two particles that constituted the doublet
at the beginning is shown by the dashed red line.

obtain the cluster size distribution (the number of clusters of size n, evaluated by the freud library [44]). The resultis shown in Fig. 5b, with the comparison to the uniformly random case. This shows that clusters (n ≥ 2) are indeed
formed more frequently than the uncorrelated random case.
Now we focus on clusters with n = 2 (doublets; see Videos S4–S7, ESI†) and n = 3 (triplets; see Videos S8–S11,
ESI†). First, note that such clusters are dynamic: clusters are not trapped at fixed locations on the substrate, butinstead each constituent particle moves similarly to isolated particles. This is confirmed by the MSD of the centroidof the particles forming a doublet or a triplet:
MSDc (τ ) =

N −m

X
1
(xc (tk+m ) − xc (tk ))2 ,
N −m

(19)

k=1

where xc := (x1 + x2 )/2 for doublets and xc := (x1 + x2 + x3 )/3 for triplets, and the lag time is τ = (m/3000) s for
300 Hz and τ = (m/2500) s for 250 Hz. This MSD shows oscillations in the short-time region and superdiffusion in a
relatively long time region, similarly to that of isolated particles (Fig. 5c).
Doublets and triplets are bound for some length of times. To confirm the binding, we define the mean squared

9FIG. 6. Experimental results of doublets (150 V). (a) An example of trajectories of a doublet (300 Hz, over 1091 periods= 3.639 s). The insets show the positions recorded every half period for the first 10 periods. (b) Distance |x21 (t)| and relative
angle ϕ(t) (azimuth of x21 (t)) between the two particles of the doublet shown in (a). The inset shows an enlargement displaying
the first 10 periods, with different magnification rates for |x21 (t)| and ϕ(t). (c) Histograms of the inner product of the selfpropelling directions of doublets, n̂1 · n̂2 . Errorbars are the standard errors evaluated from two doublets (300 Hz) and threedoublets (250 Hz).

relative distance (MSRD) of the constituent particles forming a doublet or a triplet:
MSRD(τ ) =

N −m

X
1
(|xij (tk+m )| − |xij (tk )|)2 ,
N −m

(20)

k=1

which quantifies if the relative distance grows or not. The values of MSRD(τ ) for such doublets and triplets that
kept r < rth during the observation time (1091T ) are shown in Fig. 5d (solid line). All show the sign of saturation,suggesting that the particles are bound, more visibly for doublets than triplets. However, the binding state of doublets
and triplets does not last indefinitely; they can also collapse spontaneously or by an approach of other particles. The
dashed gray line in Fig. 5d shows MSRD(τ ) of the doublet that collapsed during the same observation time. The
dashed black lines show MSRD(τ ) between a single particle and one of the particles forming a doublet, which join
during the observation time to form a triplet (see Video S12, ESI†). These datasets of MSRD(τ ) do not saturate
because the corresponding interparticle distance did not always satisfy r < rth during the observation time. In
contrast, from the same triplet that was formed by merging of a doublet and a single particle, MSRD(τ ) for thetwo constituent particles of this doublet saturated (dashed red line in Fig. 5d) as expected, because this pair wasbound during the entire observation time. Note that only two out of eleven clusters (doublets and triplets) showedspontaneous collapse during our observation time (1091T ); this suggests that their life time is in the order of 103 T orso.
In the following, we characterize structure and dynamics of a doublet and a triplet (Secs. 5 A and 5 B.), compare
the motion of a constituent particle with that of an isolated particle (Sec. 5 C), and discuss the interaction between
constituent particles (Sec. 5 D).
A.

DoubletFirst we focus on doublets, which were at least 100 µm away from other particles. We analyzed five doublets(two for 300 Hz and the other three for 250 Hz) which were kept bound during the whole observation time 1091T .
Trajectories of the two particles of a doublet (300 Hz) are shown in Fig. 6a, as an example. The insets show the10

FIG. 7. Experimental results of triplets (150 V). (a) An example of trajectories of the three particles forming a triplet (over1091 periods = 3.639 s, 300 Hz). The insets show the positions recorded every half period for the first 10 periods. (b) Distance
between each pair, |xij (t)|. The inset shows an enlargement displaying the first 10 periods. (300 Hz and 250 Hz) (c) Inner
angles of the triplet, αi (t). (300 Hz and 250 Hz) (d) Histogram of the inner products of the propulsion vectors. Errorbars aretaken for three constituents of a triplet. (e) Histograms of the polar and nematic order parameters of the propulsion vectors,
in the experiment (blue: 300 Hz, orange: 250 Hz) and in the case of three random unit vectors (pink).

positions recorded every half period for the first 10 periods. Similarly to isolated particles, the doublet constituentsalso exhibited short-time reciprocating behavior (at the frequency of the external field) and long-time persistentmotion (see Videos S4–S7, ESI†). Despite the persistent motion of each particle, the interparticle distance |x21 | wasmaintained in some range as shown in Fig. 6b, with noticeable oscillations at a few Hz in addition to the smaller onesat the external frequency (Fig. 6b inset). The relative angle ϕ(t), defined by the azimuth of x21 , also changed in time
(green curve and right axis of Fig. 6b). To characterize the correlation between the motions of the two particles, wedefine the unit vector indicating the self-propelling direction over a half period, called the propulsion vector hereafter,n̂i (tm ) := ∆xi (tm )/|∆xi (tm )| with ∆xi (tm ) := xi (tm+1 ) − xi (tm ). (tm := mT /2, m = 0, 1, 2 . . . ). Figure 6c shows thehistograms of the inner product n̂1 · n̂2 , indicating that the two particles in a doublet tend to align their self-propelling
directions, mostly parallel but sometimes anti-parallel. This suggests that they also tend to align their polarization
vectors Pσk , according to Eq. (11) of the single-particle theory. Further discussions on the interaction will be given inSec. 5 D.

B.

TripletHere we focus on triplets, which were at least 100 µm away from other particles. We analyzed two triplets (one for
300 Hz and the other for 250 Hz) which were kept bound during the whole observation time 1091T . Particle trajectories
are shown in Fig. 7a. The insets show the positions recorded every half period for the first 10 periods. Similarly
to isolated particles and doublets, the particles forming a triplet also showed short-time reciprocating behavior andlong-time persistent motion (see Videos S8–S11, ESI†). The interparticle distance |xij | of any pair i, j = 1, 2, 3 wasmaintained in some range as shown in Fig. 7b. The geometric configuration was also characterized in terms of theinner angles αi of the triangle, whose vertices are set by the positions of the three particles (see the sketch in Fig. 7c).
The result in Fig. 7c shows that the triangle tends to be obtuse most of the time. The propulsion vectors n̂i are also
defined as in the doublet case and analyzed. Figure 7d shows the histogram of the inner products n̂i · n̂j . This indicatesthat any pair of particles tends to align their self-propelling directions, mostly parallel but sometimes anti-parallel.
To further characterize the relation among three propulsion vectors, we use the following two order parameters: the11FIG. 8. Comparison of energy spectral density among isolated single particles (red), particles forming doublets (blue: 300 Hz,
purple: 250 Hz), and particles forming triplets (green: 300 Hz, orange: 250 Hz). For each spectrum, average is taken overall directions (vx and vy ) and all constituent particles of all samples. Three isolated particles, two doublets at 300 Hz, threedoublets at 250 Hz, one triplet at 300 Hz and one triplet at 250 Hz are used. The inset shows an enlargement of a low-frequency
region. To reduce statistical fluctuations, the time series were divided into 30 segments (each segment contains 36 periods),
then the power spectra were averaged.

polar order parameter ΦP and the nematic order parameter ΦN , defined by3ΦP =1 X iθk
e
,
3
k=1

3

ΦN =

1 X 2iθk
e
,
3

(21)

k=1

where θk is the azimuth of the propulsion vector n̂k . We consider four states as shown in Fig. 7e, characterized by ΦP
and ΦN being ≤ 0.5 or > 0.5, and measured the fraction of time spent in each state for 300 Hz and 250 Hz (Fig. 7e).
For comparison, the time fraction for the case of three randomly generated unit vectors is also shown. The result in
Fig. 7e shows that the triplet tends to be in the aligned state (ΦP , ΦN > 0.5).

C.

Comparison between isolated and clustering particles

As we have seen, particles show similar short-time reciprocating motion and long-time persistent motion, regardlessof whether they form a cluster or not. We discuss similarities and dissimilarities between isolated particles andclustering particles here. First, let us see the MSDs for isolated particles as well as the centroids of the doublets andtriplets as already shown in Fig. 5c. While all these results show oscillations in the short-time region and superdiffusion
in a relatively long time region, we also notice a few differences. First, the oscillations in the short-time region are
more prominent for the triplet, then for the doublet, and least prominent for the isolated particles. This suggests thatthe particle interaction contributes to the larger amplitude of the reciprocating motion. Moreover, the time region
showing the superdiffusive exponent is shorter for the clusters than for the isolated particles. This is presumably
because of a binding effect of the clusters.
We also compare the energy spectral density (Fig. 8). Similarly to that of the isolated particles, the spectra for the
doublets and triplets also show a sharp peak at the external frequency as well as a low-frequency Lorentzian mode(inset of Fig. 8). Besides, however, the spectra for the doublet and triplet also show the considerable second and third
harmonic peaks (600 Hz and 900 Hz for 300 Hz, 500 Hz and 750 Hz for 250 Hz, respectively) which were absent orinvisible in the isolated case. Since the third harmonic peak is theoretically expected to appear even in the isolatedcase (see Sec. 4 and Fig. 4d), we consider that it became visible here because of the reduction of the backgroundnoise in the high frequency region (see Fig. 8). This reduction of noise is possibly because the interaction increased
the stability of the dynamics of the polarization and the velocity. In contrast, the apparent second harmonic peak ischaracteristic of the clustering particles; for isolated particles, it did not appear in the experiment, and is expected
to be very weak in theory even if we consider the biased single-particle model in Sec. 4. More precisely, the strengthof the second harmonic peak in this model is no more than 1/1000 of that of the fundamental frequency, so that itis reasonable that we did not see it. Therefore, the considerable second harmonic peak is inferred to be due to theinteraction discussed in the next subsection.

12
D.

Discussions on the interaction between particles

Similarly to the DC case, Quincke rollers under an AC field also interact through hydrodynamic flow and electrostatic
effect. As described in Bricard et al. [21], the hydrodynamic flow generated by a roller at xj is a dipolar flow for
xij := |xij |  H, where H is the cell thickness (H = 30 µm here). In the vicinity of the particle (xij . H), aligningflow is produced [21]. The reciprocating motion under the AC field simply makes the flow field oscillate. On theelectrostatic interaction, the parallel component of the electric field at xi induced by the polarization P(j) of particle(j)j, denoted by δEk (xi , t), is given by [21]3(j)
δEk (xi , t) =2π0 x3ij

"

#
a (j)
a3
a2 (j)
P x̂ij − 2 Pk · (5x̂ij x̂ij − I) + O( 3 ) ,
xij z
xij
xij

(22)

with x̂ij := xij /|xij | and the identity matrix I. Thus, particle i is exposed to the total electric field E0 cos ωt +
P
(j)
(j)
(j)
and Pk have odd and even harmonics,
j δEk (xi , t), which should replace E0 in Eqs. (4) and (5). Since Pz(j)respectively (see Sec. 4 and Fig. 4d), δEk (xi , t) has both of them, and so does Pσ of particle i. This accounts forthe considerable second harmonic peak observed in the energy spectral density only for clustering particles (Fig. 8).
We consider that the stable formation of clusters is also a result of such hydrodynamic and electrostatic interactions
under the AC external field. Interestingly, clustering has also been reported for DC Quincke rollers with periodic
on-off switching [23] but not under a constant DC field [21, 22]. Therefore, time-dependent propulsion and interaction
may be a key to clustering. Elucidating the clustering mechanism is an important open problem left for future studies.

6.

CONCLUSION AND OUTLOOKOur experimental investigation into Quincke rollers under an AC electric field revealed the reciprocating motion at
the frequency of the external field, accompanied by ABP-like long-time persistent motion. The single-particle theory
explained how this reciprocating motion arises as a result of oscillations of the polarization, which turned out tohave different parities between the vertical and parallel components. The ABP-like long-time persistent motion isconsidered to arise by the effect of charging asymmetry at the upside and downside of the particle near the electrode.
The extreme deviation between the thermal and observed values of the rotational diffusion coefficient can be attributed
to athermal noise caused by the interplay between surface roughness of the electrodes and electrokinetic effect at thesurface of particles. We also found formation of clusters such as doublets and triplets, and characterized their structureand dynamics. We showed that clustering results in the generation of the considerable second harmonic peak in theenergy spectral density, which we accounted for by the dipole-dipole interaction between particles.
There are a number of interesting directions to take for future studies of AC Quincke rollers. Experimentally,
it is tempting to study formation of larger clusters in denser suspension, as well as dynamical processes of cluster
formation. Studying the effect of AOT concentration, as reported by Zhang et al. [24] for strong DC fields, may also
be interesting in the AC case. The effect of AOT micelles may be incorporated in our model, in terms of positiondependent conductivity as modeled by Zhang et al. [24]. Simulations of a collection of particles interacting throughthe hydrodynamic and electrostatic interactions formulated by Bricard et al. [21], but under an AC field, are also
desired.
From broader perspectives, studies of artificial bipolar SPPs have not been as developed as those for the unidirectional, polar counterparts. The system of AC Quincke rollers presented in this work may play a unique role in thiscontext, thanks to their regular reciprocation that takes place at the externally controlled frequency. It is known thatfor some types of self-propelled Janus particles the self-propulsion velocity may be reversed by changing parameterson the external energy source (such as the frequency of the applied voltage, the strength of the light irradiation, etc.)[19, 45–47], but the propulsion mechanism and the particle interaction are typically very different in the reversed
state, so that symmetric reciprocal motion cannot be realized by simply switching between the normal and reversed
states back and forth. By contrast, the AC Quincke rollers studied in the present work have advantages that thereciprocal motion is almost symmetric, owing simply to the electric field reversal. While having fundamental aspects
in common with other active matter systems, such as self-propulsion of particles and energy dissipation at multiple
scales, periodically driven active systems may serve as a new class of active matter worth exploring. We hope that
it will develop and contribute to deeper understanding of active matter and related phenomena in biology and nonequilibrium physics, and that controllable bipolar SPPs such as AC Quincke rollers will serve as a useful tool in thisdirection.

13
AUTHOR CONTRIBUTIONSANK and MS designed the experiment. KAT contributed to continue the research. ANK performed experiments andnumerical simulations. ANK and MS contributed to theoretical calculations. All authors discussed and contributed
to the analysis and the interpretation of the results. All authors wrote the manuscript.

CONFLICTS OF INTEREST

There are no conflicts to declare.

ACKNOWLEDGEMENTS

We thank D. Nishiguchi for his help to prepare the cell, and H. Wada, N. Yoshinaga, F. Peruani, H. Kitahata for
fruitful discussions. We thank Ge-Wei Chen and the analysis center of SJTU for helping to measure surface roughnessof ITO film. This work is supported by Grand-in-Aid for JSPS Fellows (Grant No. JP17J06659) and KAKENHI(Grant Nos. JP25103004, JP19H05800, JP20H01826).

[1] 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).
[2] C. Bechinger, R. Di Leonardo, H. Löwen, C. Reichhardt, G. Volpe, and G. Volpe, Rev. Mod. Phys. 88, 045006 (2016).
[3] A. Doostmohammadi, J. Ignés-Mullol, J. M. Yeomans, and F. Sagués, Nature Communications 9, 3246 (2018).
[4] H. Chaté, Annu. Rev. Condens. Matter Phys. 11, 189 (2019).
[5] T. Vicsek, A. Czirók, E. Ben-Jacob, I. Cohen, and O. Shochet, Phys. Rev. Lett. 75, 1226 (1995).
[6] I. Theurkauff, C. Cottin-Bizonne, J. Palacci, C. Ybert, and L. Bocquet, Phys. Rev. Lett. 108, 268303 (2012).
[7] J. Palacci, S. Sacanna, A. P. Steinberg, D. J. Pine, and P. M. Chaikin, Science 339, 936 (2013).
[8] I. Buttinoni, J. Bialké, F. Kümmel, H. Löwen, C. Bechinger, and T. Speck, Phys. Rev. Lett. 110, 238301 (2013).
[9] F. Ginot, I. Theurkauff, F. Detcheverry, C. Ybert, and C. Cottin-Bizonne, Nature Communications 9, 696 (2018).
[10] J. Toner and Y. Tu, Phys. Rev. Lett. 75, 4326 (1995).
[11] S. Ramaswamy, R. A. Simha, and J. Toner, Europhys. Lett. 62, 196 (2003).
[12] V. Narayan, S. Ramaswamy, and N. Menon, Science 317, 105 (2007).
[13] J. Deseigne, O. Dauchot, and H. Chaté, Phys. Rev. Lett. 105, 098001 (2010).
[14] 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).
[15] Y. Wang, R. M. Hernandez, D. J. Bartlett, J. M. Bingham, T. R. Kline, A. Sen, and T. E. Mallouk, Langmuir 22, 10451(2006).
[16] S. Gangwal, O. J. Cayre, M. Z. Bazant, and O. D. Velev, Phys. Rev. Lett. 100, 058302 (2008).
[17] H.-R. Jiang, N. Yoshinaga, and M. Sano, Phys. Rev. Lett. 105, 268302 (2010).
[18] D. Nishiguchi and M. Sano, Phys. Rev. E 92, 052309 (2015).
[19] T. Mano, J.-B. Delfau, J. Iwasawa, and M. Sano, Proceedings of the National Academy of Sciences 114, E2580 (2017).
[20] J. Iwasawa, D. Nishiguchi, and M. Sano, Phys. Rev. Research 3, 043104 (2021).
[21] A. Bricard, J.-B. Caussin, N. Desreumaux, O. Dauchot, and D. Bartolo, Nature 503, 95 (2013).
[22] A. Bricard, J.-B. Caussin, D. Das, C. Savoie, V. Chikkadi, K. Shitara, O. Chepizhko, F. Peruani, D. Saintillan, and
D. Bartolo, Nature Communications 6, 7470 (2015).
[23] H. Karani, G. E. Pradillo, and P. M. Vlahovska, Phys. Rev. Lett. 123, 208002 (2019).
[24] Z. Zhang, H. Yuan, Y. Dou, M. O. de la Cruz, and K. J. M. Bishop, Phys. Rev. Lett. 126, 258001 (2021).
[25] K. Nagai, Y. Sumino, H. Kitahata, and K. Yoshikawa, Phys. Rev. E 71, 065301 (2005).
[26] H. Ebata and M. Sano, Scientific Reports 5, 8546 (2015).
[27] C. Lozano, J. R. Gomez-Solano, and C. Bechinger, New Journal of Physics 20, 015008 (2018).
[28] B. Qian, D. Montiel, A. Bregulla, F. Cichos, and H. Yang, Chemical Science 4, 1420 (2013).
[29] H. C. Berg and D. A. Brown, Nature 239, 500 (1972).
[30] Y. Wu, A. D. Kaiser, Y. Jiang, and M. S. Alber, Proceedings of the National Academy of Sciences 106, 1222 (2009).
[31] K. Kawaguchi, R. Kageyama, and M. Sano, Nature 545, 327 (2017).
[32] S. Thutupalli and S. Herminghaus, The European Physical Journal E 36, 91 (2013).
[33] G. Quincke, Annu. Phys. Chem. 59, 417 (1896).
[34] Applying AC electric fields has the additional advantage of not damaging electrodes, in contrast applying DC electric fielddamages electrodes in about 10 minutes due to chemical reactions on the surface.

14[35] I. Sbalzarini and P. Koumoutsakos, Journal of Structural Biology 151, 182 (2005).
[36] A. Squarcini, A. Solon, and G. Oshanin, New Journal of Physics 24, 013018 (2022).
[37] A. Goldman, R. Cox, and H. Brenner, Chemical Engineering Science 22, 637 (1967).
[38] A. Goldman, R. Cox, and H. Brenner, Chemical Engineering Science 22, 653 (1967).
[39] M. E. O’Neill and K. Stewartson, Journal of Fluid Mechanics 27, 705–724 (1967).
[40] Q. Liu and A. Prosperetti, Journal of Fluid Mechanics 657, 1–21 (2010).
[41] J. R. Melcher and G. I. Taylor, Annual Review of Fluid Mechanics 1, 111 (1969).
[42] N. Pannacci, L. Lobry, and E. Lemaire, Phys. Rev. Lett. 99, 2 (2007).
[43] P. Grassberger and I. Procaccia, Phys. Rev. Lett. 50, 346 (1983).
[44] V. Ramasubramani, B. D. Dice, E. S. Harper, M. P. Spellings, J. A. Anderson, and S. C. Glotzer, Computer Physics
Communications 254, 107275 (2020).
[45] H. R. Vutukuri, M. Lisicki, E. Lauga, and J. Vermant, Nature Communications 11, 2628 (2020).
[46] J. R. Gomez-Solano, S. Samin, C. Lozano, P. Ruedas-Batuecas, R. van Roij, and C. Bechinger, Scientific Reports 7, 14891(2017).
[47] R. Suzuki, H. R. Jiang, and M. Sano, arXiv:1104.5607 (2011).

Electronic Supplementary Information for“Active colloid with externally induced periodic
bipolar motility and its cooperative motion”
arXiv:2203.12882v3 [cond-mat.soft] 14 Jul 2022

Airi N. Kato, Kazumasa A. Takeuchi, and Masaki Sano

Supplementary Note 1: Model of a single self-propelled particle
with a sinusoidal and constant velocity
Consider a particle self-propelling at speed v(t) under the effect of translational and rotational diffusion.
The translational diffusion and rotational diffusion coefficients are denoted by D and Dθ , respectively.
Let x and θ be the position and the propelling direction, respectively. The equations of motion read:ẋ(t) = v(t)n̂(t) + ζ (t)(S1)

θ̇(t) = ξ(t).

(S2)

Here, n̂(t) = (cos θ(t), sin θ(t)) is the propelling direction vector. The noises are white Gaussian:hζi (t)ζj (t0 )i = 4Dδ(t − t0 )δij , hζi (t)i = 0, hξ(t)ξ(t0 )i = 2Dθ δ(t − t0 ), and hξ(t)i = 0 for arbitrary t, t0 .
Here let us assume the self-propulsion velocity v(t) to bev(t) = v0 + v1 cos ωt.

(S3)

We calculate the energy power spectrum of this model via the velocity autocorrelation. The velocity
autocorrelation is calculated as follows.
hv(t1 ) · v(t2 )i =v(t1 )v(t2 )hn̂(t1 ) · n̂(t2 )i + 4Dδ(t1 − t2 )= v(t1 )v(t2 )hcos θ(t1 ) cos θ(t2 ) + sin θ(t1 ) sin θ(t2 )i + 4Dδ(t1 − t2 )=v(t1 )v(t2 )hcos(θ(t1 ) − θ(t2 ))i + 4Dδ(t1 − t2 )= v(t1 )v(t2 )hcos(∆t1 −t2 θ)i + 4Dδ(t1 − t2 ),(S4)

where ∆τ θ := θ(t + τ ) − θ(t) and the time translation invariance of θ(t) is taken into account. The
probability density of ∆τ θ, denoted by P (∆θ, τ ), has the characteristic function that satisfies φ(s, t) :=R∞
2dθP (θ, t)eiθs = e−Dθ s t . For τ > 0,−∞hcos(∆τ θ)i =ZP (θ, τ ) cos θdθφ(1, τ ) + φ(−1, τ )2= e−Dθ τ .
=(S5)Therefore, neglecting the translation diffusion, we obtain the following expression for the velocity
autocorrelation:
hv(t1 ) · v(t2 )i ==v(t1 )v(t2 )e−Dθ |t1 −t2 | 2
v0 + v0 v1 {cos(ωt1 ) + cos(ωt2 )} + v12 cos(ωt1 ) cos(ωt2 ) e−Dθ |t1 −t2 |(S6)which cannot be expressed as a function of (t1 − t2 ). In other words, hv(t1 ) · v(t2 )i is not invariant
under time translation.
1Generally, the energy spectral density is calculated by the Fourier transformation of the velocity
autocorrelation hv(t) · v(t + τ )i. In our model, however, hv(t) · v(t + τ )i cannot be expressed by afunction of τ only. Therefore, here we consider the periodically-averaged velocity correlation C(τ ) :=hv(t) · v(t + τ )i, where the overline denotes averaging over a period of the external field. This isevaluated as follows:C(τ ) =
v02 + v12 cos(ωt) cos(ω(t + τ )) e−Dθ τv2=v02 + 1 (cos(2ωt + ωτ ) + cos(ωτ )) e−Dθ τ
2


v12
2
=
v0 +
cos(ωτ ) e−Dθ τ .
(S7)
2Then, the spectral density S(f ) is obtained by Fourier transform of Eq. (S7), as follows:


v12
Dθ
2Dθ
Dθ
2
+
+ 2
,
S(f ) = v0 2
Dθ + (2πf )2
2 Dθ2 + (ω + 2πf )2Dθ + (ω − 2πf )2(S8)

when the translational diffusion is neglected. The translational diffusion adds a constant terms to
Eq. (S8). For low frequency f , the first term is dominant in Eq. (S8), which is Lorentzian with theθcutoff frequency f ∼ D
2π .
To estimate Dθ , we consider the experimentally observed low-frequency region to be a simple
2v 2
2v 2
2Dθ
Lorentzian: the first term of Eq. (S8) v02 D2 +(2πf= Dθ0 · 1+(2πf1 /Dθ )2 . The coefficient Dθ0 was)2θdetermined by the value at f = 0. By fitting the data for f < 40 Hz, we obtained Dθ ∼ 31.4 ± 0.6 s−1(Dθ−1 ∼ 0.03 s). Note that this effective rotational diffusive coefficient Dθ is athermal, can beoriginated from such as surface heterogeneity of the electrode and particles. (Cf. thermal rotational
kB T
−3 −1
diffusive coefficient: Dθthermal = 8πas with the room temperature T and the viscosity of3 η ∼ 10
hexadecane η.)

Supplementary figures

Figure S1: Periodic states of the single-particle model for different frequencies of the external electricfield. Trajectories from initial conditions with P̃kσ (0) > 0 are shown, in (P̃kσ , P̃zσ ) space. For highenough frequencies, P̃kσ vanishes and the particle does not reciprocate (see Eq.(10) in the main paper).

2

Figure S2: Experimental results of the energy spectral density of particles forming a doublet, fordifferent frequencies of the external electric field. The amplitude was kept at 150 V (root-meansquare amplitude). The motion of the particles was captured at the frame rate 10 times faster thanthe external frequency (e.g., 2500 fps for 250 Hz). The highest peak of each spectrum is found at
the frequency of the external field, corresponding to the reciprocating motion of the particles. The
second harmonic peak due to the dipole-dipole interaction is also seen in all cases. Panel (b) shows anenlargement of (a) in a low-frequency region. To reduce statistical fluctuations, the time series weredivided into 30 segments (each segment contains 36 periods), then the power spectra were averaged.
All spectra show a similar low-frequency mode.

Supplementary Note 2: Possible effects of surface roughness of
ITO electrode and its interaction with particle surface
To discuss the possibility that the athermal value of the rotational diffusion coefficient results fromsurface heterogeneity of the electrode, we measured the surface roughness of the indium-tin oxide (ITO)
film on the glass substrate used in the experiments (Mitsuru Optical Co. Ltd.) with an atomic force
microscope (AFM, Bruker, Dimension FastScan Bio). The standard deviation of the height fluctuations
is 3.8 nm (Fig. S3a) . Histogram of the height fluctuations is close to a Gaussian distribution andthe power spectrum of the height fluctuations is close to a Lorentzian. From the fitting, we obtained
the characteristic wavelength to be λ ∼ 0.75 µm. External mechanical perturbations to the moving
particle on a rough substrate can be estimated from potential energy calculations. The potential energy
fluctuations calculated by assuming that the particle is moving along the ITO surface landscape areless or of the same order of the thermal fluctuation kB T . However, surface roughness of the order of
4 nm is close to the size of Debye length of the surface double layer of the particle. This Debye length
level of roughness may cause even stronger fluctuations in the flow of charge around the particle and its
distribution, which may have a significant impact on the effective rotational diffusion of the particle.
Since the mean speed of the particles averaged over the period of the externally applied electric field isabout vm ∼ 60 µm/s and the characteristic time of the Lorentzian in Fig.3b is τ ∼ 0.032 s, the numberof fluctuation events experienced by a particle within τ is estimated by n ∼ vm τ /λ. We obtainedn ∼ 2.4. This number of fluctuation events is sufficient to cause loss of orientation information andresult in athermal rotational diffusion of the particle.

3

Figure S3: Measurement of surface roughness of ITO substrate. (a) Height roughness of the ITO coaton the glass substrate measured with an AFM (xy-area size: 10 × 10 µm2 ). Z-range of color-coded
image is 32.7 nm. (b) Histogram of height fluctuations. Standard deviation of the height fluctuation is
3.8 nm. (c) Radial distribution of the height fluctuation power spectrum in k-space calculated from thetwo-dimensional Fourier transformation of the height fluctuations. Lorentzian fit gives a characteristic
wavenumber of 4.8 µm−1 . (wavelength ∼ 0.75 µm)

Supplementary videos• Video S1 (S1.mov): Short-time behavior of an isolated particle. The particle reciprocates at the external frequency (300 Hz), though it is only barely visible in the movie. The movieis played at 0.01 times the real speed.
• Video S2 (S2.mov): Long-time behavior of an isolated particle (stroboscopic). This
movie shows long-time persistent motion of the particle. The frame acquisition rate is set to beequal to the external frequency (300 Hz). The movie is played at 0.1 times the real speed.
• Video S3 (S3.mov): A wide-field video (x10 objective) showing multiple clusters(stroboscopic, 250 Hz). Clusters are formed and maintained longer than 1000 periods if they
do not encounter other particles. The movie is played at 0.1 times the real speed.
• Video S4 (S4.mov): Short-time behavior of a doublet (300 Hz). The two particles
reciprocate at the external frequency in a cooperative manner. The propelling directions tend toalign, but not always. The movie is played at 0.01 times the real speed.
• Video S5 (S5.mov): Long-time behavior of a doublet (stroboscopic, 300 Hz). The twoparticles show persistent motions, with the interparticle distance fluctuating but kept in somerange. The movie is played at 0.1 times the real speed.
• Video S6 (S6.mov): Short-time behavior of a doublet (250 Hz). The two particles
reciprocate at the external frequency in a cooperative manner. The propelling directions tend toalign, but not always. The movie is played at 0.01 times the real speed.
• Video S7 (S7.mov): Long-time behavior of a doublet (stroboscopic, 250 Hz). The twoparticles show persistent motions, with the interparticle distance fluctuating but kept in somerange. The movie is played at 0.1 times the real speed.

4• Video S8 (S8.mov): Short-time behavior of a triplet (300 Hz). The three particles
reciprocate at the external frequency (300 Hz) in a cooperative manner. The propelling directions
tend to align, but not always. The movie is played at 0.01 times the real speed.
• Video S9 (S9.mov): Long-time behavior of a triplet (stroboscopic, 300 Hz). The threeparticles show persistent motions, with the interparticle distances fluctuating but kept in somerange. The movie is played at 0.1 times the real speed.
• Video S10 (S10.mov): Short-time behavior of a triplet (250 Hz). The three particles
reciprocate at the external frequency in a cooperative manner. The propelling directions tend toalign, but not always. The movie is played at 0.01 times the real speed.
• Video S11 (S11.mov): Long-time behavior of a triplet (stroboscopic, 250 Hz). The
three particles show persistent motions, with the interparticle distances fluctuating but kept insome range. The movie is played at 0.1 times the real speed.
• Video S12 (S12.mov): Long-time behavior of a triplet formation (stroboscopic,
300 Hz). An initially isolated particle joined a doublet and formed a triplet. The shape ofthe triangle formed by the three particles dynamically changes. The movie is played at 0.1 times
the real speed.

5