E ﬀ ects of thermal ﬂ uctuations and ﬂ uid compressibility on hydrodynamic synchronization of microrotors at ﬁ nite oscillatory Reynolds number: a multiparticle collision dynamics simulation study

We investigate the emergent dynamical behavior of hydrodynamically coupled microrotors by means of multiparticle collision dynamics (MPC) simulations. The two rotors are con ﬁ ned in a plane and move along circles driven by active forces. Comparing simulations to theoretical results based on linearized hydrodynamics, we demonstrate that time-dependent hydrodynamic interactions lead to synchronization of the rotational motion. Thermal noise implies large ﬂ uctuations of the phase-angle di ﬀ erence between the rotors, but synchronization prevails and the ensemble-averaged time dependence of the phase-angle di ﬀ erence agrees well with analytical predictions. Moreover, we demonstrate that compressibility e ﬀ ects lead to longer synchronization times. In addition, the relevance of the inertia terms of the Navier – Stokes equation are discussed, speci ﬁ cally the linear unsteady acceleration term characterized by the oscillatory Reynolds number Re T . We illustrate the continuous breakdown of synchronization with the Reynolds number Re T , in analogy to the continuous breakdown of the scallop theorem with decreasing Reynolds number.


Introduction
Cell motility is a remarkable accomplishment of evolution, being fundamental for a variety of cellular activities.2][3] Furthermore, agellar motion plays a major role in eukaryotes, 1,2 where they transport uid in the respiratory system in form of cilia, 4 are involved in cellular communications, 5 and even determine the morphological leright asymmetry in the embryo. 6he various types of collective agella and cilia motions require a synchronization of their beating pattern.][22][23][24][25][26][27][28] The synchronized beating of nearby swimming spermatozoa was rst modeled and analyzed in ref. 29, suggesting that hydrodynamic interactions lead to synchronization.Since then, the hydrodynamics of active systems at low Reynolds numbers has become a subject of major interest. 3,30Indeed, recent experimental results on the uid-mediated interactions between agella of two Volvox carteri cells unambiguously reveal the importance of hydrodynamic interactions on their synchronized motion. 31However, synchronization is not easily achieved for systems at zero Reynolds number, 15,16,30,[32][33][34] where the uid dynamics is described by Stokes equations.The presence of kinematic reversibility of these equations combined with swimmer symmetries may prevent synchronization. 15,16,30,32,34To overcome this fundamental limitation of Stokes equations and to generate a non-reciprocal dynamics, various alternatives have been suggested.This comprises inclusion of additional degrees of freedom such as system exibility, 7,[15][16][17][18]20,23,25,33,35,36 a specic coupling, 37 e.g., via surfaces, 38 or nonreversible driving forces. 13,36,39,40 In adition, specic system designs combined with hydrodynamic interactions lead to synchronization as has been shown for sheets 32,41 and models of Chlamydomonas.][43][44] Alternatively and more general, the linear unsteady Stokes equation can be adopted to describe the uid properties.34,45 Here, the time reversibility of the underlying dynamical equation is broken.
The described lack of synchronization due to kinematic reversibility is equivalent to Purcell's scallop theorem 46 for locomotion at zero Reynolds numbers.The same type of symmetry arguments apply for both effects.Hence, the two examples belong to a more general class of physical phenomena, which are all based on the specicities of the time independent, linear, and zero Reynolds number Stokes equations, which are identical under time reversal.Hence, any insight into a particular phenomenon has implications for a broader range of physical effects.Recently, it has been shown 45,47 that even for oscillatory reciprocal forcing of a solid body at arbitrary small Reynolds numbers Re > 0 a net translational motion is obtained.6][47] The consideration demonstrates, however, that the scallop theorem breaks down in a continuous way and no critical inertia is required for the onset of motion. 47Implications for other phenomena, e.g., synchronization need to be resolved.
For a theoretical understanding of the mechanisms governing synchronization in complex (biological) systems, computer simulations are particularly benecial, since they capture potential non-linear effects, hydrodynamic interactions, and thermal uctuations. 17,18,48,49The latter have been recently shown to markedly effect the dynamics of E. coli bacteria. 50Thereby, mesoscale hydrodynamic simulation techniques such as the multiparticle collision dynamics (MPC) method [51][52][53][54] are particularly valuable.][59][60][61][62][63][64][65][66][67][68] In this article, we apply the MPC approach to study the time-dependent, correlated dynamics of two hydrodynamically coupled microrotors.We adopt the rotor model proposed in ref. 23.As we have shown analytically, the rotational dynamics of such coupled rotors is synchronized by timedependent hydrodynamic interactions. 34MPC simulations are performed to verify the previous results, to elucidate the inuence of thermal uctuations and that of uid compressibility.Moreover, we discuss the importance of the inertia terms of the Navier-Stokes equation-the linear unsteady acceleration term and the non-linear advection term-by introducing corresponding Reynolds numbers 45,47 and demonstrate that the two can be comparable for micrometersize objects.From a simulation technical point of view, we establish a link between the Landau-Lifshitz Navier-Stokes uctuating hydrodynamic equations and the MPC approach in terms of stress tensors for angular and non-angular momentum conserving simulations.Specically, we present the respective dissipation uctuation relations.
The paper is organized as follows.The uid and the rotor model are described in Sec. 2. In Sec. 3, relevant basic hydrodynamic equations are discussed.Theoretical results for synchronization and the ow eld are presented in Sec. 4, and simulation results are shown in Sec. 5. Finally, Sec. 6 summarizes our ndings.
2 Fluid and rotor model

Multiparticle collision dynamics
The MPC uid consists of N point particles of mass m with the positions r i and velocities v i (i ¼ 1, ., N), which interact with each other by a stochastic, momentum-conserving process.The particle dynamics proceed in two steps-streaming and collision.In the ballistic streaming step, the particle positions are updated via where h is the collision time step.In the collision step the simulation box is partitioned into cubic cells of length a, in which multiparticle collisions are performed.In the SRD version of MPC, [51][52][53] the relative velocity of each particle, with respect to the center-of-mass velocity of the cell, is rotated by a xed angle a around a randomly oriented axis.Hence, the new velocities are where R(a) is the rotation matrix, 69 is the center-of-mass velocity, and N c the total number of particles in the cell.In its original version, MPC violates Galilean invariance.It is restored by a random shi of the collision grid at every step. 70he collision rule (2) violates angular momentum conservation, 71 which is associated with a non-symmetric stress tensor. 54,72,73Angular momentum conservation is reestablished on the cell level by a solid-body type rotation of relative velocities aer a collision according to 71 r j;c ðtÞ Â ½v j;c ðtÞ À RðaÞv j;c ðtÞ É # where I is the moment of inertia tensor of the particles in the center-of-mass reference frame; r i,c (t) and v i,c (t) are the respective relative positions and velocities aer streaming, i.e., r i,c ¼ r i À r cm and v i,c ¼ v i À v cm , with the center-of-mass position r cm .
In order to simulate an isothermal uid, a collision-cellbased, local Maxwellian thermostat is applied, where the relative velocities of the particles in a cell are scaled according to the Maxwell-Boltzmann scaling (MBS) method. 74

Microrotor
We adopt the rotor model of ref. 23: two beads of radius R H move along circles of radius R, each driven by an active force F i (cf.Fig. 1).The two circles are centered at , where êx is the unit vector along the x-axis and d the centerto-center distance; both beads are conned in the xy-plane.The trajectories of the bead centers can be expressed as in terms of the phase angles 4 i (t).The driving forces are of equal magnitude and point along the tangents ti (t) of the trajectories, where 6][77][78] Thereby, the beads are treated as uid particles with the mass M. Thus, in cells with a bead the center-of-mass velocity is given by v cm ðtÞ ¼ The velocity V(t + h) of the bead aer a collision follows according to eqn (2) or (4), respectively, taking into account the appropriate mass M.
Hence, the phase angles evolve as between MPC collisions.In a collision, the angular velocities change.The new values aer a collision follow from the relation

Parameters
All simulations are performed with the rotation angle a ¼ 130 , the mean number of particles per collision cell hN c i ¼ All results presented in the following have been obtained with the non-angular momentum conserving collision rule (2).However, simulations conrm that these results are independent of the applied MPC variant and they agree very well with each other.

Hydrodynamics
As is well established, the hydrodynamic properties of the MPC uid are excellently described by the linearized Navier-Stokes equations on sufficiently large length and time scales. 51,52,54,72,73oreover, the dynamics of the coupled rotors of Sec.2.2 can be described analytically within the linearized Navier-Stokes equations. 34In order to establish a link between simulation results and analytical considerations, we will briey describe the basic hydrodynamic background.

Navier-Stokes equations and MPC uid
Mass and momentum conservation of an isothermal uid are expressed by the continuity and the Navier-Stokes equation [81][82][83] respectively.Here, r + dr(r, t) denotes the mass density with its mean value r and its (small) uctuations dr(r, t) at the position r in space and the time t.v ¼ v(r, t) is the uid velocity eld, f(r, t) a volume force, and s(r, t) the stress tensor.Note that only the mean density r will appear throughout the rest of the paper.
For the MPC uid, the explicit form of the stress tensor depends on the lack or presence of angular momentum conservation during collision.In general, the stress tensor can be expressed as with the pressure p and the Cartesian indices a, b, a 0 , b 0 {x, y, z}.
(i) For angular momentum conserving uids, the standard symmetric stress tensor is obtained 81 with Fig. 1 Model system of hydrodynamically coupled rotors. 23Two beads move along fixed circular trajectories, each driven by a constant tangential force F.
in three dimensions.h denotes the shear viscosity.Note, we assume that the bulk viscosity is zero. 54,84ii) For the non-angular momentum conserving variant of the MPC approach of Sec.2.1, the non-symmetric tensor is given by 73,84,85 with 86,87 Alternative expressions for the stress tensor have been used, which differ from eqn ( 16) by a term with vanishing divergence only and thus yield the same Navier-Stokes equation. 54,71,72,84ith the stress tensor (16), eqn (12) turns into For angular momentum conserving uids h k is replaced by h in eqn (18).

Inertia and Reynolds numbers
In order to asses the relevance of the various terms in eqn (18), in particular the time-dependent and non-linear inertia terms, we scale the velocity eld by a typical value V, length by L, and time by T, as usual, 88 which yields the equation where the primed quantities are dimensionless and of O (1).Furthermore, we introduced the Reynolds numbers with the kinematic viscosity n ¼ h/r.Typically, T is dened as T ¼ L/V which yields Re T ¼ Re, and for Re ( 1 the le hand side of eqn ( 19) is neglected.In particular, Re ¼ 0 is assumed for microswimmers, which leads to peculiarities in their locomotion as expressed by the well-known scallop theorem. 46he oscillatory Reynolds number Re T can be written as Re T ¼ s n /T, with s n ¼ L 2 /n.Hence, Re T is the ratio of the viscose time scale s n for shear wave propagation over the distance L and the characteristic system time T. In order to establish proper hydrodynamic interactions, s n /T < 1 and, hence, Re T < 1.
In ref. 45 and 47, it is assumed that Re T ( Re, which may apply for the systems considered there, but may not be valid in general.Alternatively to the proposed breakdown of the scallop theorem due to the nonlinear advection term, the linear unsteady acceleration term implies such a breakdown even at zero Reynolds number Re, because the term breaks time reversibility.
A priori, the relevance of the various Reynolds numbers is not evident.As far as the rotor model of Fig. 1 is concerned, the relevant length scale is the distance between the two beads L ¼ d, the characteristic time is the rotational period, and the relevant velocity is V ¼ uR z R/T, with the rotational frequency u.Hence, Re/Re T z R/d ( 1, since typically R ( d.Thus, Re T is at least as important as Re itself.
For a quantitative estimation of the Reynolds number, we consider the beating of the agella of Chlamydomonas. 39,42With L z 15 mm, R z 4 mm, and T z 15 ms, we nd Re z 10 À3 and Re T z 10 À2 in water.
The emergent dynamical behavior due to the presence of the unsteady acceleration term has recently been studied experimentally and theoretically for a colloidal particle conned in an optical trap. 89Here, the characteristic time is the trap relaxation time s K ¼ 6phR H /K, which depends on the trap stiffness K.The studies reveal strong hydrodynamic correlations for Re T ¼ s n /s K T 0.1.More precisely, the studies show the self-interaction of a single colloid mediated by its surrounding uid.This fundamentally changes the particle's thermal uctuations, which are mediated by hydrodynamic self-interactions. 90A similar type of coupling can be achieved for two simple colloidal rotors.This might open opportunities to exploit hydrodynamic phenomena for novel sensing applications.

Fluctuating hydrodynamics
In the following, we focus on the linearized Navier-Stokes equation and neglect the advection term, i.e., we set Re ¼ 0, but keep the unsteady acceleration term.Moreover, to account for thermal uctuations inherent in the MPC uid, we introduce a random stress tensor s R (r, t).Then, eqn (18) turns into (Landau-Lifshitz Navier-Stokes equation) 54,56 The volume force f now consists of a deterministic force f D and the random force f R ¼ V$s R .The stress tensor s R is assumed to be a Gaussian and Markovian stochastic process with the moments 54 h aba 0 b 0 is either given by eqn (15) or (17) depending on the applied MPC procedure.
Taking the divergence of eqn (22) and using the linearized continuity eqn (11) together with the ideal gas equation of state with the velocity of sound c ¼ ffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi k B Q=m p for an isothermal system, we obtain the wave equation In Fourier space, eqn ( 22) and ( 24) read as 54,56 This

View Article Online
with k ¼ |k| and the argument (k, u) of all variables with a bar.In order to solve eqn ( 25) and ( 26) we dene the longitudinal and transverse projection operators P L (k) ¼ kk , with the unit vector k and the dyadic product kk , and P T (k) ¼ 1 À P L (k) along with v L ¼ P L v and v T ¼ P T v. Using these projection operators, we obtain the solution where and with h ¼ h + h k /3.The transverse part of the velocity v T describes shear waves, whereas the longitudinal part v L represents sound waves.Fourier transformation with respect to frequency yields the time-dependent tensors 54 where For the longitudinal part, we nd " for 4c 2 > k 2 ñ2 , and " , and assume t $ 0. We obtain the velocity eld v(r, t) by Fourier transformation of eqn (27) and using the convolution theorem vðr; tÞ ¼ The Green's function have been derived. 34,91For the longitudinal part, we obtain the integral representation A L ¼ I 2 and B L ¼ 3I 2 À 2I 1 , with 4 Theoretical results

Rotor equations of motion
We consider the limit of point particles and assume a balance between the driving and the hydrodynamic friction forces acting on the beads, i.e., where g ¼ 6phR H is the friction coefficient and v(r, t) is given by eqn (33).The forces F i and F R i are related with the respective force densities f D and f R according to with the volume V i of particle i.Hence, we obtain the equations of motion The restriction to a circular trajectory of a bead implies a constraining force.In the following considerations, we neglect this constraining force, as it leads to terms of quadratic order in R H /d.Moreover, for an analytical solution of the equations of motion (40), we apply the mean-eld approximation R 2 (t) À R 1 (t 0 ) z d ¼ dê x , which is strictly valid for d/R [ 1, 34 and neglect the random force, i.e., we set F R i h 0.
Using the position vectors ( 5) and forces (6), eqn (40) yield the coupled integro-differential equations for the phase angles 4 i .Here, u ¼ F/(gR) ¼ 2p/T is the intrinsic angular frequency with the rotation period T. 34 The creeping-ow limit studied previously 23,30,36,40 follows, when we assume that tj (t 0 ) under the integral is changing much more slowly in time than the hydrodynamic tensor Q.Then, integration of the transverse hydrodynamic tensor Q T (r, t) yields the asymptotic expression for t / N ð t with the unit vector r ¼ r/|r| and the dyadic product rr.The rst term on the right hand side is the well-known Ossen tensor. 88,92he long-range hydrodynamic correlations yield a leading timedependent term, which decays slowly as $ 1= ffiffi t p .The sound contribution decays much faster, 54 since it includes an exponentially decaying factor and is thus irrelevant for the asymptotic behavior.

Synchronization and ow eld
To characterize the rotor dynamics, we solve eqn (41) numerically by a combination of Euler's integration scheme and the trapezoidal rule. 93In eqn (41), Q(r, t) ¼ Q L (r, t) + Q T (r, t) is required in position space.In contrast to the transversal component Q T , for which a closed expression is given in eqn ( 35) and ( 36), the longitudinal component Q L has to be computed numerically for each time step by evaluating the integrals (37).
Fig. 2 shows numerically obtained phase-angle differences D(t) ¼ 4 1 (t) À 4 2 (t) for various initial differences D(0).The phase-angle differences exhibit an exponential decay for long times, except for the initial value D(0) ¼ p, which is a metastable state as discussed in ref. 34.Hence, the phase difference vanishes with time and the rotors exhibit an inphase synchronized rotation.
As discussed in ref. 34, the equation of motion for the phase difference is well approximated by for t [ s n ¼ d 2 /n.Hence, D(t) decays exponentially in the asymptotic limit t / N with the characteristic time Neglecting compressibility effects by setting The analytical result ( 45) agrees very well with numerically determined characteristic decay times for the full numerical solution of eqn (41), with Q ¼ Q T , as long as T [ s n . 34he synchronized rotation is energetically favorable, as can be visualized by the velocity eld.For pointlike particles, it can be computed according to Fig. 3 depicts examples for an asynchronous and synchronous state.The velocity prole decays very fast with separation from a beat.Thereby, it is anisotropic-of Stokeslet shapewith a slower decay in the tangential forward and backward direction.In the synchronous state, the overall ow eld is smoother, which minimizes dissipation.

Discussion
The linear unsteady acceleration term leads to synchronization of the rotational motion of the rotors (cf.Sec.3.2 and 4.2).Interestingly, the ratio of the synchronization time and the rotational period t s /T is proportional to the inverse of the square root of the Reynolds number Re T ¼ s n /T.This implies that synchronization appears for any non-zero Reynolds number Re T .However, a vanishing Reynolds number leads to an innitely long synchronization time.This is in agreement with the fact that there is no synchronization in the Stokes limit. 23ence, we get a continuous crossover from a synchronizing to a non-synchronizing system.This is similar to the considerations in ref. 45 and 47 for the non-linear inertia term, with the continuous breakdown of the scallop theorem with inertia due to the unsteady acceleration term.
The dependence on the radio d/R H is characteristic for synchronization and has been found for other synchronization mechanisms. 25,36,40Interestingly, in our case t s can solely be expressed by the ratio of the viscous time scale R H 2 /n and the period T. Hence, synchronization is determined by shear wave propagation over the size R H of the sphere rather than the distance between the two rotors.However, we have to keep in mind that the approximate expression (45) applies for T [ s n only, i.e., shear waves traversed the distance between the beads.In the limit t 0 [ s n , the integrand of eqn (44) reduces to (4pnt 0 ) À3/2 sin(ut 0 ), and the integral is determined by the uid long-time tail only. 54The details of the set-up become irrelevant on time-scales longer than the viscous time-scale s n .This is similar to the polymer center-of-mass velocity auto-correlation function, which is also independent of polymer properties on respective long time scales and is solely determined by the uid long-time tail. 56On short time scales (t/T ( 10 in Fig. 2), however, there is naturally a signicant dependence of synchronization on d due to the nite shear-wave propagation time.
For a quantitative estimation of the synchronization time of a biological system, we return to Chlamydomonas.In Sec.3.2, we estimated the Reynolds number to be Re T z 10 À2 .Taking as radius of a agellum R H z 0.1 mm, we nd t s /T z 50, i.e., 50 periods would be required to synchronize the beating between the two agella.Hence, the provided mechanism is denitely too slow to be effective for synchronization of Chlamydomonas agella.An important factor here is the thin agellum.
Synthetic systems, such as colloids driven along circular trajectories, may provide more control over the relevant parameters.Consider two colloids dissolved in water with a radius of R H ¼ 2 mm driven along a circle with the rotation period T ¼ 10 À2 s.The distance between the rotors and the radius of their trajectories does not really matter.We will assume values adopted in ref. 94, i.e., R z 3 mm and d z 20 mm, which yields t s /T z 2 and s n /T ¼ d 2 /(nT) ( 1. Hence, within about two cycles the colloids synchronize their rotation.As indicated in ref. 90 for colloids conned in optical traps, we expect that the rotor setup can be exploited for novel sensing applications for colloidal particles of various shape or living cells due to hydrodynamic phenomena.

Co-rotating beads
MPC simulation results for the time dependence of the phaseangle difference D(t) ¼ 4 1 (t) À 4 2 (t) of co-rotating beads are displayed in Fig. 4. The three phase-angle differences decay   initially, but there is no clear asymptotic value.In contrast, the variations due to thermal uctuations are rather large.To extract the average behavior, we performed up to 10 3 independent simulation runs for a parameter set.The average D indeed shows synchronization of the bead rotation despite strong thermal uctuations as displayed in Fig. 5 for various Péclet numbers.As is evident from the gure, the simulation results compare well with the theoretical predictions according to eqn (41) when taking into account compressibility effects, although there are small differences for Pe ¼ 140 and 160.
The distribution of phase-angle differences in the synchronized state are shown in Fig. 6 in a polar representation.The simulation data are well described by a wrapped normal distribution, where the mean m and standard deviation s of a set of phase differences D n are calculated according to for N r realizations. 95he time dependence of the standard deviations s for the systems of Fig. 5 are presented in Fig. 7.We nd a signicant dependence of s on the applied force.Thereby, the uctuations decrease with increasing Pe.For the considered range of Péclet numbers, the ratio of the standard deviations (z0.4/0.2 ¼ 2) decreases faster with increasing Pe than the ratio of the Pe decreases (120/180 ¼ 2/3).The increase of the uctuations with decreasing Pe implies very large s on the order of p.
We like to stress that synchronization reveals a strong dependence on periodic boundary conditions, as illustrated in Fig. 8 and already reported in ref. 25 for a similar, but timeindependent rotor model.The synchronization time is evidently signicantly shorter for smaller simulation boxes, increases with increasing box size L x and approaches an asymptotic value.
Within the accuracy of the simulations, the curves for L x ¼ 100a and 150a are similar and are therefore close to the asymptotic limit.The following results are obtained for the box length L x ¼ L y ¼ 100a and L z ¼ 20a.
The inuence of compressibility on the decay of the phaseangle difference is illustrated in Fig. 9, where simulations are compared with analytical results for incompressible and compressible uids.Interestingly, for the selected parameters, compressibility slows down synchronization considerably.
In order to quantify the effects of compressibility, we dene the maximum deviation of the phase-angle differences and the difference of the synchronization times between a compressible (index c) and an incompressible (index i) uid.Results for these quantities obtained from the theoretical description of Sec.IV are presented in Fig. 10 for various Péclet numbers.Evidently, there

Counter-rotating beads
The time-dependent hydrodynamic properties of counterrotating beads can also be studied.Here, the driving forces in eqn (40) are given by Similar mean eld calculations as for the co-rotating beads can be employed and suggest a slow exponential decay of 4 1 + 4 2 / 0. Numerical investigations of eqn (40) without random forces and sound modes qualitatively agree with this prediction, although the dynamics of the rotors are not well described quantitatively by the mean-eld approach.
Interestingly, in our MPC simulations, the rotors exhibit an anti-synchronized stationary state, i.e., 4 1 + 4 2 / p, in contrast to the theoretical results.Whether this is due to the multiplicative noise term in eqn (40), or due to sound modes remains unclear, since a mean-eld approach is not applicable.In any case, the stationary synchronized state is reached at markedly longer time scales, i.e., the synchronization times increase considerably.

Summary and conclusions
We have performed analytical calculations and computer simulations to elucidate the inuence of time-dependent hydrodynamic interactions, thermal uctuations, and uid compressibility on the synchronization of the rotational motion of microrotors.We nd that the presence of the linear unsteady acceleration term in the Navier-Stokes equation leads to synchronization of the rotational motion.Synchronization even prevails in the presence of thermal uctuations, as shown by our mesoscale simulations exploiting the MPC method.Thereby, the simulation results are well described by our mean-eld analytical theory.Fluid compressibility affects the synchronization time over a certain range of Péclet numbers.Thereby, compressibility implies larger synchronization times compared to incompressible uids.This is related to the instantaneous propagation of hydrodynamics in incompressible uids. 91aturally, synchronization by time-dependent hydrodynamic interactions is possible, because the time reversibility of the underlying dynamical equations is broken, despite a timereversible cyclic rotor motion.Our studies show that the synchronization time exhibits the dependence t s =T $ 1= ffiffiffiffiffiffiffiffi Re T p on the oscillatory Reynolds number Re T ¼ s n /T, with s n the viscous time scale and T the period.Hence, t s approaches innity for Re T / 0, i.e., for systems where the period is far larger than the viscous time.However, the transition is gradual, which implies a continuous breakdown of synchronization in analogy to the continuous breakdown of the scallop theorem with decreasing Reynolds number Re. 45,47 Only for Re T ¼ 0, there is no synchronization anymore, again in analogy with the scallop theorem.
Our MPC simulations emphasize the importance of timedependent hydrodynamic interactions as a synchronization mechanism.Aside from physical aspects, this is also important from a simulational point of view.In mesoscale simulations, such as MPC, a clear separation of time scales, e.g., s c ( s n ( T, is oen not possible and Reynolds numbers are rather on the order of 10 À1 to 10 À2 than zero.Thus, time-dependent hydrodynamic interactions are present and relevant, and results have to be interpreted with care when compared with results of the Stokes equation.In this context, we would like to mention that the present nite Reynolds number has very little effect on  synchronization, as is evident by the agreement between simulations and analytical theory. The simulations reveal a strong dependence of synchronization on the size of the simulation box.Thereby, the nite size effect becomes more pronounced for smaller Péclet numbers and the determination of the asymptotic behavior becomes increasingly demanding.The box-size dependence can be understood as follows.No box-size effect appears when the rotors synchronized their motion before the sound traversed the box length L, which requires the time s L c ¼ L/c, hence, t s /s L c < 1.With eqn (45) and the intrinsic frequency u $ F, we obtain t s /s L c $ (R/R H ) 3/2 /(LF 3/2 ).Thus, a reduction of the force has to be compensated by an increase of the box size in order to maintain a small ratio t s /s L c .Moreover, thermal uctuations are signicantly more pronounced for small Pe, which impedes the extraction of the average synchronization behavior.We emphasize, however, that the rotor motion synchronizes even for small Péclet numbers; only the numerical effort increases.
A critical comparison with typical values of microswimmers, e.g., Chlamydomonas, reveals that time-dependent hydrodynamic interactions are present, but are most likely not the dominant mechanisms for synchronization of their agella motion.In contrast, systems of micrometer-size objects, e.g., colloids or biological cells, can be tailored such that hydrodynamic interactions are most relevant.This might open opportunities to exploit hydrodynamic phenomena for novel sensing applications with tailored set-ups.

Fig. 3
Fig. 3 Flow fields of hydrodynamically coupled rotors for an asynchronous (top) and a synchronized (bottom) state for the P éclet number Pe ¼ 120.The flow field of a point particle diverges, hence, we introduce the maximum cut-off velocity v 0 .

Fig. 4
Fig. 4 Time dependence of phase-angle differences for three realizations with the initial value D(0) ¼ p/2 and Pe ¼ 180.

Fig. 5
Fig. 5 Average phase-angle difference D(t) for the P éclet numbers Pe ¼ 120, 140, 160, and 180 (top to bottom).Symbols represent simulation results and the solid lines are obtained by numerical integration of eqn (41) including transversal and longitudinal modes of the hydrodynamic tensor.

Fig. 6
Fig. 6 Distribution of the phase-angle difference in the synchronized state for Pe ¼ 100.The red line indicates a wrapped normal distribution with parameters estimated according to eqn (47).

Fig. 8
Fig. 8 Average phase-angle differences for the simulation box sizes L x /a ¼ 25, 50, 75, 100 and 150 (bottom to top).The box dimension L z ¼ 20a is fixed and the P écelt number is Pe ¼ 100.

Fig. 9
Fig. 9 Phase-angle difference D(t) for Pe ¼ 120.The solid line follows by numerical integration of eqn (41) including transversal and longitudinal modes of the hydrodynamic tensor, whereas the dashed line includes transverse modes only.Symbols are ensemble-averaged simulation results and the vertical line indicates the definition of s D .

Fig. 10
Fig. 10 Maximum deviations s D (eqn (48), bullets) of the phase-angle difference and deviations s s (eqn (49), diamonds) between synchronization times for various P écelt numbers of compressible and incompressible fluids.
The radius of a circle is set to R ¼ 2a and the distance to d ¼ 5a in terms of the MPC length scale.In simulations, forces in the range F/(k B Q/a) ¼ 10-100 are considered, which corresponds to the Péclet numbers Pe ¼ FR/(k B Q) ¼ 20-200.Periodic boundary conditions are applied for the uid, with a rectangular square-cuboid box of side lengths L x ¼ L y ¼ 100a and L z ¼ 20a.
80Qa 2 =m p z 0:0023 and the hydrodynamics radius R H ¼ k B Q/(6phD 0 ) z 0.28a assuming no-slip boundary conditions.80 , for both, the transversal and longitudinal part, i.e., A ¼ A T + A L and B ¼ B T + B L .For the transverse part, the closed expressions 3e Àik$r Qðk; tÞ(34)is known as the dynamical Oseen tensor.In general, the matrix elements of Q(r, t) can be expressed as (r ¼ |r|)