2 FIG. 1: Schematic picture of the rotor (top view). The active force F$_{i}$ exerted by the i -th rotor on the fluid is deviated by a fixed angle δ i from the radial direction. The reaction force drives the rotor at the intrinsic frequency ω i = F sin δ i /ζb , where b is the radius of rotation and ζ the viscous drag coef- ficient of the bead. F x y i i b i (a) (b) tally (flagellum). Motion of the bead is constrained on a circular orbit of radius b located at height h from the substrate, which we take to be the xy -plane. The posi- tion of the i -th bead is thus given by r i = r$_{0}$ $_{i}$+ h e $_{z}$+ b n i where r$_{0}$ i is its base position on the square lattice and n i = (cos φ$_{i}$, sin φ$_{i}$, 0) is the unit vector that gives the ori- entation of the arm via its phase φ$_{i}$ = φ$_{i}$ ( t ). The velocity of the bead reads v i = ˙ φ$_{i}$ t i where t i = ( − sin φ$_{i}$, cos φ$_{i}$, 0) is the unit vector tangential to the trajectory. We as- sume that the active force F i exerted by the rotor on the fluid has a constant magnitude F , and makes a fixed angle δ$_{i}$ (measured clockwise) from the radial direction; F i = F (cos δ$_{i}$ n i − sin δ$_{i}$ t $_{i}$). See Fig. 1 for the configura- tion. The reaction force − F i on the rotor arm gives the driving torque T$_{i}$ = F b sin δ$_{i}$ and the intrinsic frequency ω$_{i}$ = F sin δ$_{i}$/ζb , where ζ = 6 πηa is the viscous drag co- efficient. The mounting angles δ$_{i}$ ’s are assumed to have the Gaussian distribution with the standard deviation δ$_{0}$ . We assume that the rotors are widely spaced so that a, b, h ≪ d . Then the velocity field of the fluid created by the active forces is given by v ( r ) = ∑ i G ( r − r $_{i}$) · F $_{i}$, where G ( r ) = (3 h $^{2}$/ 2 πη ) · r ⊥ r $_{⊥}$/ | r | $^{5}$, r ⊥ = ( x, y, 0) is the asymptotic expression of the Oseen-Blake tensor [23] for h/d ≪ 1. The rotor’s angular velocity is given by ω$_{i}$ + v ( r $_{i}$) · t $_{i}$/b , or, more explicitly, Here, ω$_{0}$ = F/ζb , r ij = r i − r j = | r ij | (cos θ$_{ij}$, sin θ$_{ij}$ ), and γ = ζh $^{2}$/ηd3 = 6 πah $^{2}$/d3 is the dimensionless coupling constant. For a real bacterial carpet, a ∼ h ∼ 1 µ m and d ∼ 10 µ m give the rough estimate γ ∼ 10 − $^{2}$. Numerical Simulation. - We implemented the model on a L × L square lattice and numerically integrated eq. (2) by the Euler method. We assumed the peri- odic boundary condition and computed the velocity field every time step in the Fourier space. We set γ = 0 . 1 and varied the angle deviation δ$_{0}$ as the control parameter. The system size used was L = 128 for most of the results shown below, while L = 32, 64 and 256 are also used to check finite-size effect. Starting from random initial con- figurations of φ$_{i}$ ( t = 0), the system reached a dynamical steady state by the time t = 1 × 10 $^{4}$/ω$_{0}$ . The statisti- cal data shown below are taken from the time window 1 × 10 4 < ω$_{0}$t < 2 . 5 × 10 $^{5}$. We plot the orientational order parameter S = |〈 n 〉| = |〈 e $^{iφ}$〉 | as a function of δ$_{0}$ in Fig. 2. Also shown is the standard deviation (STD) of the actual frequency Ω$_{i}$ = 〈 ˙ φ$_{i}$ 〉 normalized by the STD of the intrinsic frequency ω$_{i}$ , Q = √ 〈 Ω $^{2}$〉 / 〈 ω $^{2}$〉 . Note that S = 1 and Q = 0 in the fully synchronized state and S = 0 and Q = 1 in the desynchronized limit. As we increase δ$_{0}$ , S and Q slowly converge to the desynchronized limit. While the change in the orientational order parameter is sharper for a larger system size, the frequency deviation has little L - dependence. For comparison, we also show the results of the mean-field theory, which will be explained in the next section. In Fig. 3, we plot the distribution function of the ac- tual frequency Ω normalized by the STD of intrinsic fre- quency, for different values of δ$_{0}$ . The distribution con- sists of a sharp delta-function like peak at Ω = 0 and broad symmetric tails for Ω > 0 and Ω < 0. For δ$_{0}$ ≤ 3 $^{◦}$, most of the rotors are coherent and contribute to the center-peak. For δ$_{0}$ = 10 $^{◦}$, the distribution is close to that of the intrinsic frequency, while the center peak still remains. the above data suggest that the synchroniza- tion transition in this system is more gradual than that found in a globally coupled system, and it is difficult to locate the transition point exactly. On the other hand, the variance of the order parameter Var( S ) = 〈 S ( t ) $^{2}$〉 − 〈 S ( t ) 〉 2 as a function of δ$_{0}$ (Fig. 4(a)) has a peak near δ$_{0}$ = 6 $^{◦}$, suggesting that there is a subtle balance between synchronization and desynchronization. We will call this the threshold angle and denote by δ$_{th}$ . In Fig.4(b), we plot the temporal correlation function of the order parameter C$_{S}$ ( t ) = 〈 S ( t + t $^{′}$) S ( t $^{′}$) 〉−〈 S 〉 $^{2}$. We find an oscillatory behavior with long correlation time at δ$_{0}$ = 6 $^{◦}$. The presence of a characteristic period is also di- rectly observed in the plot of S ( t ) in Fig. 4(c). Although the origin of the oscillation is beyond the scope of the present paper, a preliminary study shows that the oscil- latory behavior at the threshold angle is a unique feature resulting from the long-ranged nature of the interaction. We also plot the orientational correlation function in Fig.5. Here the angular brackets mean taking average over r ′ as well as the azimuthal angle of r . For 4 ≤ δ$_{0}$ ≤ 6 $^{◦}$, we observe an exponential decay of the correlation G$_{n}$ ( | r | ) = 〈 [ n ( r + r $^{′}$) − n ] · [ n ( r $^{′}$) − n ] 〉 (3) dφ$_{i}$ dt = ω$_{0}$ sin δ$_{i}$ − 3 γω$_{0}$d 3 4 π ∑ j = i 1 | r ij | 3 [ sin( φ$_{i}$ − φ$_{j}$ + δ$_{j}$ ) + cos( φ$_{i}$ + φ$_{j}$ − δ$_{j}$ − 2 θ$_{ij}$ ) ] . (2) P ( δ ) = 1 √$_{2}$$_{πδ}$$_{0}$ exp ( − δ 2 2 δ 2 0 ) (1)