Metachronal waves in magnetic micro-robotic paddles for artificial cilia

Biological cilia generate fluid movement within viscosity-dominated environments using beating motions that break time-reversal symmetry. This creates a metachronal wave, which enhances flow efficiency. Artificially mimicking this behaviour could improve microfluidic point-of-care devices, since viscosity-dominated fluid dynamics impede fluid flow and mixing of reagents, limiting potential for multiplexing diagnostic tests. However, current biomimicry schemes require either variation in the hydrodynamic response across a cilia array or a complex magnetic anisotropy configuration to synchronise the actuation sequence with the driving field. Here, we show that simple modifications to the structural design introduce phase differences between individual actuators, leading to the spontaneous formation of metachronal waves. This generates flow speeds of up to 16 μm/s as far as 675 μm above the actuator plane. By introducing metachronal waves through lithographic structuring, large scale manufacture becomes feasible. Additionally, by demonstrating that metachronal waves emerge from non-uniformity in internal structural mechanics, we offer fresh insight into the mechanics of cilia coordination. Biological cilia are extremely effective for fluid flow, due to their beating motion creating metachronal waves. Here, metachronal waves are created in micro-robotic paddles that rotate in response to a magnetic field, creating fluid flow speeds of up to 16 μm/s.

M agnetic micro-robots that can be actuated within fluids have grown into a thriving field of research due to their potential medical applications in drug delivery and microsurgery [1][2][3] . Due to the microscale dimensions of these robots, viscosity dominates over inertial motion, a regime characterized by a low Reynolds number. Within this regime, fluid flow around a fixed micro-robot depends solely on the forces that the micro-robot applies at any instant, not on the rate of force application. Consequently, fluids return to their unperturbed state if micro-robot actuation has time-reversal symmetry 4 : there is no net movement of fluid and no net displacement of the microrobot. In order to introduce asymmetry into micro-robot actuation, many researchers have found inspiration from strategies employed by biological organisms. Mimicking strategies employed in biological microswimmers [5][6][7] , free-floating robots have been demonstrated to self-propel by rotating helical flagella 1-3 , generating travelling waves along flexible flagella [8][9][10] , controlling complex interactions between a colloidal swarm 11,12 , or asynchronously modulating shape during actuation 13 . Substrate-bound micro-robots frequently model biological cilia in order to generate fluid movement [14][15][16][17] .
Cilia are hair-like microstructures that protrude from the outer surface of biological cells and flex to move fluids and biomolecules along cell walls. Individually, cilia perform complex cyclical movements to break time-symmetry, combining a rotatory motion with an alteration of shape between the straight "effective" stroke and the bent "recovery" stroke 18,19 . Despite a lack of any internal coordination system, phase differences in the oscillatory actuation of neighbouring cilia generate pulses, called metachronal waves, that travel along the cell surface to enhance the efficiency of the generated flow 19,20 . Theoretical models indicate that the metachronal waves emerge from hydrodynamic interactions between cilia [21][22][23] .
Artificially reproducing the properties of cilia could generate more efficient flow generation in microfluidic applications. Several mechanisms have been proposed to create artificial cilia, including inducing electrostatic bending 24 , illumination of photoactive polymer flaps 25 , controlling acoustic resonance 26 and pneumatic manipulation 27 . By far the most studied mechanism is magnetic actuation, in which columns typically formed from superparamagnetic particles 15,16 or a polymer/magnetic particle composite 14,17,28,29 are distorted by the presence of a magnetic field. While all these schemes mimic cilia in appearance, only a few are capable of producing metachronal waves as they actuate all cilia with the same phase. These schemes impose metachronallike actuation under a uniform magnetic field by setting the magnetizations of individual cilia during fabrication. Use of millimetre-scale designs enables either the cilia 28 or the substrate 29 to be flexed prior to saturating the magnetization along a defined direction. Alternatively, several stages of magnetization alignment may be performed during microfabrication to effectively create a spatially varying anisotropy profile over the cilia array 30 . Each of these strategies predetermines the form of the metachronal wave according to the field strength and direction, while the speed of the wave is fixed by the field frequency.
A simpler approach to generating metachronal actuation is to modify the cilia length across the array 31 , although this has the drawback of also spatially varying the hydrodynamics of the component producing flow. Here, we take a different approach, modulating a part of the structure that influences actuation, but which produces negligible flow in its own right. We demonstrate that spatially varying the axle width of paddle-based magnetoelastic micro-robots enables metachronal waves to emerge spontaneously as the paddles rotate between bistable states defined by the field direction.

Results
Mechanical properties of paddle membranes. To establish the mechanical properties of the micro-robots, they were attached across a polydimethylsiloxane (PDMS) frame of 5 mm inner diameter and imaged in water under a uniaxial magnetic field. The micro-robots consisted of 10-µm-thick PDMS membranes divided into 650 µm × 500 µm cells, so that approximately 50 freestanding cells filled the 5 mm diameter frame. Each cell was composed of four 120 µm circular paddles with centre-to-centre spacing of 150 µm (Fig. 1a). Each paddle contained a hexagonal array of 6-µm-thick, 10 µm diameter Co-Ni-P cylinders with a pitch of 20 µm and was joined to a stiff frame by an axle of width w = 3-12 µm. Samples were imaged in a well filled with water to a height, h, above the membrane (Fig. 1b) with the axles aligned along the y-axis co-ordinate (Fig. 1a). Torque from the magnetic field (applied perpendicular to the axles, along the x-axis) caused the paddles to rotate on their axles.
Before considering how metachronal action can be introduced, we first demonstrate the basic paddle actuation. In cells containing four identical 3 µm axles, paddles rotated in unison due to the torque from a 15 Hz  Co-Ni-P particles appeared as bright spots prior to the forward stroke ( Fig. 2a), but were suppressed afterwards (Fig. 2e). This shows that the paddles had two stable orientations, at 0°under forward field (in the relaxed state) and rotated by almost 180°u nder reversed field (the metastable state). Figure 2k demonstrates that the paddle remained in the metastable state until a forward field was applied, indicating that the timescale of rotation due to elastic torsion, counteracted by fluid drag, was slow compared to the field frequency.
Each paddle displayed a systematic sense of rotation, although there was some variation in the sense of rotation between different paddles. Stray field interactions between the paddles would favour an opposite sense of rotation between neighbours, so may contribute to the rotation bias, but cannot solely be responsible for the pattern seen in Supplementary Video 1. Noting that any deviation from perfect parallel alignment between the paddle magnetization and the field would be sufficient to bias the rotation sense, we suggest that differences in rotation sense were introduced either by slight variations in out-of-plane magnetization texture or minor distortions in the membrane structure.
Competition between elastic torsion and the minimal magnetic torque as rotation approached 180°prevented complete symmetric rotation. Nevertheless, all fields achieved near-180°r otation except at the lowest fields (below 4 kA/m), which had slightly smaller maximum rotation (Fig. 2l). No change in rotation amplitude occurred under field frequencies between 1 and 35 Hz (at 11.2 kA/m, Supplementary Fig. 2), indicating that resonance does not play a role in actuation. Paddles on axles wider than 3 µm exhibited similar trends (Fig. 2l). However, the reduction in the maximum rotation achieved at fields below 4 kA/ m became more pronounced as the axle width increased. Such reduced actuation was a consequence of the increased elastic torsion in wider axles. Therefore, while rotation was driven by the magnetic field, the structural design of the membrane provided an independent mechanism through which the actuation response could be modified.
Generation of metachronal actuation. To understand how timereversal symmetry of paddle rotation could be broken to produce metachronal waves and therefore generate flow in a low Reynolds number environment, we developed a simplified model of the membrane, consisting of periodically repeating two-paddle cells ( Fig. 3a and Supplementary Note 1). In the model, the membrane frame connecting the paddles was modelled as a series of nonmagnetic beads (Fig. 3a, B1-4) joined by inflexible links to form a backbone structure. Each paddle was represented by two hard magnetic beads (Fig. 3a, P1-4) connected by inflexible links about a pivoting axis at B2 or B3. All links in the model were elastically stiff, producing negligible change in length under tension or compression. Building on our previous models of magnetically actuated swimmers 32-35 , we modelled the movement of the membrane in a fluid, but pinned its translational motion. Due to the restricted translational freedom, actuation that would otherwise produce propagation instead resulted in flow 35,36 . Torques produced by a magnetic field caused the magnetic beads to rotate around their respective paddle axes, with resultant flow due to each bead determined by Stokes' law (the hydrodynamic resistance of the links was neglected). To describe the mechanical twisting of the paddle axles, rotation of each paddle around its pivot away from its rest position was resisted by a restoring elastic torsion proportional to a (non-dimensional) rotational stiffness. Each paddle was modelled with an individual rotational stiffness, l 1 and l 2 , enabling the paddles to respond independently to the applied field.
Despite its simplicity, the model was able to reproduce experimental features of actuation and to predict design paradigms for generating flow in a low Reynolds number environment. When both paddles had identical rotational stiffness, the maximum paddle rotation depended on the field amplitude, in agreement with the mechanical experiments (Fig. 2l). Additionally, the magnitude of the rotational stiffness also influenced the largest actuation achieved ( Supplementary  Fig. 3). However, no flow was generated when the paddles had the same rotational stiffness (Fig. 3b), since the paddle actuation traced a path that was reciprocal in time. To produce a flow, it was necessary for the rotational stiffness of the paddles to differ (Fig. 3b).
Since elastic torsion is proportional to rotational stiffness, the two paddles with different rotational stiffness experienced different overall torque, even under identical magnetic torque. This led to phase differences between the time taken for each paddle to achieve a given rotation angle ( Fig. 3c and Supplementary Video 2). As the magnetic torque is highly nonlinear with time (being dependent on both the paddle rotation and a time-dependent field amplitude), small initial differences in torque were amplified to produce large phase difference between paddles. Fluid interactions between paddles and spring interactions via the backbone support had little causal effect on the emergence of a phase difference between paddles. Time symmetry between forward and reverse fields was broken because the field reversed direction while the paddle remained at near-180°r otation, causing rotation to be magnetically driven during the recovery stroke as well as the forward stroke (Supplementary Video 2). Whereas elastic torsion acted against the magnetic torque in the forward stroke, it acted in the same direction as the magnetic torque during the recovery stroke. Therefore, the phase differences between paddles were enhanced during the recovery stroke.
Analysis of the conditions under which the modelled paddle rotation reproduced the experimental measurements indicated that rotational stiffness was proportional to the axle width ( Supplementary Fig. 3). Therefore, the model predicted that a graded axle width profile would be sufficient to enable flow generation. To test this experimentally, the actuation of a membrane containing cells with axles of width w = 3, 6, 9 and 12 μm was imaged under a 15 Hz, 1.9 kA/m field. During both the forward (Fig. 4a) and reverse (Fig. 4b) field cycles, the variation in the mechanical properties of the axles ensured that the paddles rotated sequentially to produce a metachronal wave throughout the membrane (Supplementary Video 3). The order of actuation was similar in the forward and reverse cycles, with the metachronal wave initiating at the 12 μm axle and propagating along the cell towards the 3 μm axle. Whereas paddles on identical axles rotated in unison (Fig. 2k), the variation of the rotational stiffness between the paddles in the graded axle width membrane introduced phase differences in the rotation angles of the paddles during actuation (Fig. 4c). To quantify the phase differences between the different paddles, we assessed the time delay to reach a specific angle, 90°(when returning to 0°) and found that the time delay decreased systematically with increasing axle width (Fig. 4d). Inhomogeneities in the sample meant there was a spread in the phase of paddles on nominally identical axles in different cells across the membrane, with around 1.5 ms standard deviation in the time delay for each axle width ( Supplementary Fig. 4). Nevertheless, each individual cell followed a similar actuation sequence, such that Fig. 4 is representative of the average trend. Therefore, use of a graded axle width profile broke time-reversal symmetry in the membrane to produce non-reciprocal actuation.
Production of fluid movement. Having established the mechanism of actuation in the membranes, we now consider the effect of the actuation on the fluid around the membrane. Fluid motion generated by the graded axle width profile membrane (w = 3, 6, 9 and 12 μm within each cell) was measured from the motion of polystyrene tracer beads on the water surface (Supplementary Video 4). Whereas negligible bead motion occurred in the absence of an applied field ( Fig. 5a and b), actuating the membrane generated persistent movement (Fig. 5c-f). Net flow (calculated from the average translation vector of the beads) was  Fig. 5h) fields generated much faster flow than occurred prior to activation due to random drift. Since the enhanced flow speed was statistically significant (using a onetailed paired t-test), it was unlikely to be caused by random fluctuations in surface flow, such as may be caused by temperature gradients, air currents or other environment factors. Indeed, the generation of a significant flow along the direction of decreasing axle width was in agreement with the modelled effect of inducing non-reciprocal actuation via a graded rotational stiffness profile (Fig. 3b).
Flow speed is a key performance indicator for many microfluidic applications, since it determines how fast fluid samples can be processed. We therefore measured how the net flow speed could be controlled using the amplitude and frequency of the applied field across a range of fluid heights. At a constant 15 Hz frequency, net flow speed was independent of the field amplitude (Fig. 5i), averaging 7.3 ± 2.5 μm/s. Since the field strength only alters the speed of disc rotation, not the sequence of actuation, the insensitivity demonstrates that the membrane actuated within the low Reynolds number regime. By contrast, increasing the field frequency generated progressive enhancement of the flow speed, following a linear trend up to 20 Hz (Fig. 5j). Despite not modifying the actuation sequence, proportional increases in flow speed with field frequency are consistent with low Reynolds number operation, assuming a similar amount of fluid was moved during each actuation cycle. Indeed, calculating the Reynolds number (Re) from the tangential velocity of the paddles 13 , flow was dominated by viscosity (Re < 1) up to 20 Hz, coinciding with the region of linear frequency dependence. At higher frequencies, the effects of inertia manifested in greater-than proportional increases in flow speed (Fig. 5j). Both Fig. 5i and j indicate that the flow speed produced by the membrane was largely unaffected by the fluid height, possibly due to the free surface at the top of our well (Fig. 1b). Therefore, the membrane is capable of generating significant flows within the low Reynolds number regime at heights compatible with microfluidic applications.

Discussion
Several features displayed by our membranes have potential uses in microfluidic point-of-care devices. The ability to generate flows with speeds up to 16 μm/s by tuning the field frequency potentially provides a mechanism to portion out fluid measures. Furthermore, insensitivity of flow speed to the fluid height above the membrane means that performance is tolerant towards variations in the volume of sample fluid provided by the user. This could facilitate the distribution of test samples among an array of biosensors for parallel diagnosis of multiple diseases.
On a more fundamental level, by tailoring the structural (elastic) profile of the axles within our membrane, we demonstrated an architecture in which metachronal waves were produced spontaneously. As such, the wave-generation mechanism was distinct from previous strategies of mimicking cilia motion, which either rely on imposing a pre-determined actuation sequence through the magnetic anisotropy profile 37,38 , or induce phase differences through progressively restricting movement of cilia along an array 31 . Considering the simplicity with which the structural profile can be controlled in a lithographically fabricated system, our method of producing metachronal waves could feasibly be adapted to induce complex flow patterns, such as a curved flow path, which would otherwise be difficult to achieve. Our membrane architecture could also provide insight into the spontaneous emergence of metachronal waves in biological cilia 20 , which has previously been attributed solely to the effects of hydrodynamic synchronization 39 . Given that our artificial system only displays metachronal waves when there is variation in the elastic properties of the actuating elements, our results highlight that internal structural mechanics may play a previously unrecognized role in the spontaneous generation of biological metachronal waves.

Methods
Sample fabrication. Magneto-elastic composite membranes were fabricated in three stages 40 . Initially, magnetically hard Co-Ni-P particles were deposited onto a glass/PMMA/Au substrate using photolithography with a 10-µm-thick AZ9260 resist and electrodeposition in an aqueous solution of 60 g/L of cobalt (II) sulfate heptahydrate (CoSO 4 ·7H 2 O), 60 g/L of nickel (II) sulfate hexahydrate (NiSO 4 ·6H 2 O), 3.9 g/L of sodium hypophosphite monohydrate (NaPO 2 H 2 ·H 2 O) and 40 g/L of ammonium chloride (NH 4 Cl) at 30°C. Secondly, a PDMS membrane was patterned by performing a further photolithographic step and spinning PDMS (Sylgard 184 PDMS base mixed in a 10:1 ratio with curing agent) at 7500 r.p.m. for 150 s. To maximize elasticity, the PDMS was cured at room temperature over 2 days 41 . Excess PDMS above the height of the photoresist was removed through reactive ion etching with SF 6 (10 sccm) and O 2 (2 sccm) under a 20-W power and 30 mTorr pressure. Finally, liberation of the membrane was achieved by using a nickel-compatible gold etchant to expose the sacrificial PMMA layer beneath the sample, adhering a supporting frame to the membrane with uncured PDMS (cured in situ) and soaking in acetone to dissolve the sacrificial PMMA layer. Residual gold on the membrane, left due to masking from the membrane itself, was removed by a second wash in nickel-compatible gold etchant. To define a consistent magnetization direction across the sample, the liberated structure was magnetized in an 80 kA/m field.
Actuation and flow speed measurement. For actuation measurements, samples were placed into 20 mm diameter well filled with water and imaged using a microscope incorporating a uniaxial Helmholtz coil that generated fields of 1.9-17.4 kA/m at frequencies of 5-35 Hz. Water height was calculated from the difference in focal position between the membrane and water surface. To resolve paddle rotation, membrane actuation was captured at 1500 frames per second (fps). Rotation angle was measured trigonometrically, as the inverse cosine of the projected paddle width during rotation, normalized to the maximum paddle width. Flow speed was measured using 15 µm diameter polystyrene beads that were added to the water surface. Bead motion was recorded at 200 fps, both prior to and during actuation. Net flow speed during field application was calculated as the magnitude of the average displacement vectors of 58-187 beads across the imaging area over a period of 15 s. Similarly, average flow prior to field application was measured from the drift of 50 beads over 4 s, enabling the net flow speeds to be drift-corrected. All measurements from images were performed manually using the ImageJ software.
Statistical analysis. To compare flow before and during actuation, the average flow speed from 50 beads was measured on n = 3 separate occasions under identical conditions (1.9 kA/m or 17.4 kA/m field at 15 Hz, h = 544 µm). On each occasion, flow was measured continuously from 4 s immediately before field application to 4 s after the field had been applied (on the same beads), ensuring that each measurement was independent and included a similar drift component before and during actuation, enabling a paired statistical test. Since induced flow could not be slower than the drift, all statistical comparisons were performed using a one-tailed paired t-test; differences between means were considered significant when p < 0.05.
Theoretical model. To investigate the mechanism of flow generation in the magneto-elastic membranes, we developed a minimal theoretical model based on our previous models of magnetic swimmers [32][33][34][35] . This is described in detail in Supplementary Note 1. Briefly, the model consists of eight beads within a unit cell, connected together via links in the arrangement shown in Fig. 3. Beads B1-4 form a backbone structure, with beads P1 and P2 connected to the backbone at bead B2 and similarly, beads P3 and P4 connected to bead B3. Therefore, the bead triplets P1-B2-P2 and P3-B3-P4 form two rotating paddles on the backbone structure. Beads P1-2 and P3-4 incorporate a magnetically hard magnetic moment aligned with their respective paddle orientations. The whole cell structure is repeated along the x-direction with periodic boundary conditions. Solving the non-dimensionalized equation of motion for the system, movement of each bead is determined from spring, bending and fluidic interactions. Additionally, the magnetic moment of beads P1-4 mean they experience a magnetic torque due to an external elliptical magnetic field, which oscillates primarily in the x-direction with a small y-component to break the initial symmetry. Each link is treated as a high-stiffness elastic spring. Bending, described by misalignment of joined links, is resisted by a rotational stiffness that creates a restoring torque between the bead triplets describing the joint. Bead triplets describing misalignment between paddle 1 (paddle 2) and the backbone links have rotational stiffness l 1 (l 2 ). To make the paddles and backbone inflexible, bead triplets defining misalignment within each paddle (P1-B2-P2 and P3-B3-P4) and within the backbone are assigned a much greater rotational stiffness, l paddle = l backbone = 50. Fluidic forces acting on each bead due to drag from the surrounding fluid and flow from other beads are found using expansions that include the Stokes drag and the leading order fluid interaction term. Comparing the system with a pinned swimmer, the net movement of the beads corresponds to a flow generated in the opposite direction 35,36 .

Data availability
The datasets generated and analysed during the current study are available in the University of Exeter Repository, ORE (https://doi.org/10.24378/exe.2923).

Code availability
The MATLAB code used to model the effects of rotational stiffness on flow generated is available in the University of Exeter Repository, ORE (https://doi.org/10.24378/ exe.2923).