Emergence of traveling waves in linear arrays of electromechanical oscillators

Traveling waves of mechanical actuation provide a versatile strategy for locomotion and transport in both natural and engineered systems across many scales. These rhythmic motor patterns are often orchestrated by systems of coupled oscillators such as beating cilia or firing neurons. Here, we show that similar motions can be realized within linear arrays of conductive particles that oscillate between biased electrodes through cycles of contact charging and electrostatic actuation. The repulsive interactions among the particles along with spatial gradients in their natural frequencies lead to phase-locked states characterized by gradients in the oscillation phase. The frequency and wavelength of these traveling waves can be specified independently by varying the applied voltage and the electrode separation. We demonstrate how traveling wave synchronization can enable the directed transport of material cargo. Our results suggest that simple energy inputs can coordinate complex motions with opportunities for soft robotics and colloidal machines. For small scale biological systems such as cilia, movement is achieved by rhythmic motor patterns that organize spontaneously within arrays of driven oscillators. The authors show that conductive spheres oscillating between biased electrodes create similar traveling wave motions which can be used to direct the transport of cargo.

T raveling waves of mechanical actuation provide a versatile strategy for locomotion and transport in both natural [1][2][3] and engineered 4-8 systems across many scales. In vertebrates such as the aquatic lamprey 1,9 , these and other rhythmic motor patterns are orchestrated by networks of neurons called central pattern generators (CPGs) 10 , which are often idealized as systems of coupled oscillators 1,9 . The rhythmic output of these oscillators is relayed to actuators (e.g., muscles) to produce complex motions without the need for sensory feedback. Similar control strategies based on CPGs are used to direct locomotion in macroscopic robots 4 . At smaller scales, however, it becomes increasingly challenging to accommodate centralized control systems capable of directing the coordinated actions of multiple actuators. Instead, microorganisms such as ciliated protozoa integrate pattern generation and mechanical actuation within a single material system. The oscillatory motions of beating cilia couple to one another through hydrodynamic interactions to produce metachronal waves that drive cellular locomotion through viscous surroundings 2,11,12 . This biological example illustrates how the coupled motions of many mechanical oscillators can organize spontaneously and autonomously to perform dynamic functions.
The realization of synthetic systems that mimic such functions requires experimental strategies for powering mechanical oscillators and for coupling their motions to achieve the desired dynamics. One approach relies on coupling reaction-diffusion patterns to the mechanical deformation of responsive gels, for example, to achieve traveling wave motions in excitable media 7,8 . Despite fascinating demonstrations of this approach on millimeter length scales, it remains challenging to miniaturize due to the need for faster reactions that compete with diffusion at smaller scales. To achieve scalable mechanisms of pattern formation, the processes that drive oscillations should scale in the same way as those used to couple neighboring oscillators. In this context, electromechanical oscillators based on contact charge electrophoresis (CCEP) 13,14 can provide a useful model on length scales spanning millimeters 15 to microns 16 (perhaps even nanometers 17,18 ).
CCEP refers to the back-and-forth motion of a conductive particle through an insulating fluid separating two electrodes subject to a constant voltage. The particle charges on contact with either electrode and moves down the applied potential gradient, thereby transporting charge between the biased electrodes. This type of electromechanical oscillator is fundamentally distinct from the weakly damped harmonic oscillators of microelectromechanical systems (MEMS) 19,20 , which rely on resonant excitation by time-varying fields. By contrast, CCEP oscillators are powered by a constant thermodynamic driving force and operate even under conditions of strong damping, which arise at small scales and in viscous environments. Similar to those of molecular motors 18,21 , CCEP motions can be rectified to perform mechanical work or to transport material cargo 18,22 . Moreover, the charge acquired by the particle and the forces driving its motion are well described by classical electrostatics, which is invariant to changes in scale. The discovery of new CCEP motions at the macroscale is therefore transferable to emerging applications at the microscale 13 .
Here, we investigate the collective dynamics of many CCEP oscillators positioned along a linear array between two (nearly) parallel electrodes (Fig. 1a). Each oscillator is comprised of a conductive sphere that moves back and forth between the electrodes along a dielectric track. Oscillatory motions are driven by the repeated charging of the particles on contact with either electrode and their subsequent movement in the applied field. The dynamics of neighboring oscillators are coupled to one another through the electrostatic interactions between the charged particles. We show how this electrostatic coupling mediates the organization of phase-locked states in which all oscillators move with a common frequency. Interestingly, the distribution of oscillator phases at steady state corresponds to traveling waves of particle motion with a characteristic wavelength comparable to the electrode separation. These experimental observations are explained by a Kuramoto-like model 23,24 that accounts for weak repulsive coupling between neighboring phase oscillators and for small systematic variations in their natural frequencies. We demonstrate how traveling wave synchronization can be used to transport material cargo along the length of the oscillator array. More generally, our approach shows how simple energy inputs can power complex patterns of mechanical actuation, which may be useful in powering the motions of soft robots 25,26 and colloidal machines [27][28][29][30][31] .

Results
Experiment. Copper spheres (a = 1 mm in radius) immersed in mineral oil were positioned within an array of dielectric tracks connecting two plane electrodes separated by a distance L (Fig. 1a). The tracks were spaced evenly with a period W = 3a and aligned perpendicular to the electrode surfaces and to the direction of gravity. Each track contained a single particle, which was free to move back and forth between the two electrodes. Application of a constant voltage (typically, V = 10 kV to 20 kV) caused the particles to oscillate continuously between the electrodes via CCEP 13,14 . The conductive particles acquired an electrostatic charge on contact with the biased electrodes and moved under the influence of the applied field. This periodic cycle of contact charging and electrostatic actuation continued for as long as the voltage was applied. Figure 1b shows the reconstructed trajectory of a single sphere oscillating between the two electrodes. Each time the particle contacts an electrode, its charge changes sign and the particle reverses direction under the influence of the field. To facilitate the analysis of multiple particles over many oscillation cycles, we record the times at which each particle contacts one of the electrodes. From this data, we approximate the phase of each oscillator by interpolating between successive contacts as φ(t) = 2π(t − t k )/(t k+1 − t k ) where t k denotes the time of the kth contact and t k ≤ t < t k+1 . By definition, the oscillator phase increases at a constant rate equal to the natural frequency, dφ/dt = ω, which is approximated by averaging over many oscillation cycles as ω = 〈2π/(t k+1 − t k )〉 k . Repeating this analysis for each particle in isolation (i.e., one track at a time), we observed small systematic variations in the natural frequency ω n with respect to the oscillator position n along the array (Fig. 1c, solid markers). The spatial gradients in the oscillator frequency were caused by small deviations in the electrode alignment, which was controlled only to within ca. 1°of parallel. Particles oscillated faster where the electrodes were closer together due to an increase in field strength at those locations.
Despite variations in their natural frequencies, linear arrays of N particle oscillators moving simultaneously evolved in time to a phase-locked state, in which each particle moved with a common frequency (Fig. 1c). Interestingly, the synchronized particles did not move in phase with one another but rather organized to form a single traveling wave, which remained stable for hundreds of oscillation cycles ( The stability of the synchronized state and the distribution of oscillator phases therein depended on the number of oscillators N in the array. For N = 2 oscillators, the particles moved in antiphase as reported previously 15 (Fig. 2a). Such antiphase synchronization suggests that neighboring oscillators are coupled together by repulsive interactions such as the Coulombic forces between like-charged particles. As the number of particles was increased, the average phase difference between successive oscillators decreased giving rise to stable traveling waves with wavelengths spanning many oscillators ( To better understand the synchronized state, we quantified the distribution of oscillator phases for N = N* as a function of the electrode separation L and the oscillator spacing W. For each electrode configuration, we started with N > N* oscillators and removed one particle at a time from the end of the array until reaching a stable synchronized state, at which the phase difference between neighboring oscillators was constant in time. Figure 2d   difference χ ss n ¼ φ nþ1 t ð Þ À φ n t ð Þ t at the stationary state for a typical experiment. The phase difference was smallest (in magnitude) near the center of the array and largest near the edges. For each such state, we defined a characteristic wavelength in terms of the average phase difference as λ ¼ 2πW= χ ss n n . This wavelength increased linearly with the electrode spacing L but was independent of the oscillator spacing W over the range explored (Fig. 2e).
Minimal model of traveling wave synchronization. The experimental observations are reproduced by a Kuramoto-like model 23 that accounts for the local repulsive coupling between neighboring oscillators and the systematic variations in their natural frequencies. In the model, we adopt the following simplified description of CCEP dynamics 18 . On contact with either electrode, a conductive sphere acquires a charge q ¼ ± 2 3 π 3 εa 2 E, where E is the electric field at the electrode surface, and ε is the permittivity of the surrounding dielectric. This expression-first derived by Maxwell 32 -assumes that charge flows to/from the particle until its potential equals to that of the contacting electrode 33 . The electrostatic force on the particle is approximated as F = qE, which drives motion with velocity U = F/γ, where γ is a constant friction coefficient. Figure 3a shows how the charge q and position h of a single oscillator depend on its phase φ = ω o t, where ω o = πq o E o /γL is the natural frequency defined in terms of the applied field E o = V/L and the Maxwell charge q o ¼ 2 3 π 3 εa 2 E o . By comparing the measured frequency in Fig. 1c to the prediction of the model, the friction coefficient can be estimated to be γ = 1.8 × 10 −3 N s m −1 , which is ca. 4 times larger than the Stokes drag, γ s = 6πηa. The increased drag is attributed to the solid boundaries formed by the patterned tracks and the planar electrodes ( Supplementary Fig. 5) 34 .
The presence of neighboring oscillators influences both the charge that a particle acquires and the speed at which it moves. To describe these interactions, we decompose the electric field as E = E o + E′, where E o is the applied field and E′ is a disturbance field due to neighboring particles, which are approximated as point charges (see Methods). In the limit of weak interactions (i.e., when E ′ ( E o ), the moving particles are well approximated as phase oscillators with weak repulsive coupling between nearest neighbors. The phase of the nth oscillator evolves in time as where ω n is the natural frequency, and the function f( ) describes the phase-averaged interactions between neighboring oscillators as a function of their phase difference χ n = φ n+1 − φ n . The boundaries of the array are open such that oscillators at the edges (n = 1, N) interact with only one neighbor 35 .
We assume a uniform gradient in the natural frequency: where ω o is the mean oscillator frequency, and Δ is the frequency difference between successive oscillators due to a small angle θ between the electrodes Δ=ω o % 3ðW=LÞθ ( 1 ð Þ . Interestingly, this model was investigated previously as a possible explanation for traveling wave oscillations in the CPG of the aquatic lamprey 1 . In the present context, the primary interaction between neighboring oscillators is electrostatic in origin; other interactions are neglected. In particular, the neglect of hydrodynamic interactions was supported by experiments in which neighboring particles were separated by solid walls without altering their collective dynamics (Supplementary Fig. 6). Approximating the particles as point charges, we compute the electrostatic interaction averaged over one oscillation cycle for a constant phase difference χ (see Methods). These repulsive interactions are described by an odd function characterized by the location χ max and height f max of its maximum (Fig. 3b). For large electrode separations L ) W ð Þ , these quantities are well approximated as f max % 1 2 ffiffi 3 p ðq 2 o =εW 2 LγÞ and χ max % π ffiffi 2 p ðW=LÞ. Our assumption of weak coupling implies that the phase velocity due to interactions is small compared to the natural frequency-that is, f max ( ω o or, equivalently, a=W ( 0:72. Additional simulations incorporating the full amplitude dynamics provide further support for the phase oscillator approximation under the experimental conditions of a/ W = 0.33 (Supplementary Note 1 and Supplementary Fig. 7).
The competition between the repulsive interactions and the frequency gradient leads to stable stationary solutions described by f ðχ n Þ ¼ À 1 2 ΔnðN À nÞ (see Methods). This solution exists provided that the number of oscillators is below some critical value N Ã ¼ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 8f max =jΔj p . Figure 3c-e shows the stable solution in terms of the phase difference and the sine of the phase for N = 15 oscillators-just below the chosen critical value of N* = 15.5. The addition of more particles (N > N*) causes the waves to break periodically in the center of the array (Fig. 3f). Physically, faster oscillators pile up behind the slower ones and are prevented from passing by the local repulsive interactions. In this way, the frequency gradient acts to compress the oscillator phases together to create longer waves that travel always from slower to faster oscillators. When compression by the frequency gradient exceeds the repulsive barriers between neighboring oscillators, global synchronization is lost and the waves break.
At the critical oscillator number (N = N*), repulsive interactions are at their maximum (χ ss n $ χ max ), and the characteristic wavelength scales linearly with the electrode separation in agreement with experimental observations (Fig. 2e)-that is, λ $ 2πW=χ max $ L for W ( L. Moreover, the critical oscillator number observed in the experiment implies a certain angle between the electrodes, which can be estimated from the model as For the conditions of Fig. 2, this angle is predicted to be θ = 1.1°, which agrees well with that measured from the experimental images (Fig. 1b). Smaller angles allow for stable waves containing more particles.
The average phase within the wave evolves in time at a constant rate equal to the average frequency ω o , which is specified independently of the wavelength. In the experiment, the oscillator frequency could be altered by changing the applied voltage; however, the range of accessible frequencies was limited by dielectric breakdown at higher voltages and by particle sticking at lower voltages 14 . Notably, the frequency of the phase-locked state was slightly faster than the average natural frequency (Fig. 1d). Dipolar interactions among the particles (neglected here) break the odd symmetry of the interaction function thereby altering the frequency of the synchronized state.

Discussion
We have shown how arrays of electromechanical oscillators can organize spontaneously to form synchronized traveling waves of particle motions powered by a constant input voltage. The direction of wave propagation is determined by small gradients in the natural frequencies of the oscillators and can be controlled by introducing a small angle between the otherwise parallel electrodes. The characteristic wavelength is approximately equal to the electrode separation L and corresponds to the largest possible wave that can be stabilized by repulsive interactions among the charged particles. The traveling wave motions are robust to perturbations and can be harnessed to direct the transport of material cargo. In particular, Fig. 4a shows how traveling waves can direct the motion of gas bubbles floating at the interface just above the oscillating particles (see also Supplementary Movie 6).
Bubbles are transported in the direction of wave propagation at speeds comparable to the wave velocity (Fig. 4b). In contrast to previous strategies for rectifying CCEP motions based on ratcheted channels 22 or asymmetric particles 16 , the present approach relies on the self-organization of multiple particles working in concert.
Beyond this initial demonstration, traveling wave synchronization of CCEP oscillators may be useful for peristaltic pumping within microfluidic systems. Unlike standard pressure pumps, those based on traveling wave motions allow for recirculating flows 36 and would complement existing applications of CCEP in microfluidic cargo transport 22,37 , separations 22 , and fluid mixing 38 . These CCEP-powered unit operations rely on constant voltages at low input power, which makes them attractive for use in portable, battery-powered microfluidic devices 13 . One important limitation of these oscillators is their reliance on the dielectric environment provided by non-polar fluids; CCEP motions cannot be sustained in even weakly conductive liquids such as deionized water 38 . However, recent advances in soft robotics suggest one strategy for circumventing this limitation by encapsulating non-polar liquids in stretchable, impermeable compartments 26 . These soft composite materials can be deformed by applying voltages to stretchable electrodes patterned on their surfaces. By incorporating arrays of CCEP oscillators within such dielectric compartments, it should be possible to create self-organized motions that drive transient deformations and thereby locomotion of soft robotic materials.

Methods
Experiment set-up. Periodic arrays of dielectric tracks were 3D printed in acrylonitrile butadiene styrene (ABS) with a period of W = 3-5 mm. The array was sandwiched between two copper plates separated by a distance L = 10-30 mm (McMaster-Carr 3350K201) and immersed in mineral oil (Sigma Aldrich CAS No. 8042-47-5). The tracks were aligned perpendicular to the electrode surfaces and to the direction of gravity. Each track contained a single copper sphere (McMaster-Carr No. 64715K16, radius a = 1 mm), which rolled freely between the two electrodes (Fig. 1a). Prior to use, the system was heated in a 60°C oven for several hours to remove any moisture. The copper electrodes were connected to an amplifier (Trek 20/20C) with an output voltage V = 0-20 kV. The particles were illuminated from below by a light emitting diode (LED) and their motions captured by a high-speed camera (Phantom V310). During each experiment, the electrodes were energized to a specified voltage and the resulting particle motions captured. The voltage was switched off for at least 1 min between successive experiments to allow for the dissipation of any space charge accumulated on the surfaces of the tracks and/or the electrodes. Particle location data were extracted using standard image tracking routines in MATLAB.
Bubble transport. Bubbles were generated within the oscillator array by a continuous flow of air supplied by a syringe pump at a rate of 0.2 ml min −1 . The air was delivered through a tube that connected to a hole in the base of one of the tracks (Fig. 4a). The bubbles (ca. 4 mm in diameter) were transported down the length of the array by the coordinated motions of N = 17 particles of radius a = 1.5 mm (Supplementary Movie 6). In these experiments, the applied voltage was V = 18 kV, the electrode separation was L = 20 mm, and the height of the mineral oil above the base of the track was 4 mm. Bubbles were always transported in the direction of wave propagation, which was controlled by introducing a finite angle between the two electrodes. Control experiments with no applied voltage showed no bubble motion in either direction.
Electrostatic interactions. We first consider a single point charge q positioned at a height z = h between two grounded planes at z = 0 and z = L. The resulting electrostatic potential at point x is where r ¼ ffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi is the radial distance from the charge, and K 0 ( ) is the zeroth order modified Bessel function of the second kind. For large arguments, the Bessel function decays exponentially as K 0 ðsÞ ! e Às ffiffiffiffiffiffiffiffiffi π=2s p ; the infinite series can be truncated for some m ) L=πr to obtain an accurate approximation. The corresponding electric field in the z-direction is given by E z = −∂ϕ/∂z.
We now consider how the disturbance field E′ due to one oscillator i influences the dynamics of another oscillator j in the limit of weak coupling. At zeroth order in E′, the phase of each oscillator increases at a constant rate ω o such that φ i = ω o t and φ j = ω o t + χ, where χ = φ j − φ i is the constant phase difference. The charge q = q(φ) and position h = h(φ) of each oscillator depend on the phase as shown in Fig. 3a. At first order in E′, the disturbance in the phase of oscillator j evolves as where the factor πn(φ j )/Lγ relates the electric force to the corresponding phase velocity with n(φ j ) = ±1 indicating the direction of travel. The bracketed terms describe two types of electrostatic interactions. First, the disturbance field due to particle i drives particle j to move faster or slower between the electrodes. Using the point charge solution (2), this disturbance E′(φ i , φ j ) is given by where W is the oscillator spacing. Second, the disturbance field due to particle i alters the charge acquired by particle j on contact with either electrode; the disturbance charge q′(φ i , φ j ) is given by Physically, the charge on particle j is determined by its most recent contact with either electrode (φ j = 0 or π); the field due to particle i at the time of that contact determines the disturbance charge.
We can now average these two interactions over one oscillation cycle to obtain the phase-averaged interaction function, Carrying out the integration, each of the two electrostatic interactions produce contributions of the same mathematical form with the second term contributing twice that of the first, This final expression is plotted in Fig. 3b for the case of W/L = 0.12.
Stationary solution. Starting from Eq. (1), we recast the oscillator dynamics in terms of the phase difference χ n = φ n+1 − φ n and the average phase Φ ¼ 1 N P n φ n . Taking the difference in the phase dynamics of successive oscillators, we obtain the following equation for the phase difference: which makes use of the fact that f( ) is an odd function. At the open boundaries of the array, the phase difference evolves as In addition to these N − 1 equations, the dynamics of the N oscillators is described by the that of the average phase, ∂Φ/∂t = ω o , which is fully decoupled from the phase differences. Setting the time derivatives in Eqs. (8) and (9) equal to zero, the resulting recurrence equation can be solved to obtain the stationary solution, f ðχ n Þ ¼ 1 2 Δnðn À NÞ, presented in the main text. This solution exists provided that f max > 1 8 N 2 jΔj (Supplementary Note 2) and is stable when f′(χ n ) < 0 (Supplementary Note 3). The characteristic relaxation time for approaching the stationary state is given by the diffusive-like scaling relation τ $ χ max N 2 =π 2 f max .

Data availability
The data that support the findings of this study are available from the corresponding author upon reasonable request.