Hydrodynamic synchronization and collective dynamics of colloidal particles driven along a circular path

TL;DR

This paper theoretically investigates the collective behavior of colloidal particles driven along a circular path by an optical vortex, revealing synchronization phenomena and dynamic rearrangements that align with experimental observations.

Contribution

It derives phase equations considering hydrodynamic and repulsive interactions, providing new insights into the collective dynamics of colloidal particles on a circular trajectory.

Findings

Particles synchronize and form doublets that move faster than singles.

Periodic rearrangements of particle groups occur for certain numbers of particles.

The system exhibits oscillating or stable synchronized states depending on initial conditions.

Abstract

We study theoretically the collective dynamics of particles driven by an optical vortex along a circular path. Phase equations of N particles are derived by taking into account both hydrodynamic and repulsive interactions between them. For N = 2, the particles attract with each other and synchronize, forming a doublet that moves faster than a singlet. For N = 3 and 5, we find periodic rearrangement of doublets and a singlet. For N = 4 and 6, the system exhibits either a periodic oscillating state or a stable synchronized state depending on the initial conditions. These results reproduce main features of previous experimental findings. We quantitatively discuss the mechanisms governing the non-trivial collective dynamics.

Tables

Table 1. Table 1: The maximum, minimum and time-averaged phase velocities v max , v min , v ave subscript 𝑣 max subscript 𝑣 min subscript 𝑣 ave v_{\rm max},v_{\rm min},v_{\rm ave} and period of oscillation T 𝑇 T . For N = 4 𝑁 4 N=4 and 6 6 6 , we show the values for both stable and oscillating (osc.) modes.

	$N=2$	$N=3$	$N=4$(osc.)	$N=4$(stable)	$N=5$	$N=6$(osc.)	$N=6$(stable)
$v_{\max }$	1.58	1.67	1.68	1.46	1.65	1.64	1.54
$v_{\min }$	1.58	0.93	0.89	1.46	1.03	1.07	1.54
$v_{\mathrm{ave}}$	1.58	1.42	1.36	1.46	1.46	1.45	1.54
$T$	-	31.6	18.1	-	31.4	16.4	-

Equations

$${\text{r}}_i=R_i{\text{e}}_i^r=(R_i\cos {\varphi }_i,R_i\sin {\varphi }_i,0),$$

$${\text{e}}_i^r$$

$${\text{e}}_i^{\varphi }$$

$${\Delta }_{ij}={\varphi }_j-{\varphi }_i.$$

$${\text{F}}_i=F_i^r{\text{e}}_i^r+F_i^{\varphi }{\text{e}}_i^{\varphi }.$$

$$F_i^{\varphi }=F=const.,$$

$$F_i^r=-k_r\delta R_i,\delta R_i=R_i-R.$$

$$\kappa =\frac{k_{r}R}{F}$$

$$\eta {\nabla }^2\text{v}-\nabla p$$

$$\nabla \cdot \text{v}$$

$${\text{g}}_i=\gamma [\text{v}({\text{r}}_i)-{\dot{\text{r}}}_i],$$

$$\text{v}({\text{r}}_i)\simeq -\sum _{j\ne i}\text{G}({\text{r}}_i,{\text{r}}_j)\cdot {\text{g}}_j.$$

$$\text{G}(\text{r})=\frac{1}{8\pi \eta r}\left(\text{1}+\frac{\text{r}\text{r}}{r^2}\right)+\frac{a^2}{12\pi \eta r^3}\left(\text{1}-\frac{3\text{r}\text{r}}{r^2}\right).$$

$${\hat{\mathbf{G}}}_{ij}=\gamma \mathbf{G}({\mathbf{r}}_i-{\mathbf{r}}_j)$$

$$\alpha =\frac{a}{R}$$

$$\mathbf{v}({\mathbf{r}}_i)\simeq \sum _{j\ne i}{\hat{\mathbf{G}}}_{ij}\cdot {\dot{\text{r}}}_j.$$

$${\text{F}}_{ji}^{rep}$$

$$F_{ji}^{rep}$$

$$U(\Delta )=U_0\text{quantity}(\frac{R|\Delta |}{2a}-1{)}^{-12},$$

$$u=\frac{U_0}{Fa}.$$

$$F_i^r{\text{e}}_i^r+F_i^{\varphi }{\text{e}}_i^{\varphi }+\sum _{j\sim i}{F}_{ji}^{rep}{\mathbf{e}}_i^{\varphi }=\gamma \sum _{j\ne i}\left({\dot{\text{r}}}_i-{\hat{\mathbf{G}}}_{ij}\cdot {\dot{\text{r}}}_j\right),$$

$$-k_r\delta R_i=\gamma {\dot{R}}_i-\gamma \sum _{j\ne i}{\text{e}}_i^r\cdot {\hat{\text{G}}}_{ij}\cdot {\dot{\text{r}}}_j$$

$$\gamma R_i{\dot{\varphi }}_i=F+{\text{e}}_i^{\varphi }\cdot \sum _{j\sim i}{\text{F}}_{ji}^{rep}+\gamma \sum _{j\ne i}{\text{e}}_i^{\varphi }\cdot {\hat{\text{G}}}_{ij}\cdot {\dot{\text{r}}}_j,$$

$$\omega =\frac{F}{\gamma R}.$$

$$\frac{\delta R_i}{R}\approx {\kappa }^{-1}\sum _{j\ne i}{\text{e}}_i^r\cdot {\hat{\text{G}}}_{ij}\cdot {\text{e}}_j^{\varphi }.$$

$$\frac{{\dot{\varphi }}_i}{\omega }=1$$

$$-\text{quantity}[f^{rep}({\Delta }_{i,(i+1)\mathrm{mod}N})-f^{rep}({\Delta }_{(i-1)\mathrm{mod}N,i})].$$

$$f^{rep}(\Delta )=6u{\left(\frac{|\Delta |}{2\alpha }-1\right)}^{-13}$$

$${\text{e}}_i^r\cdot {\hat{\mathbf{G}}}_{ij}\cdot {\text{e}}_j^{\varphi }=-\frac{9\alpha \sin {\Delta }_{ij}}{16\left|\sin \left(\frac{{\Delta }_{ij}}{2}\right)\right|}+\frac{{\alpha }^3\sin {\Delta }_{ij}}{32{\left|\sin \left(\frac{{\Delta }_{ij}}{2}\right)\right|}^3},$$

$${\text{e}}_i^{\varphi }\cdot {\hat{\mathbf{G}}}_{ij}\cdot {\text{e}}_j^{\varphi }=\frac{\alpha \text{quantity}(3+9\cos {\Delta }_{ij})}{16\left|\sin \left(\frac{{\Delta }_{ij}}{2}\right)\right|}-\frac{{\alpha }^3(3+\cos {\Delta }_{ij})}{32{\left|\sin \left(\frac{{\Delta }_{ij}}{2}\right)\right|}^3}.$$

$$\frac{{\dot{\varphi }}_i}{\omega }=1+\sum _{k\ne i}\Omega ({\Delta }_{ik}).$$

$$\frac{{\dot{\Delta }}_{ij}}{\omega }=\sum _{k\ne j}\Omega ({\Delta }_{jk})-\sum _{k\ne i}\Omega ({\Delta }_{ik}),$$

Full text

I Introduction

Hydrodynamic synchronization is considered to be a key mechanism that controls the collective dynamics of flagella and cilia. For example, bacteria such as E. coli shows run-and-tumble motion induced by dynamical bundling/unbundling of their flagella, and ciliates such as Paramecium regulates the direction of motion using metachronal waves[1]. Metachronal waves of cilia also play vital roles in moving fluid in mammalian bodies (airways, embryo, ventricles, etc.) These phenomena suggest the importance of indirect hydrodynamic interaction between flagella and cilia in their coordinate motion. Complexity due to dynamical regulation by molecular motors is removed by considering an artificial model system. Experiments have been pushed forward using colloidal particles driven by optical tweezers, and verified theoretical scenarios in the emergence of synchronization [2, 3, 4]: dynamical switching of the driving potential (rower model)[5, 6], deformability of the trajectory [7], and periodic modulation of the driving force [8, 9, 10, 11]. The last mechanism enables flexible control of the phase difference between two rotors [2]. Fluid inertia may also induce synchronization in a finite-Reynolds number system [12, 13].

In this paper, we theoretically address the dynamics of colloidal particles driven along a circular path. Non-trivial collective motion have been found in experiments using an optical vortex [14, 15, 16, 17, 18, 19] or optical tweezers [20, 21, 22, 23, 24, 25]. An optical vortex can drive multiple particles simultaneously along a single circular trajectory [26], while optical tweezers manipulate the particles individually. For three particles, a limit-cycle behavior predicted in Ref. [27] is confirmed and effects of non-constant force profiles have been discussed [20, 14]. Sokolov et al. [15] showed that a pair of particles attract each other due to radial displacements from the guided circular trajectory. Increasing the particle number $N$ up to $N=11$, they also found that the collective mobility changes non-monotonically: for odd $N$, one unpaired particle remains and reduces the collective mobility. Sassa et al. [16] performed experiments for up to $N=9$ particles and observed multiple collective states for even $N$: the particles either form stable doublets (stable states) or exhibit periodic rearrangement of doublets and (a) singlet(s) (oscillating states), and switching between these states are observed. More complex dynamics are observed for particles of different sizes[17] and particles in confined geometry [28]. Theoretically, linear stability analysis [27, 16] and direct numerical simulations [15, 17] have been done to capture the essential mechanism of the collective motion. In the present paper, we pursue an analytical approach to gain more quantitative description of the phenomena. Employing the phase reduction method, which was successfully applied to hydrodynamic synchronization in a system with deformable trajectories [7, 2], we obtain the dynamical equations in a compact form. Numerical solutions of the phase equations for up to $N=6$ reproduce the main features of the experimental findings. Construction of the paper is as follows: In Section II, we describe the model and derive the phase equations. In Section III, numerical solution of the phase equations are presented for $N=2-6$ particles. In Section IV, we discuss the results in comparison with the experiments.

II Model

II.1 Geometry and Forces

We consider $N$ colloidal particles of radius $a$ that are driven along a circular trajectory of radius $R$ on the $z=0$ plane. The position of the $i$-th particle ($i=1,2,...,N$) is given by

[TABLE]

where $R_{i}=R_{i}(t)$ is the instantaneous orbital radius of the $i$-th particle, which can deviate from the target radius $R$ set by the optical vortex, and $\phi_{i}=\phi_{i}(t)$ is the phase. Here and hereafter we use the radial and tangential unit vectors

[TABLE]

and

[TABLE]

We label the particles along the counter-clockwise direction so that $\phi_{1}<\phi_{2}<\cdots<\phi_{N}$, and define the phase difference

[TABLE]

The motion of the particles are governed by three types of forces acting on them: the forces due to the optical vortex, the hydrodynamic forces, and the repulsive interaction among the particles.

The force by the optical vortex is decomposed into the radial and tangential components as

[TABLE]

We assume that the tangential driving force is constant,

[TABLE]

while the radial restoring force is linearly proportional to the radial displacement:

[TABLE]

The spring constant $k_{r}$ is determined by the curvature of the optical trapping potential along the radial direction. We will use the dimensionless spring constant

[TABLE]

to describe the relative strength of trapping. We assume strong trapping ($\kappa\gg 1$) and regard $\kappa^{-1}$ as a small perturbative parameter.

The particles move in an incompressible fluid of shear viscosity $\eta$. Since we consider a micrometer-sized system, the Reynolds number is negligibly small (typically $Re\sim 10^{-5}$) and the flow velocity field ${\bf v}({\bf r})$ obeys the Stokes equation,

[TABLE]

The hydrodynamic drag force exerted on the $i$-th particle is given by

[TABLE]

where $\gamma=6\pi\eta a$ is the drag coefficient. The flow velocity is determined by

[TABLE]

Here, $\textbf{G}(\textbf{r}_{i}-\textbf{r}_{j})$ is the off-diagonal elements of the grand mobility tensor and describes the hydrodynamic interaction between the particles. We use the Rotne-Prager-Yamakawa approximation [29, 30],

[TABLE]

The first term on the RHS of (13) is the Oseen tensor, which is the Green function of the Stokes equation, and the second term is the correction due to finite volume of the spheres. Introducing the dimensionless tensor

[TABLE]

which scales as $a/|{\bf r}_{i}-{\bf r}_{j}|$, we find that

[TABLE]

is the dimensionless parameter describing the relative strength of the hydrodynamic interaction. Substituting (12) into (11) and then back into (12), we get an expression of the velocity field to the second order with respect to $\gamma{\bf G}$. This can be continued recursively to give a series expansion of the velocity field in terms of $\alpha$. Neglecting $O(\alpha^{2})$ terms, we get the linear relation

[TABLE]

The intermolecular forces cause excluded volume interaction between the particles. Assuming that the radial displacements of the particles are small compared to their radii, we can neglect the radial component of the repulsive force. Thus the repulsive force exerted by the $j$-th particle on the $i$-th is written in the form

[TABLE]

where $U(\Delta)$ is the potential as the function of the phase difference.

We choose the function [27]

[TABLE]

where $U_{0}$ is a positive constant, and $\Delta$ is the phase difference folded into $[-\pi,\pi]$. Because this potential diverges when the particles are in contact $(\Delta=2a/R)$, we can avoid overlapping of particles in numerical analysis. Since the force decays rapidly as a function of the distance, we retain only the interaction between the nearest neighbor pairs ($i$ and $i\pm 1$ mod $N$). The relative strength of the repulsive force compared to the driving force is given by the dimensionless parameter

[TABLE]

II.2 Phase Equation

Now we derive the time evolution equation for the phase $\phi_{i}(t)$. The equation of force balance is

[TABLE]

where $\sum_{j\sim i}$ means that the sum is taken over particles next to the $i$-th particle. On the RHS, we have the hydrodynamic forces from (11) and (16), in which we substitute the particle velocity $\dot{\bf r}_{i}=\dot{R}_{i}{\bf e}_{i}^{r}+R_{i}\dot{\phi}_{i}{\bf e}_{i}^{\phi}.$ Taking the inner products of $\textbf{e}_{i}^{r}$, ${\bf e}_{i}^{\phi}$ and (21), we get

[TABLE]

and

[TABLE]

respectively, where we also used (6) and (7). In the absence of hydrodynamic and repulsive interactions ($\alpha=u=0$) and in the limit of strong trapping ($\kappa^{-1}=0$), Eq.(23) gives the intrinsic phase velocity

[TABLE]

We also see from (22) that the radial displacement is an $O(\kappa^{-1}\alpha)$ quantity. Then we can use $\dot{\bf r}_{i}\simeq R\dot{\phi}_{i}{\bf e}_{i}^{\phi}\simeq R\omega{\bf e}_{i}^{\phi}$ in the last term on the RHS of (22), to obtain

[TABLE]

For Eq. (23), we use the same approximation on the RHS and divide by $\gamma\omega R_{i}=\gamma\omega R(1+\delta R_{i}/R)$, in which (25) is substituted. Neglecting terms of $O(\kappa^{-2}\alpha^{2},\kappa^{-1}u,\alpha u)$, we arrive at the phase equation

[TABLE]

Here, $\hat{\bf G}_{ij}$ is evaluated by setting $\delta R_{i}=\delta R_{j}=0$ and becomes a function of the phases $\phi_{i}$ and $\phi_{j}$, and

[TABLE]

is the dimensionless repulsive force. The hydrodynamic coupling is represented by the second and third terms on the RHS of (26). Their physical meanings are illustrated in Fig. 1 and are understood as follows. If two particles are moving in tandem, the particle in the rear pushes the front particle via the fluid and the front one pulls the rear. This effect is shown by the second term (tangential-tangential coupling) and increases the phase velocity of the pair. However, because the particles are moving on a circle, the hydrodynamic force exerted by the rear causes a radial displacement of the front particle. The increased orbital radius means decrease of the phase velocity as the driving force is constant. This effect is shown by the third term (tangential-radial coupling). The two terms are explicitly given by

[TABLE]

and

[TABLE]

Note that the RHS of the phase equation (26) are functions of the phase difference only, and written in the form

[TABLE]

Therefore we can get a closed set of equations in terms of the phase difference, as

[TABLE]

In particular, for $N=2$, the phase difference $\Delta_{12}$ obeys

[TABLE]

The latter equality defines the effective potential $V(\Delta)$, which is a useful tool to visualize the stability of the synchronized states [8]. It is calculated by integrating the RHS of (32), as

[TABLE]

II.3 Settings for Numerical Analysis

The model contains three small parameters $\kappa^{-1}=F/(k_{r}R)$, $\alpha=a/R$, and $u=U_{0}/(Fa)$. In numerical analysis, we used $\kappa=30$, $\alpha=0.1$, and $u=10^{-12}$ unless stated otherwise. The first two parameters are in the same order of magnitude as the experimental values $\kappa=30$ and $\alpha\simeq 0.24$ or $0.16$ [15]. The value of $u$ is chosen so that the repulsive force balances with the driving force ($U(\Delta)\simeq Fa$) at the minimum of the effective potential $\Delta=\Delta_{m}\simeq 0.22$. The effective potential is plotted in Fig. 2 with these parameters.

The repulsive potential for $\Delta<\Delta_{m}$ is very steep and prevent the particles to contact each other, while the hydrodynamic coupling causes a long-range attraction.

In the following, we take $\omega^{-1}$ as the unit of time so that the intrinsic phase velocity $\omega$ is unity. We used Mathematica to numerically integrate the phase equations.

III Numerical Results

III.1 The case $N=2$

For an $N=2$ system , we show the time evolution of $\Delta_{12}$ in Fig. 3 with $\kappa$ varied and the initial condition $\Delta_{12}=\pi/2$. The phase difference monotonically decrease to the minimum of the effective potential, where a doublet is formed. As the radial confinement by the optical vortex becomes stronger, it takes more time to form a doublet. The master plot with the scaled time $\kappa^{-1}t$ in the inset of Fig. 3 shows that the timescale of doublet formation is proportional to $\kappa$. This is understood from Eq.(32) and the fact that the repulsive force is negligibly small except near contact.

The time evolution of the radial displacement and phase velocity are shown in Fig. 4. As we see from the figure, the radial displacement of the front particle (particle 2) increases as the two particles approach each other and shows a small drop before it converges to a constant. The displacement of the rear particle (particle 1) is $\delta R_{1}=-\delta R_{2}$ as seen from (25). As the radial displacements are much smaller than the particle radius ($|\delta R_{i}/R|\simeq 0.004\ll\alpha=0.1$), we can justify the assumption that the repulsive force acts only in the tangential direction.

The phase velocities of the two particles monotonically increase and converge to a constant

[TABLE]

which is the velocity of the doublet. Because the velocity of a singlet in the $N=1$ case is set to unity ($v_{s}=1$) by assumption, the ratio between the doublet and singlet velocities is $v_{d}/v_{s}\simeq 1.58$.

A doublet is stable because the rear particle pushes the front particle out from the target trajectory by hydrodynamic interaction, which makes the driving torque smaller for the front particle [15].

III.2 The case $N=3$

The solution of the phase equation (26) for an $N=3$ system becomes a periodic function of time, as shown in Fig. 5. At the time indicated by the left arrow, the particles form a triplet (labeled 1-2-3) as we have $\Delta_{12}=\Delta_{23}\simeq 2\Delta_{m}$. At the time indicated by the right arrow, the particles are separated into a doublet and a singlet as we have $\Delta_{12}\gg\Delta_{23}\simeq\Delta_{m}$. The doublet is faster than a singlet and goes one more cycle than the singlet to form a triplet again. Three collisions take place in a period $T$ to shuffle the triplets, 1-2-3, 3-1-2, and 2-3-1. The period is found to be

[TABLE]

Time evolution of the reduced radial displacements $\delta R_{i}/R$ are shown in Fig. 5 (b). In a triplet, the front (rear) particle has positive (negative) radial displacement, respectively, while the center particle is located on targeted trajectory, as shown schematically in the figure.

The phase velocity of each particle is plotted in Fig. 5 (c). It has three peaks and one deep valley in one cycle. Let us focus on the particle 1 in the plot. The first peak is obtained when the particle is collided at its rear end (by the particle 3) and form a doublet. The second and the highest peak is reached when the doublet collides with a singlet and the particle 1 becomes the center of the triplet 3-1-2. The triplet is split into the doublet 1-2 and the singlet 3 and the phase velocities decrease. The third and the lowest peak is attained when the doublet 1-2 approaches the singlet 3. The minimal phase velocity is marked by a singlet when it is separated from the doublet and is on the opposite side of the circle.

The doublet phase velocity $v_{d}$ is calculated by averaging over the time window during which the particle 1 forms a doublet with the particle 3, as

[TABLE]

The period (35) is estimated as the time for a doublet to go one more cycle than a singlet, multiplied by three to account for the three shuffles, as $T_{N=3}\simeq 2\pi/(v_{d}-v_{s})\times 3\simeq 33.1$, which is only 10% larger than the value (35).

Now let us consider the reason why a triplet is unstable and split [27, 16]. In Fig.5 (d)(e), we show the contributions to $\dot{\Delta}_{ij}$ from the tangential-tangential coupling (the second term on the RHS of Eq.(26)) and the repulsive force, respectively. Compared to them, the contribution of the tangential-radial coupling (the third term on the RHS of Eq.(26)) is negligibly small (not shown). Thus the splitting of a triplet is caused by the tangential-tangential coupling, which is the dominant contribution to $\dot{\Delta}_{ij}$. As seen from Fig. 5(d), the phase differences satisfy $\dot{\Delta}_{12}\sim-\dot{\Delta}_{23}>0$ when the triplet $1-2-3$ is formed. This is because the center particle 2 is hydrodynamically screened and receives a drag force weaker than those on the other particles [27]. Thus the triplet is split into a doublet and a singlet. On the other hand, the repulsive force shows a sharp peak when the particles forming a doublet get closest. However, its contribution to $\dot{\Delta}_{ij}$ is an order of magnitude smaller than the tangential-tangential coupling, as we see from Fig. 5(e).

The phase portrait of the dynamics is given in terms of the trajectory in the $\Delta_{12}$-$\Delta_{23}$ plane, which is shown in Fig.6. In the long-time limit, the trajectory converges to a triangle-shaped limit cycle. In other words, the system exhibits the same rhythmic motion regardless of the initial conditions.

III.3 The case $N=4$

The time evolution of an $N=4$ system is shown in Fig. 7. Depending on the initial configurations, we obtained either (i) a periodically oscillating state consisting of a doublet and two singlets that are shuffled with each other, or (ii) a stable synchronized state consisting of two doublets.

The period for the oscillating state is found to be

[TABLE]

which is much smaller than the period for $N=3$. This is because two singlets are positioned in the opposite side of the circle when a doublet moves between the two. Therefore, the distance traveled by a doublet in an $N=4$ system between two collisions is roughly half of the one for $N=3$. Taking into account the four shuffles per period, and using the doublet velocity $v_{d}=1.63$, we get the estimate $T_{N=4}\simeq\pi/(v_{d}-v_{s})\times 4\simeq 19.9$.

For the stable synchronized states, the phase difference between two doublets is $\Delta_{23}=\Delta_{41}\simeq 0.92\pi$, which means that the doublets are placed in the opposite side of the circle.

In Fig. 8, we show the state diagram in terms of the initial conditions. The initial conditions are varied by setting $\phi_{i}$’s $(i=1,2,3,4)$ as integer multiples of $\epsilon=0.1\pi$. The state diagram is drawn in the $(\Delta_{12},\Delta_{23},\Delta_{34})$ space with $\Delta_{ij}$’s ranging in $[\epsilon,16\epsilon]$. It shows that the system reaches the stable state if two of the particles are close to each other in the initial state. For example, the corner points of the state diagram, $(\Delta_{12},\Delta_{23},\Delta_{34})=(16\epsilon,\epsilon,\epsilon),(\epsilon,16\epsilon,\epsilon),(\epsilon,\epsilon,16\epsilon)$, give the stable state. The oscillating state is obtained if the particles are equally spaced or form a triplet and a singlet in the initial state. The number fraction of the initial configurations that give rise to the oscillating state was 54%. Therefore, it is slightly more probable to get the oscillating state if we randomly choose the initial configuration.

III.4 The case $N=5$

In Fig. 9, we show the time evolution of the phase difference, radial displacements and phase velocities for an $N=5$ system. As in the case $N=3$, the system converges to a periodic oscillating state independent of initial conditions. The value of the period is almost the same as in the case $N=3$:

[TABLE]

The phase difference exhibit double peaks in a period, which is a feature not found in the cases $N=3,4$. This is understood as follows. At the first peak marked by the arrow a) in the figure, the triplet 4-5-1 and the doublet 2-3 are formed. Then the triplet moves faster than the doublet as seen in Fig. 9 (middle), because the center particle is hydrodynamically screened [27, 16]. The triplet is decomposed into a doublet and a singlet at the time indicated by the arrow b). The new doublet (5-1) is slower than the other doublet (2-3) because it has a larger radial displacement, as seen from Fig. 9 (bottom). This causes the slight decrease of the phase difference $\Delta_{12}$ and the dip between the two peaks. Then the triplet 2-3-4 is formed at the time indicated by the arrow c), and $\Delta_{12}$ increases again to second peak.

III.5 The case $N=6$

Time evolution of the phase differences for an $N=6$ system is shown in Fig.10. As in the case $N=4$, the system converges to either an oscillating or a stable synchronized states, depending on the initial conditions. The oscillating state consists of two doublets and two singlets. The period is almost the same as in the case $N=4$:

[TABLE]

The stable state consists of three doublets that are equally spaced.

IV Discussion

Now let us discuss the results in comparison with the previous studies and consider their physical meanings.A quantity of interest is the ratio between the phase velocities of a doublet and a singlet, $v_{d}/v_{s}$. It is always larger than unity because the hydrodynamic resistance is reduced as two particles get close. For an $N=2$ system, we can estimate the ratio analytically by considering two spheres that are in contact with each other and moving in tandem along a straight line. The Rotne-Prager-Yamakawa approximation (12) gives $v_{d}/v_{s}=8/5=1.6$, which is close to our numerical result $v_{d}/v_{s}\simeq 1.58$. The small discrepancy is resulting from the non-vanishing radial displacement and tangential separation in the numerical analysis. The Oseen tensor, which is obtained by neglecting the correction terms in the Rotne-Prager-Yamakawa tensor due to finite particle radius, gives the much smaller value $v_{d}/v_{s}=4/3\simeq 1.33$. A more exact value is obtained from the lubrication theory [31] as $v_{d}/v_{s}\simeq 1.578$. This is fairly close to the value from the Rotne-Prager-Yamakawa approximation as well as to our numerical result. The finite temperature simulations by Sokolov et al. [15] gave $v_{d}/v_{s}=1.57\pm 0.05$, which is also close to our result. Experimentally, the ratio $v_{d}/v_{s}$ is reported as $1.10\pm 0.02$ in Ref.[15] and $1.32$ by Sassa et al. [16]. The reason why the experimental values are smaller than the theoretical ones is not clear at this stage. However, it is reported in Ref. [16] that fluctuations in the singlet velocity is as large as 27%, which is attributed to nonuniformity of the optical driving force along the circular path. In an earlier experiment [14], the fluctuations of the optical intensity along the trajectory is also about 20% from the mean. If fluctuations in the radial displacement is large, they would cause a larger viscous drag force on a doublet because the hydrodynamic screening by the front particle is weakened.

For $N\geq 3$, we confirmed the existence of oscillating modes with mutually exchanging doublets and singlets. The limit cycle is robust and obtained from any initial conditions. This is in tune with the experimental finding that a triplet collapses in some occasions but is soon recovered [16].

For $N=3$, the oscillating state is described by a limit cycle in the phase portrait. Similar phase portraits are obtained in the experiment [18] and in the direct numerical simulation [17]. However, the trajectories in the experiment show noisy fluctuations, which are probably due to the fluctuations of the driving force. The trajectories obtained from the simulation also show oscillations, which should be due to periodic variation of the driving force assumed in the model.

For $N=4$ and $N=6$, both oscillating and stable states were observed in the experiment [16] and in the simulation [15], in agreement with our results. According to our analysis, one of the two states is finally chosen depending on the initial conditions. However, in the experiment, transition between the two states in a time course was observed, and finally only the stable states remained. This transition may be also attributed to fluctuations of the optical driving force [16]. There are not only quenched heterogeneities but also temporal fluctuations of the optical force. In the simulation with thermal noise [15], the oscillating mode for $N=4$ disappeared at temperature higher than the experimental temperature. Temporal fluctuations of the optical force may be described by high effective temperature and explain the transition to the stable mode.

For $N=5$, we found an oscillating state with two doublets and a singlet, in accordance with experiments [16]. In the experiment, another oscillating mode with one stable doublet and three singlets is also observed sometimes and in a short time period. An oscillating mode with one doublet is also observed for $N=6$. These modes containing one doublet is not reproduced in our analysis, and is possibly due to fluctuations of the optical driving force.

In Table 1, we show the maximal, minimal and time-averaged phase velocities as well as the periods of the oscillating modes, all obtained in the dynamical steady states. The phase velocities for the stable modes ($N=2,4$ (stable), $6$ (stable) in Table 1) are those of a doublet. For oscillating modes for $N=3,4$ (osc.), $5,6$ (osc.), the maximal phase velocity is realized by the front particle of a doublet when the doublet collides with a singlet and a triplet is formed. The minimal phase velocity is realized by a singlet when it is separated from other particles. For example, for $N=3$, the singlet is slowest when the doublet is in the opposite side of the circle, because the doublet causes a hydrodynamic drag force that is anti-parallel to the moving direction of the singlet. Therefore the deviation of the minimal velocity from the reference value (unity) is in the order of $\alpha$.

The mean phase velocity can be compared with those reported in the preceding studies [15, 16]. The mean phase velocity for an $N=2n$ system in the stable mode is larger than that for $N=2n-1$ ($n=1,2,3$) . This even-odd effect is reported in the experiments [15, 16] and is interpreted by the number of doublets: for $N=2n$, The mean phase velocity for $N=2$ is larger than that for $N=4$ (stable). This also agrees with the experimental results [15, 16]. We interpret it as effect of the hydrodynamic drag force between the doublet and the two singlets that are moving in the opposite side of the circle, just as in the above explanation of the minimal singlet velocity. On the other hand, the mean phase velocity for an $N=4$ system in the oscillating mode is smaller than that for $N=3$, which is consistent with the previous simulation result [15]. As a function of $N$, the mean phase velocity took minimum value at $N=3$ in the experiments [15, 16]. However, we found that the minimum is achieved by an $N=4$ system in the oscillating mode, which is not stable in the experiments as mentioned above.

In conclusion, we analyzed theoretically the collective dynamics of colloidal particles driven by an optical vortex. We derived the phase equation for $N$ particles and quantitatively investigated the mechanisms behind the rich dynamical behaviors, such as: (i) the radial degree of freedom for doublet formation, (ii) hydrodynamic screening for triplet splitting, and (iii) symmetry of particle configurations determining the basin of attractors. The results are in good agreement with the experimental findings, including the even-odd effects in the phase velocity. A remaining problem is an understanding of the transitions between the oscillating and stable modes observed in the experiments, which might be related to spatio-temporal fluctuations of the driving force. In future experiments, a more precise control of the fluctuations of the optical force will be useful in unveiling the mechanism of the transitions, and also in quantitative assessment of the phase velocities. Theoretically, both spatial and temporal fluctuations of the driving force will be incorporated using the existing framework [8, 2]. The effects of confining walls [28] and polydispersity of the particle size [17] also remain to be analyzed in future work.