Optimal Hydrodynamic Synchronization of Colloidal Rotors

Synchronization of driven oscillators is a key aspect of flow generation in artificial and biological filaments such as cilia. Previous theoretical and numerical studies have considered the ‘‘rotor’’ model of a cilium in which the filament is coarse grained into a colloidal sphere driven with a given force law along a predefined trajectory to represent the oscillating motion of the cilium. These studies pointed to the importance of two factors in the emergence of synchronization: the modulation of the driving force around the orbit and the deformability of the trajectory. In this work it is shown via experiments, supported by numerical simulations and theory, that both of these factors are important and can be combined to produce strong synchronization (within a few cycles) even in the presence of thermal noise.

Hydrodynamic synchronization can occur in a broad range of microscale systems, interacting through low Reynolds number (Re) flows [1].An open question in biological systems, such as ciliated cells and tissues, is whether synchronization of the filaments is set purely (or at all?) by their hydrodynamic interaction [2,3], or if other mechanical forces, for example, acting through the cell base are important [4,5].In a simple model of two actively driven spheres linked by a scaffold to a third passive sphere, representing the swimmer algae Chlamydomonas, it has been shown that hydrodynamic interaction between the flagella may be negligible during swimming [4]; however, synchronization is also seen between physically immobilized cells [6].These are important issues, because flows driven by motile cilia allow a local circulation to be maintained within various organs and confer motility to a number of organisms [7]; severe human pathologies result from malfunction of cilia activity at various levels [8].The biological systems are so complex that there is a need to work with simplified systems in order to perform controlled experiments and gain a clear understanding of the underlying physics [1,6,[9][10][11][12][13][14].These experimental model systems might themselves be seen as prototypes for future microfluidic technology, for example, to make micromixers and micropumps with carpets of driven filaments [15].In addition to aiming to establish the nature and stability of the coupling, questions of optimization and efficiency can be addressed [16].
With current technology it is possible to build micronscale phase oscillators that exhibit hydrodynamic synchronization and are simple to describe theoretically, allowing quantitative studies of collective dynamics.Two simple models have been studied recently as actively driven units: (a) a configuration-dependent geometric-switch system [17][18][19], which has been realized experimentally with feedback-controlled optical tweezers [13,[20][21][22], and (b) a ''force-controlled'' model in which particles are driven along predefined orbits with defined force profiles, which has been studied numerically and theoretically [23,24], and has been used to describe collective states of many oscillators [25].Here, we consider an experimental realization of this latter, ''rotor'' model.
A general principle is that synchronization at low Re is possible only by breaking time reversal symmetry [22,26,27].In this Letter we present an experiment on the synchronization of two rotors, addressing the open question of the relative importance of the two main factors which can lead to this symmetry breaking: The radial flexibility, which was identified in [23] as an essential component, allowing trajectories to deform, and in contrast to this, the modulation of the driving force, which was proposed in [24] as a sufficient (and possibly generically more relevant) component in enabling rotor synchronization.The experimental findings are compared with fully stochastic Brownian dynamics simulations [28], including hydrodynamic interactions through the Oseen and Rotne-Prager descriptions, and to an analytical approximation.
The motion of two externally driven spherical particles at low Reynolds number is described by the force balance where i 2 f1; 2g indexes the bead.F i represents the driving force acting on bead i, and _r i its velocity.The Oseen tensor H [29] [see the Supplemental Material (SM) [30] for a discussion of near-field and solid boundary effects] describes the drag, and f i is the stochastic Brownian force on bead i.For the time scales considered here, the noise is adequately described by hf i ðtÞi ¼ 0 and hf i ðtÞf j ðt 0 Þi ¼ 2k B TH À1 i;j ðt À t 0 Þ [31].The rotors with circular trajectories of radius r are implemented experimentally by using feedback-controlled 228103-1 Ó 2013 American Physical Society optical tweezers, in which the position of the optical trap is updated based on the position of the particle, and maintained a distance Áx ahead of the particle (Fig. 1).Two very different optical tweezer systems are used, allowing different regimes of time scales to be explored.A timeshared trapping laser based on acousto-optical deflection (AOD) of a beam is used at the University of Cambridge, which results in faster rotors (typical period, with constant driving force, 0.7 s), whereas holographic optical tweezers (HOT) are employed at the University of Bristol (typical period 5 s); the setups are described in detail, respectively, in [12,20,32].In both cases the particle positions are obtained by online image processing; a constant force can be maintained on the colloidal particle if the trap positions are updated on a time scale much shorter than the relaxation time of the particle in the trap ¼ =k, where k is the trap stiffness and ¼ 6a is the drag coefficient on a sphere of radius a in a liquid of viscosity .Experimentally there is a limit implementing this ''feedback''; two time intervals have to be considered: the finite sampling time and a further delay time to update the laser position.It will be seen that these limitations do not affect the qualitative behavior of the coupled systems, but they change the synchronization strength compared to the theoretical expectation.The trap potentials from a single beam are to a good approximation harmonic, and one clearly can obtain the same force by changing Áx or k, keeping their product constant.Interestingly, varying one against the other results in the same orbit period (for a given trajectory) but a different stiffness k r in the radial direction.
In the AOD rotors experiment, a ¼ 1:74 m (silica beads, Bangs Labs), r ¼ 3:2 m, the distance between rotors is d ¼ 25:4 m, and ¼ 6 mPa Á s.In the HOT rotors experiment, a ¼ 1:67 m (polystyrene beads, Corpuscular), r ¼ 2 m, d ¼ 10 m, and ¼ 85 mPa Á s.In both setups, the beads are diluted in a water-glycerol solution to set the desired viscosity.Except when stated otherwise, the experiments described in this Letter are the AOD experiments.Trapped beads are maintained at least 50 m from any surface of the microscope slides in order to minimize effects of wall interaction.
Significant efforts are made to achieve the desired driving force FðÞ and radial stiffness.Before a coupling experiment is carried out, the driving force is calibrated on each rotor individually; see Fig. 1.Force profiles of the form FðÞ ¼ F 0 ½1 À A 2 sinð2Þ, with F 0 % 5:5 pN, were chosen since they correspond to the most efficient profiles to generate synchronization in circular trajectories [26].The force profile is controlled by changing the distance between the bead and the trap depending on its angular position .For each set of parameters, a first profile is measured on a single bead for which is set proportional to 1 À A 2 sinð2Þ.This is then fine-tuned, minimizing the relative difference between the expected and measured profiles.
The trapping force along the radial direction is tuned by creating a line of 21 traps, scanned in a random order in about 2 ms and perpendicular to the tangent direction.This creates a gradient of laser intensity and a potential landscape that is well approximated by a harmonic potential, of spring constant k r , in the range of the fluctuations in the position of the particles along the perpendicular direction.
The relative phase of the two rotors can be measured; this remains constant over successive cycles of oscillation if there is phase locking.In the systems studied here, synchronization is always in phase.Two measurements of the strength of synchronization are possible: the time required to reach synchronization from an arbitrary initial condition or, exploiting the presence of Brownian noise, the relaxation time scale in the autocorrelation of the phase difference at steady state.In both cases, the ''natural'' unit of time is the period of one rotation.When using the arbitrary initial condition method, if A 2 ¼ 0 the phase difference decays according to Eq. ( 28) in [23], which, for small phase differences, describes an exponential decay with the same time scale as the autocorrelation function.
Figure 2 shows the effect of both parameters k r and A 2 .From Refs.[23,24], it is expected that increasing the flexibility (by decreasing k r ) or increasing the asymmetry A 2 of the driving potential should lead to stronger synchronization (lower relaxation time).The experiments confirm these two theoretical predictions.
Although the trends in Fig. 2 are clear, the raw experimental relaxation times, observed in both the AOD and FIG. 1 (color online).(a) Two colloidal ''rotors'' are made by driving particles along predefined closed circular trajectories with optical trapping potentials.The driving force F, pointing towards the trap's minimum (marked by red crosses), may be maintained constant or modulated as a predefined function of the phase angle .The time dependence of is not predefined and evolves under the net action of the driving force and any other forces acting on the particle.(b) The orbital velocity can be either held constant or modulated so that it is anisotropic: the functional form FðÞ ¼ F 0 ½1 À A 2 sinð2Þ has been chosen in this work, studying the effect of the modulation parameter A 2 following [24].The velocity of a single orbiting particle is shown here for A 2 ¼ , 0.3 and 0.85 (increasing modulation).Markers are measured velocities (proportional to force at low Re) and lines are the expected shape, proportional to ½1 À A 2 sinð2Þ.HOT experimental setups, were found to be systematically higher than predicted in the simulations; i.e., the rotors synchronize slower than expected.This discrepancy was investigated carefully, and it appears that several factors act in this direction: The main contribution comes from the delay in moving the traps ahead of the particles.Simulations (shown in the SM) show that a realistic value of the experimental delay (5-20 ms) leads to an increase in the relaxation time by nearly a factor of 2. Smaller contributions are attributed to wall effects, to a detuning between the intrinsic periods of the two oscillators, and to the uncertainty in the measurement of k r (see the SM for a detailed analysis).We have found empirically that these different effects can be accounted for by adjusting the value of k r : the experiments performed for A 2 ¼ 0:0 were compared with simulations, determining that simulations run with a stiffness higher by a factor 2.21 than the experimentally measured k r matched the experimental data.This corrective factor to k r is then fixed, and applied in all the simulation data in Fig. 2, obtaining very satisfactory agreement with the other experimental data sets.
The experiments show clearly that both the radial flexibility and the modulation of FðÞ contribute to the strength of synchronization.A rigorous calculation of the strength of synchronization including both parameters is presented in the SM, adapted from [26].Here, we simply account for both the Lenz (À L ) and Golestanian (À G ) dimensionless control parameters (known respectively from [23,24]) by linear superposition.Then the strength of synchronization À is the sum of the components À L and À G and describes the evolution of the phase difference by a discrete equation of the form Áði þ 1Þ À ÁðiÞ ¼ ÀÀÁðiÞ; (2) where Á ¼ 2 À 1 , i indexes the cycle (counted on either rotor), and In Eq. ( 3), the term 3a=ð4dÞ sets the scale of the hydrodynamic coupling.The two terms in brackets correspond respectively to the contribution of the Lenz model (L) (involving radial stiffness), in which r is the rotor radius, and the Golestanian model (G) (based on the driving force profile) to the synchronization.Note that the Lenz contribution itself also depends on A 2 , since the period of the oscillators depends on the force profile.From Eq. ( 2), the relaxation time is then related to À L and À G by This theoretical model combining both the L and G models is plotted in Fig. 2 (solid lines) to match experiments, and in Fig. 3 in a wider dimensionless parameters range.The model agrees very well with simulations (with the same k r ) and experiments (with k r corrected as described above).Depending on the choice of parameters F 0 , k r , r, and A 0 , the synchronization can be dominated by either the Lenz or Golestanian contribution.A threshold between these two regions can be defined when the two terms in brackets in Eq. (3) are the same; this boundary line is plotted in Fig. 3.We now investigate the effect of thermal fluctuations.Brownian motion leads to fluctuations of Á over time.In the limit of small fluctuations, an analogy with a particle in a confining harmonic potential can be made.The fluctuations in displacement x of such a particle satisfy hx 2 i ¼ k B T=k, with k the harmonic trapping constant.We can then match the relaxation time of the confined particle with the relaxation time of the phase difference: =k $ t 0 =À, where FIG. 2 (color online).Relaxation time depending on parameters k r and A 2 .Experimental results (markers), simulations (dotted line), and theoretical value (solid line) from Eq. ( 4), well defined only for weak coupling.In both simulations and the analytical formula, a correction has been applied to k r , as described in the main text.In (a) the experimental k r is varied for A 2 ¼ 0 (plusses), 0.4 (circle), and 0.7 (square).In (b) A 2 is varied for k r ¼ 1:7 (plusses), 4.0 (circle), and 9.9 (square) pN=m.Decreasing k r and increasing A 2 produces stronger synchronization.The inset shows the transient to phase locking between 1 and 2 , for different initial conditions, in the HOT setup (shown here with A 2 ¼ 0, k r $ 10 pN=m, and F 0 $ 8 pN).It typically takes a few orbits for the phase difference to decay into the stable in-phase synchronized state.When both k r is small and A 2 is close to 1, the strength of synchronization À is not small compared to 1: theory becomes inaccurate and the argument of the log in Eq. ( 4) can even become negative, as in some sections for k r ¼ 1:7 pN=m in (b).t 0 is the intrinsic period of the rotors.We also have hx 2 i $ hðrÁÞ 2 i.Rearranging, this leads to The simulation data in Fig. 4(a) show hÁ 2 i plotted against the quantity above, for a large variety of parameters and at low temperature (ensuring small fluctuations).The data collapse onto lines, with a slope that depends only on A 2 .Since synchronization will be lost above a certain threshold of hÁ 2 i, the scaling above gives insight into the range of parameters that will allow synchronization in the presence of noise.The slopes from 4(a) are plotted in the inset graph and show a nontrivial dependence on A 2 , initially constant and then growing at larger A 2 .We do not have a full explanation of this, but we note that with nonzero A 2 the synchronization strength varies as the beads go around their orbits; it would then be expected that the variance of fluctuations is affected in a way that is not simply described by the cycle average of the potential, as discussed in [33].The simulation data in Fig. 4(b) show the applicability of Eq. ( 5) at higher temperatures, where fluctuations are large.As the temperature increases, the variance becomes superlinear, implying that the restoring strength is less than linear for large phase shifts.At very high temperatures, synchronization is lost and the distribution of Á becomes uniform in the finite range [0, 2], so that the variance converges to a constant value.
To investigate experimentally the effect of noise, we varied the distance between the rotors, as shown in the inset of Fig. 4(b).This is equivalent to varying the temperature, provided that the fluctuations are small enough such that Eq. ( 5) holds.There is good agreement between experiment and simulation (as before, with corrected k r ) for rotor separations & 30 m (first three points), and an increasing discrepancy at larger separations.Possible reasons for this discrepancy include wall interactions, which become more important at large d, or the other experimental problems discussed in the SM, emerging more significantly at low coupling.Experimentally, synchronization is lost between d ¼ 44 and 48 m.
In the future, it should be possible to parametrize and coarse grain complex biological cilia, or flexible artificial systems driven by complex mechanisms, into a  3) and ( 4) (mesh surface).The thick line separates the regions in which synchronization is dominated by the flexibility of the rotor (Lenz synchronization, LS) or by the angular modulation of the driving force (Golestanian synchronization, GS).For À L þ À G close to 0, synchronization can be lost because of thermal fluctuations; simulations were performed in a range in which rotors are synchronized, for which the relaxation time can be easily extracted.FIG. 4 (color online).Thermal noise leads to fluctuations of Á around 0 (the in-phase state).(a) Simulations for variable values of A 2 , k r , a, r, F 0 , d, and in a regime of small fluctuations (T ¼ 10 K).The variance of Á exhibits a linear dependence on t 0 k B T=ðÀr 2 Þ, with a coefficient that depends only on A 2 (as shown in the inset).(b) For a ¼ 1:74 m, d ¼ 35 m, r ¼ 3:2 m, F 0 ¼ 5:5 pN, k r ¼ 12:97 pN=m, and A 2 ¼ 0:0, the simulation temperature is varied from 0 to 1:5 Â 10 5 K.The line shows the variance from Eq. ( 5) multiplied by the coefficient obtained from (a).The inset is an enlarged portion of the main graph (plusses), with experiment results (square, experimental k r ¼ 5:9 pN=m) obtained by varying the rotor separation between 20 and 44 m.Varying d is equivalent to varying T when the fluctuations are well described by Eq. (5).Synchronization was lost between d ¼ 44 m (last point on the graph) and d ¼ 48 m (not represented: phase slips were observed and hÁ 2 i ¼ 1:7 rad 2 ).

1 )FIG. 3 (
FIG. 3 (color online).Relaxation time depending on the Lenz (À L ) and Golestanian (À G ) dimensionless control parameters.Simulations (circles) are in excellent agreement with the theoretical model in Eqs.(3) and (4) (mesh surface).The thick line separates the regions in which synchronization is dominated by the flexibility of the rotor (Lenz synchronization, LS) or by the angular modulation of the driving force (Golestanian synchronization, GS).For À L þ À G close to 0, synchronization can be lost because of thermal fluctuations; simulations were performed in a range in which rotors are synchronized, for which the relaxation time can be easily extracted.