## Abstract

Coordinated nonreciprocal dynamics in biological cilia is essential to many living systems, where the emergentmetachronal waves of cilia have been hypothesized to enhance net fluid flows at low Reynolds numbers (*Re*). Experimental investigation of this hypothesis is critical but remains challenging. Here, we report soft miniature devices with both ciliary nonreciprocal motion and metachronal coordination and use them to investigate the quantitative relationship between metachronal coordination and the induced fluid flow. We found that only antiplectic metachronal waves with specific wave vectors could enhance fluid flows compared with the synchronized case. These findings further enable various bioinspired cilia arrays with unique functionalities of pumping and mixing viscous synthetic and biological complex fluids at low *Re*. Our design method and developed soft miniature devices provide unprecedented opportunities for studying ciliary biomechanics and creating cilia-inspired wireless microfluidic pumping, object manipulation and lab- and organ-on-a-chip devices, mobile microrobots, and bioengineering systems.

## INTRODUCTION

Cilia broadly exist in a wide range of organisms on the micrometer scale and can induce notable net fluid flows at very low Reynolds number (*Re*; ~0.001 to 0.01). Their ability of producing net fluid flows plays an essential role in many living systems, such as the self-propulsion of *Paramecium* (*1*); inducing vortical fluid flows for enhancing the feeding function in coral reefs (*2*); and pumping biofluids in the respiratory (*3*), brain ventricle (*4*), and reproductive systems (*5*) in the human body. In addition to the nonreciprocal dynamics of a single cilium, metachronal coordination of cilia arrays is another fascinating behavior of cilia observed in various biological systems (*6*–*9*), which has been hypothesized to be critical to enhance net fluid flows (*10*–*15*). However, experimental validation is still missing, and the fundamental relationship between metachronal coordination and the induced fluid flow has not yet been resolved. This is because biological ciliary systems are not fully viable for controlled experiments, limiting the flexibility of investigating biological hypotheses.

One promising solution is to use small-scale soft devices with programmable motions (*16*–*20*) as a bioinspired platform to model and study their biological counterparts by more controlled experiments at a larger length scale while keeping the essential nondimensional parameters (e.g., *Re*) the same (*21*–*23*). These bioinspired soft devices can also enable unprecedented microfluidic functionalities for lab-on-a-chip (*24*, *25*) and other bioengineering applications (*26*). So far, many small-scale artificial cilia (*27*–*33*) have been proposed for emulating motions produced by biological cilia arrays. However, at small-length scales, designing and controlling both the complex individual ciliary motion and the overall coordinated dynamics within collective artificial cilia arrays pose an enormous conceptual challenge. Metachronal coordination has been mostly investigated in numerical studies (*10*–*15*, *27*). For example, Khaderi *et al*. (*27*) did not experimentally create metachronal waves in artificial cilia, but numerically studied metachronal waves and experimentally demonstrated pumping water by a paramagnetic artificial cilia array with synchronized motions. This may be due to the challenge of controlling locally distributed magnetic fields at small scales. Recent works have proposed several ways to experimentally create metachronal coordination including using the length-dependent buckling motions of paramagnetic sheets (*30*) and magnetic anisotropy in paramagnetic rods (*31*) or ferromagnetic rod arrays (*33*). However, these existing works have focused on the fabrication methods and only demonstrated simple fluid pumping with non-optimized designs. In contrast, we developed design and fabrication methodologies to encode different nonreciprocal kinematics profiles into single-cilium motions and arbitrary phase shifts in a cilia array, even with two-dimensional (2D) metachronal waves and on curved surfaces like their biological counterparts. Moreover, we can optimize both the single-cilium dynamics and cilia array coordination independently to achieve the highest pumping performance by artificial cilia arrays with optimal metachronal coordination. We systematically investigate the quantitative relationship between metachronal waves and their induced fluid flows. These findings also enable us to create cilia-inspired fluidic soft devices with unique capabilities of pumping and mixing viscous synthetic and biological fluids at low *Re*. For example, we can pump even viscous syrup, synthetic mucus, and mouse blood in narrow channels.

## RESULTS

### Design artificial cilia arrays with programmable nonreciprocal motion and metachronal coordination

Here, we present a class of soft miniature devices consisting of multiple submillimeter-scale ferromagnetic-elastic sheets [see fig. S1 (A to C) and the “Fabrication of ferromagnetic cilia,” “Simultaneous encoding of cilia magnetization profiles,” and “Assembly of artificial cilia arrays” sections in Materials and Methods] with identical geometries and nonidentical magnetization profiles. Figure 1 and movie S1 show that these magnetically actuated artificial cilia arrays are capable of creating both programmable nonreciprocal cilia motion and collective metachronal coordination. The complex nonreciprocal oscillating dynamics of a ferromagnetic-elastic sheet (referred as “cilium” later on) in viscous fluids emerges from the nonlinear relations between internal macroscopic elastic stress, hydrodynamic forces, and exerted external magnetic torques at low *Re*, which emulates the regulation of molecular motors in a real cilium. By proper design of the material elastic properties and individual cilium magnetization profiles **M**(*s*) (*s*, material coordinate) as well as the time-varying external actuation magnetic field **B**(*t*) (*t*, time), ferromagnetic-elastic cilia arrays can produce complex but still deterministic ciliary dynamics (Fig. 1A). In particular, with an encoded **M**(*s*) and a rotating **B**(*t*) (see fig. S2 and the “Magnetic actuation of ferromagnetic cilia” section in Materials and Methods for creating such a magnetic field), complex nonlinear oscillations emerge from a rotational buckling instability (*30*) when a cilium undergoes a so-called phase-locking motion (*14*) until it reaches a critical deformed state (*t*/*T*, 0.3 to 1; *T*, beating cycle; Fig. 1B) and then rapidly snaps back in the reverse direction (*t*/*T*, 0 to 0.3; Fig. 1B). This buckling-based motion with two strokes i.e., the power and recovery strokes, results in a nonzero swiping area **r**^{tip}. Specifically, Fig. 1B illustrates a cilium dynamics with a positive swiping area by programming a specific **M**(*s*). Please note that the net fluid flow is independent of the timing of cilium motion at low *Re* but rather depends on the trajectory of an artificial cilium (*8*).

By encoding different **M**(*s*) of a cilium (see fig. S3 and the “Simultaneous encoding of cilia magnetization profiles” section in Materials and Methods for encoding a specific magnetization profile), cilia dynamics with either a positive or negative swiping area can be produced as shown in fig. S4, which leads to a net fluid flow in different directions. The swiping area is a common metric to quantify the nonreciprocal motion and the induced net fluid flow (*11*, *27*, *32*). We define the swiping area ratio (SAR) as a metric to quantify the effectiveness of producing a net fluid flow within a unit time by a single cilium of length *L*, where *Re*. A computer-aided method based on a fluid-structure interaction (FSI) model (see the “Summary of the used coordinate systems” and “FSI model of single cilium dynamics” sections in Materials and Methods) was developed to optimize SAR. The single-cilium dynamics mainly depends on the material elastic modulus, **M**(*s*), and **B**(*t*). We specifically optimized **M**(*s*) and **B**(*t*) to maximize SAR, as shown in fig. S5 and described in the “Optimizing a single cilium by maximizing the SAR” section in Materials and Methods. The notable advantage of our design method and system is that the individual cilium beating dynamics can be designed separately from their metachronal waves, giving substantial flexibility and viability of the collective motion.

Figure 1 (C and D) shows that arbitrary 1D metachronal waves can be further produced within an array of magnetic artificial cilia, by encoding phase shifts Δϕ* _{x}* = ϕ

_{i + 1}(

*s*) − ϕ

*(*

_{i}*s*) in

**M**(

*s*) (Fig. 1E) and applying a rotating

**B**(

*t*). The dynamics of each cilium remains the same for one period (see the “Mechanics and design methodology of metachronal waves” section in Materials and Methods for the mechanics). These 1D metachronal waves can be quantified by the lags in the motions of neighboring cilia, such as the lag time Δ

*t*in the periodic beating dynamics when each cilium tip reaches the maximal

*y*position as shown in Fig. 1D. The phase shift in motions can be defined as

*t*∈ [0,

*T*)). The

*x*component of the wave vector is thus given by

*k*= Δψ

_{x}*/*

_{x}*d*for quantifying the spatiotemporal coordination. Consequently, Fig. 1F shows that Δϕ

_{c}*can be linearly mapped to Δψ*

_{x}*and*

_{x}*k*, which suggests that arbitrary 1D metachronal waves can be encoded into neighboring cilia. Figure 1F also shows the so-called symplectic (Δψ

_{x}*> 0) and antiplectic (Δψ*

_{x}*< 0) metachronal waves (*

_{x}*7*), where the direction of the wave propagation and the power stroke are in the same and opposite directions, respectively.

Similar to biological cilia, Figure 1G and movie S2 show that our method also allows designing and fabricating artificial cilia arrays to emulate the propagation of 2D metachronal waves in both the *x-y* and *x-z* planes with arbitrary 2D wave vectors, by encoding both nonzero Δϕ* _{x}* and Δϕ

*into a 2D array. Different from the emergence of metachronal waves due to hydrodynamic interaction (*

_{z}*12*) and distal coupling (

*34*) in biological systems, the metachronal waves in our system emerge from the rotational shift in the magnetization phase profile of each magnetic cilium subject to a rotating uniform magnetic field. Nonetheless, the resulting metachronal waves resemble those of their biological counterparts and can be used for investigating the function of metachronal coordination in producing fluid flows at low

*Re*.

### Optimal metachronal waves for pumping fluid flows efficiently

To find the optimal metachronal waves for pumping directional fluid flows, Fig. 2 (A to C) compares the directional fluid flow produced by our artificial magnetic cilia arrays with different *k _{x}* (or Δψ

*) (Fig. 2B) and intercilium spacing*

_{x}*d*(Fig. 2C), based on the particle image velocimetry (PIV) measurements inside glycerol (dynamic viscosity, 0.876 Pa·s; see the “Fluid flow measured by PIV” section in Materials and Methods for more details). Because of a small number of cilia at low

_{c}*Re*, fluidic boundary/wall effects (

*11*,

*35*) could be substantial if the cilia are placed inside confined fluidic channels or other confined spaces. Here, we avoided such effects by using nearly semi-indefinite boundary configuration as shown in Fig. 2A and fig. S2F. To quantify the performance of pumping directional fluid flows, we use the widely used average volume flow rate

*11*–

*13*) in the direction of transportation (+

*x*direction). In Fig. 2B,

*k*is varied from −2π/(3

_{x}*L*) to 5π/(9

*L*) when Δψ

*is varied from –π to 5π/6, while*

_{x}*d*equals to 1.5

_{c}*L*and other parameters are kept identical. It is found that antiplectic waves produce a maximal directional fluid flow when Δψ

*∈ [−π/2, −π/6]. For example, artificial cilia arrays with Δψ*

_{x}*= −π/3 achieve a*

_{x}^{3}·s

^{−1}, which is about 1.6 times of that produced by arrays with synchronized motions and 4 times of the minimal net fluid flow produced by arrays with a symplectic wave (Δψ

*= π/3). In addition, Fig. 2C shows that artificial cilia arrays with denser spacing yield a larger*

_{x}*d*from

_{c}*L*to 2

*L*while keeping Δψ

*= − π/3 and other parameters identical. The quantitative relationship between*

_{x}*k*and

_{x}*L*) and −π/(9

*L*), which yields the maximal

*k*from π/(6

_{x}*L*) to π/(3

*L*), by reducing

*d*from 2

_{c}*L*to

*L*, an increase of

*and intercilium spacing*

_{x}*d*on the observed fluid flow. Our experiments prove the key role of antiplectic metachronal waves in enhancing fluid transportation. The finding is in accordance with that in biological systems (

_{c}*7*,

*36*) and numerical works (

*12*,

*13*) and introduces extra insights into the special function of antiplectic waves by reporting the specific range of wave factors for producing optimal net fluid flows.

While *x* direction (Fig. 2D). Tracer particles (300 to 330 μm in diameter; density, 1.2 g·cm^{−3}) are transported in the +*x* direction by artificial cilia arrays with *k _{x}* varying from −2π/(3

*L*) to 5π/(9

*L*) when Δψ

*∈ [− π,5π/6] and*

_{x}*d*equaling to 1.5

_{c}*L*(see the “Particle transportation using cilia arrays” section in Materials and Methods for experimental details). The Lagrangian coherent structures of the local time-periodic fluid flow guarantee consistent trajectories of neutrally buoyant tracer particles (

*37*). Figure 2E shows that the artificial cilia arrays with

*k*∈ [−π/(3

_{x}*L*), −π/(9

*L*)] when Δψ

*∈ [−π/2, −π/6] could pump particles with a maximal average speed (~0.45 mm/s), agreeing with the trend in Fig. 2B. At the same time, artificial cilia arrays with denser spacing result in larger particle transportation speeds as shown in Fig. 2F, consistent with the results in Fig. 2C.*

_{x}Comparing the particle transportation trajectories shows that another function of the optimal antiplectic waves is attracting free particles. Figure 2G shows that metachronal waves with Δψ* _{x}* varying from

*d*(

_{p}*x*). Such a result could support the finding that biological cilia provide a key functionality of predation by attracting specific microbiomes toward their surfaces (

*38*).

### Mechanism of antiplectic waves enhancing fluid flows

The reason of the superior ability of enhancing fluid flows by antiplectic waves is further investigated. Figure 3 shows that metachronal waves can modulate local fluid flows by producing different time-varying boundary conditions to the fluid, thereby varying the overall fluid transportation performance. This observation is obtained by the PIV data of artificial cilia arrays with representative metachronal waves. The PIV results in Fig. 3A show that the local fluid flow around cilium no. 3 is strongly influenced by its neighboring cilia nos. 2 and 4. Compared with an artificial cilia array showing a symplectic wave, such as *t*_{1} = 0.12*T*; Fig. 3A). This leads to an increased and more continuous local flow in the power stroke direction. Meanwhile, the flow is more blocked from neighboring cilia during the recovery stroke (*t*_{2} = 0.45*T*; Fig. 3A), resulting in a reduced and blocked local fluid flow in the recovery stroke direction.

To quantitatively see such a difference in local flows, Fig. 3B compares the instant volume flow rate *Q _{x}* passing through a selected cross section within a full period. The observation is that although an artificial cilia array with synchronized motion induces a large and single peak value of

*Q*due to the superposition effect of all cilia, artificial cilia arrays with metachronal coordination yield smaller but multiple peak values of

_{x}*Q*. Therefore, to make a fair comparison, we compare

_{x}*Q*produced by two artificial cilia arrays with opposite phase shifts (

_{x}*Q*during the power stroke and a smaller value of

_{x}*Q*during the recovery stroke, resulting in a larger net fluid flow within a period.

_{x}Simple metrics are further proposed to quantify the boundary conditions and predict the net fluid flow by artificial cilia arrays with different metachronal waves. To compare different instant boundary conditions, in Fig. 3C, we choose the instant inter-tip spacing *i* = 2 to 5) between a selected cilium and its two neighbors to measure the no-slip boundary conditions. Figure 3D shows that antiplectic waves yield a larger *d*_{power} and *d*_{recovery}, respectively, which quantify the average boundary conditions produced by different cases of metachronal coordination. In Fig. 3E, experimental data show that the overall trend of the quantitative relationship between * _{x}* presented in Fig. 2B could be captured by a proportional law given by

### Functional artificial cilia as fluidic devices by encoding optimal coordination

The established quantitative relationship and mechanism between fluid flows and metachronal coordination further allow creating bioinspired cilia arrays as versatile fluidic devices with unique and important functionalities at low *Re*. We demonstrate multiple functional fluidic devices by encoding diverse metachronal coordination and nonreciprocal motions.

In Fig. 4 and movie S4, we first demonstrate a bioinspired cilia array with optimal metachronal waves propagating on curved surfaces inspired by coral reefs (*2*). Existing artificial cilia are only capable of pumping on flat surfaces (*27*–*29*, *32*, *33*), which does not resemble their biological counterparts in terms of boundary morphologies (*7*). In contrast, our developed artificial cilia array is capable of transporting particles along arbitrary curved surfaces, such as an S-shaped surface, as shown in Fig. 4A. To create bioinspired cilia planted on a boundary wall with much more complex curvatures, we generalize our method to include the cilium-attached coordinates {*L*_{si}} and the individual dynamics of each cilium as two extra sets of design variables (see the “Encoding metachronal waves in cilia arrays on curved surfaces” section in Materials and Methods). The key design rule is to compensate the rotation of the cilium-attached coordinates when encoding the magnetization profile of each cilium such that a desired constant phase shift is kept between neighboring cilia even on curved surfaces (fig. S7).

With this design methodology, we develop a bioinspired cilia array with optimal metachronal waves, which can pump viscous fluids along curved surfaces. Figure 4B shows that such a bioinspired cilia array could transport particles tangentially to the attaching surface more efficiently, compared with an array with the same boundary wall but with synchronized motions due to the locally enhanced vortical flow, as visualized by dye and PIV in Fig. 4 (C and D, respectively). More bioinspired cilia arrays could be constructed with various complex boundary morphologies, such as large-polyped, mounding, and encrusting colony morphologies (*39*). Important scientific questions, such as how geometrical details of colony morphology affect the transportation performance of coral reefs (*2*) and similar cilia-covered organisms, can be potentially answered in the future to help researchers better understand the active mass transportation in ambient aquatic environments.

Pumping viscous fluids in narrow channels at low *Re* is challenging because of the no-slip boundary conditions (*40*) but is important in many biological systems, such as tubal transportation of ova in female fallopian tubes, where ciliary motion contributes to efficient transportation in addition to muscular contractility (*5*). In Fig. 5 and movie S5, we demonstrate efficient transportation of particles in narrow channels filled with glycerol, by mimicking the natural cilia in tubal transportation. As shown in Fig. 5A, we jointly designed two cilia arrays with opposite SARs (see fig. S4), as their pumping directions are in the same direction (in {**L*** _{a}*}) subject to the same rotating

**B**(

*t*). Within each array, an optimal antiplectic metachronal wave is further encoded. In this way, the device transports fluids and particles faster (

*41*).

In addition to the functionality of transportation, in Fig. 6 and movie S6, we demonstrate that optimal metachronal coordination can be encoded into two arrays of artificial cilia for efficient fluid mixing function. Mixing adjacent laminar flow completely and rapidly at low *Re* is difficult while important for many applications (*25*). At low *Re* (<0.1), passive mixing devices usually require specific design of channel geometries, using chemical reactions, limiting the applicable scenarios. Other active mixing devices, such as magnetic micropillars, have limited mixing thoroughness between vertical layers due to their cone-shaped stirring motions (*28*). Here, we demonstrate a mixing device capable of mixing viscous fluids thoroughly and rapidly by encoding optimal coordination into two artificial cilia arrays. The cilia in the bottom array in Fig. 6A have an optimal metachronal wave (*25*) less than 0.05 within 35 s in a channel across a volume of 43 mm^{3} at a local *Re* < 0.03, as shown in Fig. 6 (B and C).

Last, the wireless magnetic actuation of our fluidic devices combined with ultrasound medical imaging could allow remote pumping and imaging of body fluids in enclosed small spaces with no physical contact, potentially applicable in the human body in the future for biomedical applications. As a proof of concept, in Fig. 7 and movie S7, we demonstrate the ability of pumping body fluids, such as fresh whole mouse blood in narrow channels (see the “Preparation and characterization of fluids” section in Materials and Methods), showing the promising function of pumping biofluids ex vivo and in vivo in future biomedical applications. They could potentially be deployed by a medical catheter or locomote by remote actuation (*16*) to target sites in the human body for delicate and minimally invasive noncontact fluidic or object/cell manipulation (*42*). These devices have potential to be applied in health care for patients with specific ciliary dysfunctions, such as assisting in pumping mucus in human respiratory systems and transporting ovum in female fallopian tubes for improving fertilization (*5*).

## DISCUSSION

We have developed bioinspired ciliary devices, which not only provide remarkable insights of how the time-varying boundaries of ciliary surfaces by different metachronal coordination could lead to different cooperatively generated ciliary flows but also use these findings to enable unprecedented optimized engineering applications. The reported design methodology is applicable in principle to scale the artificial cilia down to the size of real cilia (~10 μm in length). A scaling analysis (see the “Scaling analysis of single cilium dynamics” section in Materials and Methods) shows the feasibility of scaling down the overall size of ferromagnetic sheets to the micrometer scale. Microfabrication of ferromagnetic materials using recent advance of fabrication methods (*43*) could potentially allow manufacturing smaller artificial cilia arrays with metachronal waves in the future, although challenges of weaker magnetization, limited geometry, and limited material properties still need to be overcome.

Biological cilia are densely packed rod structures capable of nonreciprocal 2D or 3D motions (*44*), with a length scale of ~10 μm, beating at 5 to 30 Hz (*7*). In our system, instead of mimicking the biological cilia, we aim at developing artificial cilia with metachronal coordination similar to their biological counterparts, by keeping the critical dimensionless number (e.g., *Re*) similar. As the fundamental physical law is the same, we can study a scientific question in cilia mechanics: the function of the metachronal coordination in inducing fluid flows at low *Re*. With such a platform, we find that antiplectic waves can indeed enhance fluid flows within the scope of our investigation space. This finding agrees with that in biological systems (*7*, *36*) and numerical studies reported before (*12*, *13*), bringing extra insights into the ciliary biological systems.

In addition, the developed artificial cilia array with programmable nonreciprocal motion and metachronal coordination could be used as a scientific tool to further investigate many interesting scientific questions, such as how the fluid flow patterns would be different when varying the cilia beating dynamics (*45*), boundary geometries (*38*), location of topological defects (*6*), and the conditions of the surrounding fluidic environment (*46*). We anticipate our design method and developed soft miniature devices that could open a wide range of unprecedented opportunities for studying ciliary biomechanics and creating cilia-inspired wireless microfluidic pumping, object manipulation and lab- and organ-on-a-chip devices, mobile microrobots, and bioengineering systems.

## MATERIALS AND METHODS

### Fabrication of ferromagnetic cilia

The ferromagnetic-elastic sheet–based artificial cilia were made of composite materials of silicone rubber (Ecoflex 00-30, Smooth-On Inc.) and NdFeB hard magnetic microparticles (average diameter, 5 μm; MQFP-15-7, Neo Magnequench) with a weight ratio of 1:1. The two materials were mixed together and then poured on to a poly(methyl methacrylate) substrate with tapes as spacers for controlling the thickness to be 100 μm. The top surface of the mixed material was scraped by a single-edge razor blade (fig. S1A). The composite material was cured in room temperature for about 4 hours or cured in 1 hour by putting it on a hot plate at 60°C. After curing, as illustrated in fig. S1B, the ferromagnetic cilia were cut by a laser machine (LPKF ProtoLaser R3, LPKF Laser & Electronics AG) according to a designed dimension of 1 mm by 550 μm by 100 μm (*L* × *w* × *t _{b}*) marked in fig. S1C. Other dimension parameters in fig. S1C are given by

*L*

_{1}= 300 μm,

*w*

_{1}= 850 μm, and

*w*

_{2}= 200 μm. The material has a density of (2.00 ± 0.056) × 10

^{3}kg·m

^{−3}and an average Young’s modulus of 144 kPa measured by a tensile testing machine (5940 series, Instron GmbH). In the experiments, we observed that the performance of the device is quite robust without any degradation issue. The dynamics of the artificial cilia array do not change even several months after its first use, when actuated with the same control signals in the same fluids. According to a previous study, structures made of Ecoflex 00-30 can maintain their mechanical properties after being subjected to cyclic loads for over 10

^{6}cycles (

*47*).

### Simultaneous encoding of cilia magnetization profiles

We used a jig-assisted method reported in our previous work (*18*) but extended it to simultaneously encode the desired different **M**(*s*) of multiple ferromagnetic cilia using one magnetizing jig. Figure S1D illustrates the principle and process to encode the desired different **M**(*s*) on multiple cilia simultaneously. To obtain a desired magnetization profile, multiple ferromagnetic cilia are prebent into a jig with multiple cutout parts in designed shapes. They are then magnetized by applying a large magnetizing **B** (magnitude, 1.2 T) in the +*x* direction inside a vibrating sample magnetometer (EZ7 VSM, MicroSense LLC). The shape of the cutout part corresponding to each cilium is parameterized with a local slope angle profile * _{i}*(

*s*) = ϕ

^{0}(

*s*) + (

*i*− 1)Δϕ for a given cilium, we have

**M**

_{i}(

*s*) in an array simultaneously, the inverse design of the cutout part of the magnetizing jig is thus given by

Figure S3 shows that we encoded the desired **M**(*s*) by comparing the measured and predicted magnetic fields produced by a magnetized single cilium and array of cilia at their surfaces.

### Assembly of artificial cilia arrays

Figure S1 (E and F) illustrates the procedure of assembling ferromagnetic cilia into an array. The magnetized cilia were manually attached to an assembly fixture and fixed by glue (Loctite 401, Henkel AG) under a stereomicroscope (Stemi 508, Carl Zeiss AG). Motorized micromanipulators could be potentially incorporated for automating such a task in the future.

### Magnetic actuation of ferromagnetic cilia

As shown in fig. S2 (A and B), a rotating Halbach array composed of 12 cubic ferromagnets (12 mm by 12 mm by 12 mm; N42, Supermagnet.de) was used for generating a rotating uniform magnetic field. The rotational motion was actuated by a dc motor (Parallax Continuous Rotation Servo, Parallax Inc.) and controlled by an embedded controller (Arduino Uno, Arduino.cc). **B**(*t*) was characterized in fig. S2 (C to E). The magnetic field strength was about 40 mT (at 15 mm above the surface of the magnet array). On this plane, the uniformity of the magnetic field in the workspace was ~95% within a circular area of 10 mm in diameter, which covered the whole artificial cilia array in the experiments. The artificial cilia array was fixed in a container (43 mm by 32 mm by 11 mm; fig. S2F) filled with glycerol (99.8% volume ratio) and supported on a customized sample holder above the Halbach array. We used glycerol as a viscous liquid in the experiments because it has a high viscosity to create a low *Re* environment and does not cause swelling of the elastic cilia. The fluid had a density of 1.257 ×10^{3} kg·m^{−3} and a dynamic viscosity of 0.876 Pa·s at the room temperature 25°C.

### General design methodology

*Summary of the used coordinate systems*. Three coordinate systems in this work are presented in fig. S2 (A and G) and summarized as below. First, the global coordinate system {*G*_{b}} (*x*_{b}, *y*_{b}, *z*_{b}) is located at the center of the Halbach array for expressing the global actuation magnetic field and location of the artificial cilia arrays. Second, the array-attached coordinate system {*L*_{a}} (** x**,

**,**

*y***) is located on the base of the artificial cilia array, which is used to express the dynamics of artificial cilia, fluid flow, and particle transportation. Last, the cilium-attached coordinate system {**

*z*

*L*_{si}} (

**,**

*X***,**

*Y***) of the**

*Z**i*-th cilium in an array is located on the body of the cilium, which is used to describe its relative orientation and position within an array and to express the magnetization profile of each cilium.

*FSI model of single-cilium dynamics*. An FSI model in 2D is presented to guide the design of single ferromagnetic-elastic cilia. The motion of a given cilium is modeled using the Euler-Bernoulli beam theory for large deflections (*18*) while considering fluid drags. At quasi-static states, the governing equations are given below for a cilium with one end fixed to the boundary wall. For an infinitesimal element [*s s* + d*s*] of the *i*-th cilium (*i* = 1,2, …, *n*) in an array, the moment balance equation at a quasi-static state in the Eulerian space is given by* _{i}*(

*s*,

*t*) is the rotation angle at time

*t*,

*E*is the Young’s modulus,

*A*

_{cross}=

*wt*is the cross-sectional area,

_{b}**n**

*is the unit vector along the +*

_{z}*Z*axis,

**B**(

*t*) =

*B*[cos (ω

_{m}*t*) sin (ω

*t*) 0 ]

*is the rotational external magnetic field in the*

^{T}*X*-

*Y*plane,

**M**

*(*

_{i}*s*) = [

*M*(

_{ix}*s*)

*M*(

_{iy}*s*) 0]

*is the magnetization profile function, and*

^{T}**F**(

*s*,

*t*) = [

*F*]

_{h}F_{v}*is the internal force of the cilium.*

^{T}The force balance equation is given by*A*_{normal} = *wds* is the side area of the infinitesimal element and **F*** _{d}* is the fluid drag force per volume on the infinitesimal element.

**F**

*is a function of the local drag coefficient*

_{d}*C*assuming a cylinder geometry based on local curvatures (

_{d}*48*) and the relative velocity between the element and fluid flow.

**F**

*is further given by*

_{d}*is the density of the fluid and β is a scaling factor to calibrate the damping coefficient, as it is generally difficult to exactly match the local drag coefficient*

_{f}*C*(

_{d}*48*,

*49*). The local

*Re*is given by

*is the dynamic viscosity of fluids. In addition, the Cartesian position*

_{f}**u**

*(*

_{s}*s*,

*t*) of the cilium is given by

The boundary conditions are given by a set of equations

The initial conditions are given by another set of equations

The cilia motions are assumed to be planar by imposing a large width-to-length ratio (*w*/*L* = ~0.6) to minimize the out-of-plane twisting and bending motions. The single-cilium oscillating beating dynamics under a rotating **B**(*t*) can be predicted using the above model as shown in fig. S8.

*Optimizing a single cilium by maximizing the SAR*. The SAR is optimized on the basis of the discussed FSI model in the above section. The dynamics of a single cilium in viscous fluids is a function of **M**(*s*), its physical dimensions, and **B**(*t*). The 2D magnetization profile along the *i*-th cilium with a length *L*, is given by

We first fix the magnitude of **M**(*s*), as 40 kA·m^{−1} and assume that the magnitude of external magnetic field is *B _{m}* = 40 mT. Then, for simplicity, ϕ

*(*

_{i}*s*) is represented in a general continuous polynomial function with a first-order approximation as

*(0) − ϕ*

_{i}*(*

_{i}*L*) of a single cilium. We found that a maximal positive SAR could be obtained when ϕ

*(0) − ϕ*

_{i}*(*

_{i}*L*) = 1.75π (used in all experiments if not explicitly mentioned) and a maximal negative SAR could be obtained when ϕ

*(0) − ϕ*

_{i}*(*

_{i}*L*) = 0 (used in Fig. 5). These two magnetization profiles are optimal within a practical range (−1.75π to 1.75π) constrained by the template-based magnetizing method.

The magnitude *B _{m}* and frequency

*f*of

**B**(

*t*) are another two control inputs of the dynamics of a single cilium. Increasing

*f*will increase SAR, but the maximal

*Re*will also increase (fig. S5C). Therefore, we constrained

*f*to ensure that our cilia arrays operate at the low

*Re*regime in the experiments. In addition, a larger

*B*yields a larger SAR (fig. S5D) until the cilium’s deformation is limited by the boundary walls.

_{m}In addition, the dimensions of the artificial cilium also affect the motion. We defined two ratios, the thickness-to-length ratio (*r _{tl}*) and the width-to-length ratio (

*r*), to explain their effect based on our models.

_{wl}*r*has a major effect on the magnitude of the cilia motion. With the same actuation magnetic fields, artificial cilia with a larger

_{tl}*r*will exhibit a smaller maximal bending angle and a smaller SAR, as shown in fig. S9, while the buckling motion is more evident, yielding a larger output force during the power stroke after being deformed to the same state. In addition,

_{tl}*r*has less effect on the planar cilia motion according to Eq. 3 in the “FSI model of single cilium dynamics” section, while mainly affecting the out-of-plane motion of a single cilium. By increasing the width to length ratio (~0.6), the undesired out-of-plane motions could be minimized.

_{wl}*Mechanics and design methodology of metachronal waves*. Each member within an array of ferromagnetic cilia is a distributed oscillator subject to a rotating magnetic field **B**(*t*) = R* _{z}*(ω

*t*)

**B**(0). The phases of these oscillators can be represented by their body slope angles θ

*(*

_{i}*s*,

*t*). For two neighboring cilia, we have the following moment-balancing equations

The temporal differences in the time-varying and spatially distributed external magnetic torques *i*-th and (*i* + 1)-th cilia are given by

The phase shift Δϕ in **M**(*s*) can be equivalently represented as Δϕ in **B**(*t*) due to the mathematical properties of the cross product in 2D as shown in Eq. 15. Therefore, the phase shifts

*CFD simulations in comparison with experimental data*. To model the fluid flow created by the artificial cilia arrays, we use a 3D bidirectional CFD model (*15*), which can accurately describe the magnetic torque on the cilia, the deformation of the cilia, and the deformation-induced fluid flow at low *Re* [see fig. S6 (A and B)]. The model accounts for the fully coupled bidirectional solid-fluid interaction between cilia and fluid, including the cilia deformation induced by the fluid. The computational framework is briefly summarized here. In the CFD model, the cilia are represented by an assemblage of shell elements, which act as an internal boundary to the fluid. The FSI is considered by implicitly coupling the fluid dynamics and solid mechanics equations, where the Stokeslet method is used to account for the viscous drag and implemented using a boundary-element approach. To model the ciliary deformation, we use shell elements with the bending and membrane behavior accounted for by using triangular Kirchhoff elements (*50*) and constant strain triangles with drilling degrees of freedom (*51*), respectively. The large and geometrically nonlinear deformation of the cilia is modeled by adopting an updated Lagrangian approach. We focus on flows at low *Re* and use Green’s functions (*52*) to simulate fluid motion in the Stokes regime. The solid and fluid are implicitly coupled through a no-slip boundary condition.

We compared the effect of the phase shift Δψ* _{x}* (fig. S6C) and intercilium spacing

*d*(fig. S6D) on inducing fluid flow in the simulated and experimental scenarios. The simulations also show that the antiplectic waves produce a maximal directional fluid flow when

_{c}*Encoding metachronal waves in cilia arrays on curved surfaces*. The metachronal coordination can be encoded into cilia arrays planted on a curved boundary wall using the design rule of Δβ = Δϕ + Δα. Here, Δα = α_{i + 1} − α* _{i}* describes the variation of the slope angles along the curved surface by defining the rotational difference between the cilium-attached coordinates {

*L*_{s(i+1)}} and {

*L*_{si}} of two neighboring cilia. The distributed magnetic torques induced on these cilia have a relationship of

To create metachronal waves on arbitrary curved surfaces, ϕ* _{i}*(

*s*) of each cilium is inversely designed together with {

*L*_{si}} as

*X**axis of {*

_{i}

*L*_{si}} relative to the surface tangent vector direction along the boundary wall surface. Specifically, in the bioinspired cilia array in Fig. 4,

*x*(

_{ci}*i*= 1,2, …,8) was chosen with an equal spacing from −4

*L*to 4

*L*.

*Scaling analysis of single-cilium dynamics*. We first define *VE* to quantify the relative scaling of viscous drag and elastic stress on an infinitesimal element, which is given by*L*^{−2}. This equation shows that the viscosity of the fluid needs to be decreased for reducing the damping effect, *L* scales down.

We then define *ME* to quantify the relative scaling of magnetic torque and elastic stress on an infinitesimal element, which is given by

This equation shows that the same **M**(*s*) and **B**(*t*) can induce similar elastic stress within the cilium when *L* scales down.

One practical limitation of our system is that the intercilium spacing *d _{c}* cannot be less than a body length. As

*d*further decreases, the dynamics of artificial cilia may be affected by magnetic local interactions. The critical distance between neighboring cilia is a complex function of the magnetization profile

_{c}**M**(

*s*), the material stiffness, and the magnitude of external magnetic fields

**B**(

*t*). To alleviate such an effect, it is found that we could reduce the magnitude of

**M**(

*s*), increase the material Young’s modulus

*E*, or increase the magnitude of

**B**(

*t*). More research on this critical distance will be carried out in the future.

### Fluid flow measured by PIV

A PIV system (Dantec Dynamics, Dantec Inc.) was integrated with the magnetic actuation setup for visualizing and quantifying the induced fluid flows (fig. S1A). The fluorescent tracing particles (Cospheric Inc.) used in the PIV experiments had an average diameter of 2 μm. A high-speed camera (Phantom MicroLab 140, Vision Research, Ametek Inc.) with a framerate up to 5000 frames per second (fps) was used for recording the dynamics of cilia and the motion of tracing particles. The acquired image sequences were further analyzed using the software DynamicStudio 6.1 (Dantec Dynamics, Dantec Inc.) and customized codes implemented in MATLAB 2017a (MathWorks Inc.). As illustrated in Fig. 2A, we obtained the average volume flow rate *x* direction at *x* = *x*_{0}. Three layers of identically designed arrays are stacked in the *z* direction to emulate a ciliary carpet that has a larger fluid flow than a single layer. The image plane is selected as the middle plane of three layers (see fig. S1B and Materials and Methods for more details). The key parameters for calculating *y*_{0} = 0.5*L*, Δ*y* = 3*L*, Δ*z* = 3*w*, ϵ = 0.1*L*, and *m* = 4. All variables were represented in the array-attached coordinate system {*L*_{a}}.

### Particle transportation using cilia arrays

Fluorescent polyethylene microspheres (1.20 g·cm^{−3}, 300 to 355 μm in diameter) were used as neutrally buoyant tracer particles because their densities are close to glycerol (1.257 g·cm^{−3}), so that the particle’s motion reflects the fluid velocity by minimizing the effect from the gravity and buoyancy. Two cameras (Flea 3 USB3, FLIR Systems Inc.) were used for recording the trajectories of particles at about 30 fps. Particles trajectories were extracted by segmenting particle contours from the video sequences using color-based filtering algorithms implemented in MATLAB R2017a (MathWorks Inc.). The average speeds of particles were calculated along a trajectory with a displacement of 7 mm in the + *x* direction (Fig. 2G). The SDs were for *n* = 5 trials using each artificial cilia array. For each trial of the experiment, the particle was released from the same position.

### Preparation and characterization of fluids

The fluids used for pumping demonstrations included deionized water, mouse blood, synthetic mucus, and syrup. Fluids of low viscosity include fresh mouse blood (Innovative Grade US Origin Mouse CD 1 Whole Blood; Fig. 7A) and deionized water (Fig. 7B). The mouse blood is shear-thinning with a dynamic viscosity of ~10 mPa·s at a shear rate of 2.4 s^{−1} (*53*). Viscous fluids in addition to glycerol include synthetic mucus (synthetic nasal mucus, Kryolan Professional Make-up GmbH; Fig. 5C) and syrup (rice syrup, Serapis Culinary; Fig. 5D). Figure S10 shows the rheology properties of the synthetic mucus and syrup. The rheology test is carried out in a Discovery HR-2 rheometer from TA Instruments Inc.

### Statistical analysis

The SDs are indicated in the figure captions.

## Supplementary Materials

Supplementary material for this article is available at http://advances.sciencemag.org/cgi/content/full/6/45/eabc9323/DC1

This is an open-access article distributed under the terms of the Creative Commons Attribution-NonCommercial license, which permits use, distribution, and reproduction in any medium, so long as the resultant use is **not** for commercial advantage and provided the original work is properly cited.

## REFERENCES AND NOTES

**Acknowledgments:**We thank all members of the Physical Intelligence Department at the Max Planck Institute for Intelligent Systems for their comments.

**Funding:**This work is funded by the Max Planck Society and European Research Council (ERC) Advanced Grant SoMMoR project with grant no. 834531.

**Author contributions:**X.D., G.Z.L., W.H., and M.S. proposed the research. X.D. designed the experiments. X.D., W.H., and Z.R. performed the experiments. X.D. analyzed the data. R.Z. and P.R.O. performed the CFD simulations. M.S. planned and supervised the research. X.D., W.H., and M.S. wrote the manuscript with input from all authors. All authors commented on or edited the manuscript.

**Competing interests:**The authors declare that they have no competing interests.

**Data and material availability:**The modeling and design of the single-cilium dynamics was performed using customized codes implemented in MATLAB R2017a (MathWorks Inc.). The CFD analysis was performed using customized codes using FORTRAN (IFORT 13.1.2). These codes are available upon reasonable request to the corresponding author. All data needed to evaluate the conclusions in the paper are present in the paper and/or the Supplementary Materials. Additional data related to this paper may be requested from the authors.

- Copyright © 2020 The Authors, some rights reserved; exclusive licensee American Association for the Advancement of Science. No claim to original U.S. Government Works. Distributed under a Creative Commons Attribution NonCommercial License 4.0 (CC BY-NC).