Dynamic mode locking in a driven colloidal system: experiments and theory

In this article we examine the dynamics of a colloidal particle driven by a modulated force over a sinusoidal optical potential energy landscape. Coupling between the competing frequencies of the modulated drive and that of particle motion over the periodic landscape leads to synchronisation of particle motion into discrete modes. This synchronisation manifests as steps in the average particle velocity, with mode locked steps covering a range of average driving velocities. The amplitude and frequency dependence of the steps are considered, and compared to results from analytic theory, Langevin Dynamics simulations, and Dynamic Density Functional Theory. Furthermore, the critical driving velocity is studied, and simulation used to extend the range of conditions accessible in experiments alone. Finally, state diagrams from experiment, simulation, and theory are used to show the extent of the dynamically locked modes in two dimensions, as a function of both the amplitude and frequency of the modulated drive.


I. Introduction
Synchronisation is one of the most diverse fundamental physical phenomena [1]. From Huygens' pendulum clocks 350 years ago [2,3] to fireflies [4], applause [5,6], and animals' circadian rhythms [7], frequency entrainment occurs all over the natural and technological world. The phenomenon occurs when weakly coupled competing oscillators adjust their rhythms to match each other [4]. Synchronisation on the micro-scale is of technological importance, as the decreasing size of electronic and mechanical systems demands ever-smaller frequency references [8]. Recent developments include electromechanical [9,10] and optomechanical oscillators [11,12], but such systems are limited in their scalability [12].
Dynamic mode locking is a synchronisation phenomenon that occurs when systems with a natural internal frequency are driven by an external modulation. Competition between the two frequencies leads to coupling, causing the system to synchronise into repeating modes of motion. Previous work has sought to understand dynamic mode locking in superconductor vortex lattices [13][14][15][16][17], but the difficulty in visualising such systems makes model systems necessary [18]. Other systems showing such resonance behaviour include driven adatoms on atomic surfaces [19,20], and Josephson junctions [21][22][23]. The AC Josephson effect occurs when the tunneling electron pairs at an insulated superconductor junction are driven with an AC and DC current [24]. Regions appear where resistance does not increase with increasing DC current [23,24], and the shape of the resulting graph is known as 'Shapiro Steps'. Charge density waves are another technologically significant system demonstrating dynamic mode locking, and have been extensively studied experimentally and numerically [25][26][27].
Model systems composed of colloidal particles in periodic potentials have been studied for a number of years [28,29], from simple double well potentials [30][31][32][33] to directed motion [34][35][36][37][38][39], particle sorting [40], and kink generation [41] in two-dimensional (2D) optical lattices. Colloidal systems are easy to manipulate, and have accessible length and time scales, making them attractive models for the study of synchronisation at the micro-scale. Noise in Brownian systems has been found in theory to induce anomalous diffusion [42] and stochastic resonance [43][44][45], and rocking-ratchet like potentials have been used in optical and magnetic systems [46,47]. The possibility of resonance has also been explored in systems with feedback [48] and random pinning potentials [49]. Recent work studied the transport properties of a system of magnetically driven colloidal particles [50]. Recent theoretical work also examines the possibility of producing mode locking steps in 2D colloidal monolayers [51,52].
Here, a system of Brownian particles is driven over a sinusoidal optical potential energy landscape by a driving force consisting of constant and modulated parts. The natural frequency of the particle driven over the optical potential energy landscape by the DC component of the driving force couples to the frequency of the AC component. As we have shown previously [53], this coupling leads to dynamic mode locking. This work considers the frequency and amplitude dependence of the synchronisation, through experiments, Langevin Dynamics simulations and Dynamic Density Functional Theory (DDFT). The three complementary approaches are used together to build a comprehensive picture of dynamic mode locking. Firstly, the theoretical and simulation approaches are introduced in Section II, including an analytical approximation. The experimental methods are described in Section III. Results from all of the approaches are described and discussed in Section IV, including mode locked steps, state diagrams, and critical driving forces.

II. Theory and computer simulations
A. Langevin dynamics To describe a Brownian particle driven by the sum of a constant and a modulated force across a periodic optical potential energy landscape, the overdamped Langevin equation is written as: where the particle velocity, v, at position x and time t depends on the force from the optical potential energy landscape, F T , the Brownian force, ξ(t) (modelled as Gaussian white noise with a mean of zero and variance of 2ζk B T , where k B T is thermal energy), the constant driving force, F DC , the friction coefficient, ζ, and the oscillating driving force, Here, F AC is the amplitude of the modulated driving force, and ω = 2πν is the angular frequency, where ν is the frequency of the applied oscillation. Note that in this paper, 'DC' and 'AC' are used only in analogy to direct-and alternating-current, and refer to constant-and oscillating-velocity drives respectively.
The optical potential energy landscape, U T (x), is taken to be sinusoidal, as described in references [53][54][55]: where k is the trap stiffness, V 0 is the trap strength, and λ is the wavelength of the landscape. Equation 3 leads to an optical force [55]: where the critical driving force, F C , is given by the following equation [55]: Thus, the full equation of motion for a particle driven by DC and AC driving forces over a sinusoidal optical potential energy landscape is given by: Note that the critical driving force, F C , is a property of the landscape, and is equal to the DC critical driving force found in reference [55]. As the total driving force in equation 6 is a sum of the DC and time dependent AC contributions, the critical DC driving force required to overcome pinning depends on the amplitude and frequency of the modulated component of the driving force.

B. The 'high frequency' theory
While the equation of motion in equation 6 is not analytically soluble, useful insight can be obtained in the limit of high driving frequency (ν ≫ F C λζ) in the absence of noise. Within this approximation, it is possible to obtain an effective Adler equation [56,57] similar to that found for the case of constant drive alone [55]. Thus an expression for the average velocity may be written (see Appendix A for full details):  where ∆F DC = F DC − rλνζ, r = 0, ±1, ±2, . . ., and J m is the mth order Bessel function of the first kind. This 'square root law' expression is the AC driven counterpart to the simpler form found for the case of DC drive alone [55]: The condition for the approximation, ν ≫ F C λζ (see equation 19, Appendix A), means that on a landscape with a trap spacing of λ = 3.5 µm (F C ζ ≈ 1.8 µm s −1 ), the high frequency regime is valid when ν ≫ 0.5 Hz.
The dependence of the average particle velocity, v, on the driving velocity, F DC ζ, according to equation 7 is shown schematically in figure 1(a). Equation 7 describes mean particle velocity 'above' and 'below' critical points, with two critical points found for every absolute value of r, in contrast to the DC only case, which has only a single F C . Between each pair of critical points, a 'subcritical' regime exists, where the particle velocity is constant, corresponding to mode locked steps. The form of this dependence is analogous to the Shapiro Steps seen in Josephson Junctions [23], and also in charge density wave systems [25] and vortex lattices [15][16][17]58].
The two critical points F CRIT at the ends of resonant step r are found by determining F DC at the condition where the two different solutions in equation 7 coincide, ∆F DC = ± F C J −r (F AC (λνζ)) , at which point the square root vanishes. As a result, by recalling the definition of ∆F DC and by replacing F DC with F CRIT , the following is obtained: The amplitude (F AC ) and frequency (ν) dependence of F CRIT defines state diagrams, with regions containing locked modes enclosed by pairs of critical lines. Figure 1(b) shows such a state diagram as a function of F AC , for a particle driven with a frequency of ν = 3 4 Hz across an optical landscape with a trap spacing of λ = 3.5 µm. Each colour and value of r represents a single mode locked velocity. The state diagram is formed from twisted 'Arnold Tongues' [1], where each separated region of the same colour actually represents a different dynamic mode with the same average velocity [53]. The second critical point of the 'zeroth' step appears as the effective critical driving velocity, F C,eff , below which the particle is pinned to the landscape and does not slide.

C. Dynamic Density Functional Theory
The Langevin picture is stochastically equivalent to the Smoluchowski picture, in which the temporal evolution of the probability density distribution, p(x, t), of the particle position is studied rather than the stochastic trajectories of individual particles. The Smoluchowski equation can be seen as a special case of the Dynamical Density Functional Theory (DDFT) in the absence of interparticle interactions [59][60][61]. The governing equation for the probability density distribution is given by where D is the diffusion coefficient and F (x, t) = F DC + F AC cos(2πνt) − F C sin(2πx λ) is the total force acting on the particle. Equation 9 is solved numerically using a finite volume partial differential equation solver [62]. As an initial condition p(x, t = 0), a very narrow Gaussian distribution is chosen. See Appendix B for more details.
Within the Smoluchowski picture, averages of statistical quantities are defined by weighting these quantities with the particle probability distribution p(x, t), i.e. ⟨a⟩(t) = ∫ ∞ −∞ dx a(x)p(x, t). These averages are stochastically equivalent to noise averages performed in the Langevin picture. Thus, the mean particle position is ⟨x⟩(t). The mean velocity is further defined as the change in the mean particle position in time: where overbar denotes a time average. As a measure of the fluctuations around the mean particle trajectory, the variance of the particle probability distribution is considered: In the context of this work, if the standard deviation, σ(t), is much smaller than the trap spacing then almost all possible particle trajectories end up in the same trap as the mean particle position after time t. If the standard deviation is larger than λ then possible particle trajectories end up distributed in potential wells surrounding the mean. Particle fluctuations around the mean position may be quantified using an effective long-time diffusion coefficient, defined from the variance:

A. Colloidal model system
The colloidal system is composed of Dynabeads M-270 carboxylic acid (diameter 3 µm), in 20% EtOH aq , held in a quartz glass sample cell (Hellma) with internal dimensions of 9 × 20 × 0.2 mm. Particles are much more dense than the solvent, and sediment into a single layer near the bottom of the sample cell. The coefficient of friction, ζ, is found from diffusion to be 9.19 × 10 −8 kg s −1 , slightly higher than expected from Stokes friction (ζ Stokes = 6πηa with η the viscosity), due to the proximity of the particles to the wall. Particle concentration is low, so that only a single particle is visible in the field of view.

B. Experimental setup and parameters
The experimental setup consists of an infra-red (1064 nm) laser, controlled using a pair of perpendicular acoustooptical deflectors, and focused using a 50×, NA=0.55 microscope objective [54]. The one-dimensional periodic optical landscape, with trap spacing λ = 3.5 µm, is generated in Aresis Tweez software controlled from a LabView interface. A landscape with this trap spacing may be treated as sinusoidal, as shown in reference [55]. The traps are time-shared at 5 kHz, such that on the time scale of the particles (with a Brownian time of ∼ 50 s, and at least ∼ 1 3 s to be driven one trap spacing at a given F DC ), the traps form a constant potential energy landscape. The laser power and the total number of traps are held constant throughout the experiments, so that the laser power per trap is consistent. A laser power of 350 mW is set and 46 traps are used, corresponding to ∼ 0.75 mW per trap at the sample. This gives typical values of trap stiffness, k = 3.8 × 10 −7 kg s −2 , and trap strength, The driving force is provided by a PI-542.2CD piezo-stage, controlled using the LabView interface, at driving velocities of 0.05 ≤ F DC ζ ≤ 8 µm s −1 . AC driving velocity is added to the DC drive, with an amplitude 0.4 ≤ F AC ζ ≤ 14 µm s −1 , and a frequency of 1 10 Hz ≤ ν ≤ 2 Hz. Images are focused onto a Ximea CMOS camera using a 40×, NA= 0.50 microscope objective, and the particle position is recorded live at 40 Hz from the camera image.

C. Average velocity experiments
To obtain plots of average particle velocity against driving velocity, six repeats across the whole potential landscape are made at each driving velocity for each amplitude and frequency of the AC drive. Average velocity, v, is found by linearly fitting the particle trajectory, x(t), over an integer number of periods of the oscillation.

D. Critical driving velocity experiments
The critical DC driving velocity is defined as the DC driving velocity at which the particle starts to slide irreversibly across the optical potential energy landscape. It is found by iterating the DC driving velocity, with a maximum resolution of 0.05 µm s −1 . A particle is said to be pinned if it still returns to its starting lattice position after the stage has moved 100 µm, or three minutes has elapsed, whichever happens first. The region in which the particle does not irreversibly slide is essentially the zeroth mode locked step, and the effective critical driving velocity is therefore the second critical point of this step (see section II B).
The critical driving velocity found here is not the critical driving velocity found in our previous work ( [55]), F C ζ, because the total driving force in equation 6 is a sum of the DC and time dependent AC contributions. Therefore the critical driving velocity measured here is an effective critical driving velocity, F C,eff ζ, as it is not purely a property of the landscape. However, clearly when F AC = 0, F C,eff ≡ F C .

IV. Results and discussion
Results are presented which show the amplitude and frequency dependence of the mode locked steps and state diagrams illustrated in figure 1. Experimental results and Dynamic Density Functional Theory (DDFT) computations are compared to Langevin Dynamics (LD) simulations and the analytic approximation for the high frequency limit (Section II B) as appropriate. In general, good quantitative agreement is found between the various approaches.

A. The mode locking steps
The effect of introducing the oscillating force term to the equation of motion (equation 6) on v as a function of F DC ζ is shown in figure 2(a). Here, data with a modulated force of amplitude F AC ζ = 5.2 µm s −1 and frequency ν = 3 4 Hz (△) is compared to the F AC = 0 (◯, i.e. DC drive only) case for a landscape of trap spacing λ = 3.5 µm (see [53,55]). The case of F AC ≠ 0 follows the 'Shapiro Steps' form illustrated in figure 1(a). The effective critical driving velocity for the modulated case is almost zero, after which v increases until it is significantly larger than that expected for a free particle (grey line), implying that the particle is moving on average more quickly than the piezo stage. The average velocity then plateaus on the first resonant step, at v = 1λν = 2.625 µm s −1 . The step extends over a range of driving velocities, and then v increases after the second critical point, to meet another step, at twice the average particle velocity of the first.
The solid and dashed lines on figure 2(a) show results from Langevin Dynamics simulations, both with and without the noise term. The steps found from the experiments are faithfully reproduced by the simulations, with the inclusion of noise obviously important in this system of Brownian particles. The effect of noise is important in the vicinity of the critical points, where it is seen to round the edges of the steps. Figure 2(b) compares results from DDFT (upper panel), and the high frequency approximation (lower panel, see section II B) to the experimental data. The high frequency theory results, from equation 7, are calculated with respect to each critical point, and it is notable that although the step positions are largely captured, the results between steps, from adjacent critical points, do not necessarily agree.

Variance and diffusion
Next, we perform DDFT calculations and consider fluctuations around the mean particle position, which strongly depend on whether or not the system is mode locked. In figure 3(a) the variance is shown as a function of time for the conditions considered above (F AC ζ = 5.2 µm s −1 , ν = 3 4 Hz). The displayed numerical data correspond to states on the mode locked steps around the first step (see inset). For unlocked states (1.0 ≤ F DC ζ ≤ 2.0 µm s −1 ) the variance grows rapidly, corresponding to an effective diffusion much larger than the free diffusion (dashed line). In the mode locked states (2.5 ≤ F DC ζ ≤ 3.25 µm s −1 ) the variance reaches a long-lived plateau where the diffusion is (nearly) zero, before eventually crossing over. Similar intermediate plateaus have also been observed in underdamped systems [63][64][65] and static systems (F AC = 0) [66,67].
The effective long-time diffusion coefficient (equation 12) is the limit of σ 2 (t) 2t as t → ∞. Plotting D eff as a function of average driving velocity offers additional insight into the mode locked steps. Figure 3(b) shows that effective diffusion is close to zero at the mode locked steps, and much higher between. The typical double-peak signature for D eff , as described in [44], is recovered. The low D eff values in the locked regions result from a vanishing influence of thermal noise which can also be found in related systems [68,69]. This is a symptom of the predictability of the locked state: when the particle is locked into a particular mode of motion, its position on the periodic landscape after a certain time depends purely on the driving conditions. D eff is highest in unlocked states on the cusp of synchronisation conditions, as a small perturbation may cause the particle to jump to the next potential well, or stay in the present one. This corresponds to the discontinuities at the critical points in the schematic in figure 1(a). Figure 3(c) shows D eff for the first step on a log scale, showing that it decreases by ∼ 8 orders of magnitude between the unlocked and locked states. To put this into context, the lowest effective long-time diffusion coefficient D eff ≈ 2.1 × 10 −9 µm 2 s −1 corresponds to the particle being one lattice spacing away from the predicted position after approximately 45 years.

B. Dependence on the amplitude
As the synchronisation condition depends only on the trap spacing and the modulation frequency, changing the modulation amplitude alone does not alter the step velocities. Figure 4(a) shows mode locking steps obtained from both experiment and Langevin Dynamics simulations for four different amplitudes, at a frequency of ν = 3 4 Hz. This shows that there is, however, a strong dependence of the width of the locked step on the oscillation amplitude. For very low amplitude there is a small visible first step, giving a deviation from the zero oscillation data, but no second step is observed; the points lie on top of the F AC = 0 line. As amplitude increases, the first step increases in width, and a second step appears and widens. The first step then appears to narrow. There is generally a good agreement between the experimental data and the LD simulations, with small deviations possibly due to experimental uncertainties such as the variability of the laser power during the experiment. Figure 4(b) shows data at a lower frequency of ν = 1 4 Hz, for four amplitudes. The lower frequency means that a larger number of steps appear in the same range of particle velocities. Four steps are visible in the range shown (which is smaller than that in figure 4(a)), with step width varying widely. Notably, the F AC ζ = 5.2 µm s −1 line has three steps at n = 1, 3 and 4, but no step is visible at n = 2.
It is pertinent at this juncture to compare the results for ν = 3 4 Hz to the high frequency theory (equation 7). Figure  4(c) shows each of the four sets of conditions in panel (a), with the results for each critical point from equation 7 (solid lines) compared to the LD results (dashed lines) and experimental results (symbols) from panel (a). The first observation is that the high frequency approximation appears to better match the data at higher amplitudes. The most likely reason for this is that the Bessel function in equation 7 gets smaller as the argument, which is proportional to F AC , increases (see equation 19 in Appendix A). This means that at a given frequency the condition setting the validity of the high frequency approximation is fulfilled better for higher F AC . Note also that the theory consistently overestimates the step width, as it is deterministic, whereas the critical points in experiments and simulations are  somewhat rounded by noise. Finally, it may be seen that as in figure 2(b), the lines between steps determined from different critical points do not overlap, as the 'square-root law' sections are only valid in the close vicinity of their critical points.

State diagram: Low frequency regime
As was shown in figure 1(b), state diagrams may be constructed which show the extent of dynamic mode locking as a function of F AC ζ and F DC ζ. In our previous work [53], we used such a plot to locate numerous dynamic modes for a driving frequency of 1 4 Hz. In figure 5(a) we compare these prior experimental results (circles coloured according to the integer step number, n) with locked regions from DDFT (background colour) and critical lines from LD simulations in the absence of noise. Note that points corresponding to unlocked states are not shown, for ease of interpretation. There is very good agreement between the DDFT and the LD results, with the only difference being that the regions calculated from DDFT are smaller, due to the presence of noise. Both are in good agreement with the experimental data, except that locked states appear to be found at a slightly higher range of driving velocities in the simulations. Figure 5(b) again shows critical lines from LD simulations, plotted over the effective diffusion coefficient, D eff , obtained from DDFT. Locked states appear as dark blue regions, with fairly broad boundaries where step edges are smoothed by noise. Between the locked states the effective diffusion is higher, as was seen in figure 3(b). D eff is particularly high between closely spaced locked modes, as a small perturbation may lead to the particle becoming temporarily trapped in one mode or the other. The related enhanced increase in variance manifests in the experimental data as wider error bars between the modes in figure 4(b). That the diffusion is so high in these regions probably contributes to the mismatch between the experimental and simulation data in panel (a) -a small change in the experimental conditions could cause the particle to cross mode boundaries.
It is nice to observe in figure 5 that the mode boundary lines oscillate, with the first and second critical lines crossing and swapping identity between the modes. The upshot of this oscillation is that there are conditions in which certain locked modes do not appear, for example it is clear that F AC ζ = 5.2 µm s −1 lies between two regions with n = 2, which corresponds exactly to the missing step observed in figure 4(b). This effect is of course mirrored in the critical driving velocity line, being the upper mode boundary of the zeroth step, resulting in some conditions where the critical driving velocity is zero, for example F AC ζ = 5.2 µm s −1 again.

State diagram: High frequency regime
In section II B, an analytical expression (equation 8) was obtained which could predict the locations of the first and second critical points for each mode locked step. It was found that this expression should be valid in the region where ν ≫ 0.5 Hz. It is not possible to probe this region in detail in the experiments, as at higher frequencies, increased particle velocities are required to obtain higher modes, with a resulting loss in resolution. However, using DDFT it is possible to examine this regime, and obtain a state diagram similar to that from experiment. Figure 6(a) shows a state diagram calculated via DDFT for ν = 3 4 Hz, with the mode locked regions being represented by colours, as in figure 5(a). Also plotted on figure 6(a) are lines calculated from equation 8. There is a remarkably good agreement between the DDFT results and the analytical prediction, showing that the approximation is valid to surprisingly low frequencies. The main deviation occurs at lower amplitudes, as was seen in figure 4(c). As as in figure 5(a), the locked regions from DDFT are slightly smaller than the space between the critical lines, due to the noise term which must necessarily be ommited from the theory.
The mode locking footprint can also be seen in the effective diffusion coefficient: figure 6(b) shows the same calculated lines as panel (a), overlayed on D eff , represented by a colour scale. The locked states are clearly visible as the dark blue regions on the state diagram, where the effective diffusion coefficient drops dramatically as seen in figure 3(c). The unlocked states range from blue, where the particle position is largely predictable, to the yellow and red regions between the mode locked steps. An interesting observation may be made in the region between the zeroth and first modes, where the theory predicts no gap between the second critical line of the zeroth mode and the first critical line of the first mode. The effective diffusion in this region is especially high, indicating that in the stochastic system the particle trajectory is highly unpredictable. It is probable that in this region, the particle rapidly jumps between periods of being pinned to the landscape, and being in the first mode. Dashed lines indicate the position of the first step for each frequency, and some higher n steps for completeness. Dark grey circles and line indicate case for no oscillation [55], light grey line indicates case for no traps.

Critical driving velocity
The state diagrams in figures 5 and 6 show that the critical driving velocity, F C,eff , oscillates as a function of modulation amplitude. In figure 7 the critical driving velocities are shown in isolation, and experimental results are compared to LD simulations and the high frequency theory. Both sets of data, for ν = 1 4 Hz and ν = 3 4 Hz show a Bessel-function-like form, with the range of the oscillations determined by the frequency. Each peak actually represents a different pinned mode [53], and regions occur between the modes where the critical force is close to zero and thermal motion is sufficient to overcome the barriers. The theoretical prediction for the high frequency regime (equation 8) predicts the shape of the experimental data reasonably well, but the LD simulations provide a somewhat better quantitative fit.

C. Dependence on the frequency
The frequency dependence of the step velocities is expressed in the synchronisation condition, v = nλν. Figure 8(a) shows v as a function of F DC ζ at three different frequencies, for an amplitude of F AC ζ = 2.0 µm s −1 . At a lower frequency more steps are seen over the same range of particle velocities as more harmonics are attainable. Indeed at the lowest frequency, ν = 1 4 Hz, the mean particle velocity shows three steps corresponding to the first three integer multiples of νλ. The second and third steps at this frequency therefore coincide with the first steps for the other two frequencies, of ν = 1 2 Hz and ν = 3 4 Hz. As was seen in figure 4(a) and (b), there is generally a good agreement between the LD simulations and the experiments. Figure 8(b) shows the frequency dependence of the average particle velocity for an amplitude of F AC ζ = 5.2 µm s −1 , for frequencies ranging from ν = 1 4 Hz to ν = 3 2 Hz. It is worth noting that, as was seen in figure 4(b), not all possible steps appear; for example there is no visible step at v = 1 2 λ s −1 for ν = 1 4 Hz. Both figures 8(a) and 8(b) show that step width is frequency dependent, in addition to the amplitude dependence shown above. Furthermore, the effective critical driving velocity clearly depends on the frequency: in both cases it is seen to increase with ν, approaching the oscillation-free critical driving velocity of F C ζ ≈ 1.8 µm s −1 .

State diagrams and critical driving velocity
DDFT computation is used to produce frequency-dependence state diagrams, which are not feasable experimentally as very low frequencies require extremely long run times. Figure 9(a) shows four such state diagrams, at a range of amplitudes from F AC < F C to F AC > F C . Also included on two of the state diagrams are effective critical driving velocity lines determined from LD simulations. As with the experiments, the run times required to produce further LD data, particularly at low frequencies, are prohibitively long. From the data presented, however, the origin of the increasing effective critical driving velocity seen in figures 8(a) and 8(b) is clear, although the LD line on the F AC ζ = 5.2 µm s −1 plot in particular highlights that the picture is more complex. At F AC > F C (where F C is the DC critical driving velocity), a series of bumps appear at low frequency, which are too fine to be resolved by the DDFT data. Below F C , however (i.e. at F AC ζ = 1.7 µm s −1 ), these bumps disappear, and the depinning transition is defined as a single monotonic increase. This effect has been noted previously for the Frenkel Kontorova model, for a chain of interacting particles [39,45,70]. Figure 9(b) shows the frequency dependent effective critical driving velocity for F AC > F C in more detail. Experimentally determined values for F C are shown, along with the LD simulation result. The experimental data shows a single smaller bump similar to that seen in the DDFT data, as it is also unable to resolve the numerous smaller bumps shown by the LD data. The large bump at higher frequency, which plateaus at the DC-only critical driving velocity F C , represents the truly pinned state, where the particle does not move at all during a cycle of the oscillation, whereas the smaller bumps and lower frequency represent states where the net particle motion is zero, as it moves but returns to the same potential well at the end of every cycle. Also presented in figure 9(b) is the prediction of the high frequency approximation, which captures the form of the experimental data, but has a higher magnitude due to the absence of a noise term in equation 8, and the LD line for F AC ζ = 1.7 µm s −1 , for comparison.

V. Conclusions
Colloidal particles driven by the sum of constant and oscillating forces through a quasi-one-dimensional periodic optical potential energy landscape have been shown to exhibit rich non-linear dynamical behaviour. Experiments showed that when an oscillating drive is applied, the average particle velocity has a staircase-like dependence on the average driving velocity, where the steps represent states of synchronisation between the particle motion and the substrate potential. These results could be faithfuly reproduced using Langevin Dynamics simulations and Dynamic Density Functional Theory (DDFT), and in conditions of a high driving frequency the data mapped surprisingly well to an analytic approximation. Probing the variance in the particle position and the trajectory diffusion using DDFT showed that the effective diffusion coefficient drops dramatically at the resonant mode-locked steps, explaining why the analytic theory (calculated at 'zero temperature' with no fluctuations) is so successful.
The use of simulation and computation in addition to experimental results allowed a full exploration of the amplitude and frequency dependence of dynamic mode locking. State diagrams showing both the amplitude and frequency dependence of the extent of the locked modes exposed the oscillating nature of the critical lines which define the mode locked steps. These critical lines enclose regions which have been previously shown to represent different dynamic modes with the same net particle motion. Finally, the effective critical driving velocity below which a particle is pinned to the potential landscape has been studied, and it has been shown to have an oscillating dependence on the modulation amplitude, with some conditions having no effective critical driving velocity. The frequency dependence has been shown to be more complex, depending on whether the amplitude of the oscillation is above or below the critical driving velocity defined by the landscape.
By using a combination of experiments, computation, and analytic theory, it has been possible to explore the effect of a very wide range of conditions on dynamic mode locking, thereby giving a solid experimental and theoretical foundation to this dynamic synchronisation phenomenon.
In order to solve equation 6 for a certain range of driving frequencies, the work of Cotteverte et al. [71], Chow et al. [72], and Reichhardt et al. [58] is followed. Ideas developed in these previous works are used to draw an analytical approximation in this work. The first step in solving equation 6 is to neglect noise and split the particle trajectory, x(t), into a part due to the terms independent of x(t) and a deviation from it, caused by terms dependent on x(t): Accordingly, the equation for x 0 (t) is taken to contain the DC and AC parts of the driving force, which can be integrated to yield For the second part, δ(t), we obtain from equation 6: where we have accounted for equations 14 and 15. Using the identity: sin [A sin (ωt) + B] ≡ ∑ ∞ m=−∞ J m (A) sin(B+mωt), where J m is the mth order Bessel function of the first kind, A = 2πF AC (ζωλ) and B = (2π λ)(F DC t ζ + δ), equation 16 becomes: Equation 17 is difficult to solve, so an approximation is made that only the leading term of the sum is retained, where mω = −rω ≈ − 2πF DC λζ (r = 0, ±1, ±2, . . .), determines the mode locking. This approximation is valid at high enough frequencies. Indeed, the variation with time of the leading term with m = −r given by equation 18 is the slowest relative to ωt, 2ωt, . . . of the next to leading terms with m = r ± 1, r ± 2, . . .. Provided that the dependence on δ(t) can be neglected in all next to leading terms, their averages over a period of the external modulation vanish. The scale of δ(t) can be estimated by noticing from equation 17 that dδ(t) dt ∼ F C ζ (cf. [72]) or, more accurately, dδ(t) dt ∼ J −r (A)F C ζ and hence δ(t) ∼ J −r (A)F C t ζ.
By considering the term with the next to slowest variation, m = −r ± 1, we require that ωt ≫ δ(t) to arrive at the condition for the validity of our high-frequency approximation: Equation 17 therefore becomes: where ∆F DC = F DC − rλνζ is a small change in the constant part of the driving force, F DC . Introducing a variable q(t) = (∆F DC ζ)t + δ(t) further reduces equation 20 to the form of an Adler equation [56,57], equivalent to that found for the case of constant drive alone [55]: The expression for the average velocity may then be written, by noting that v = ⟨ dx(t) dt ⟩ = F DC ζ +⟨ dδ(t) dt ⟩ = rλν +⟨ dq(t) dt ⟩, as the oscillating force term (the modulated part of the driving velocity) becomes zero after time averaging: for r = 0, ±1, ±2, . . ..

B. Dynamic Density Functional Theory -further details
Implementation details: The finite volume partial differential equation solver FiPy 3.1 [62] is being used to perform the integration of the Smoluchowski equation. The grid of the computer system consists of 10,000 cells and has a total length of 50 µm (≈ 14 potential wells) with periodic boundary conditions. The computations were terminated whenever the probability distribution was so widely spread that effects of periodicity could not be neglected.
Initial conditions: The mean particle trajectory enters in general a short transient state before sychronising with the external AC driving force. This synchronised state is characterised by a periodic phase that modulates the linear drift of the mean trajectory. In order to suppress effects of the transient state we first estimate the mean position in the synchronised state of the respective system. Then we start the computation with a very narrow Gaussian function located in the determined position as the initial probability density p(x, t = 0).
Effective diffusion coefficient: In the mode locked states the limit of equation 12 could not be reached within the time span of our computations due to intermediate plateaus as shown in figure 3(a). In these cases the linear increase of the plateaus for 500 oscillations was calculated and used to determine D eff .