1. Introduction
Cilia are thin hair-like cellular protrusions that serve a variety of fundamental roles in many eukaryotes. The internal structure has a characteristic ‘9+2’ arrangement of microtubules. Driven by the distributed sliding forces acting on neighbouring microtubules by molecular motors, cilia beat asymmetrically with a distinct power and recovery stroke to break the time-reversal symmetry and generate net propulsion at low Reynolds number (
$Re$) (Blake & Sleigh Reference Blake and Sleigh1974; Purcell Reference Purcell1977). One hypothesis for the spontaneous beating is based on dynamic instabilities developed when the motor activity exceeds a threshold (Camalet, Julicher & Prost Reference Camalet, Jülicher and Prost1999; Bayly & Dutcher Reference Bayly and Dutcher2016; Oriola, Gadêlha & Casademunt Reference Oriola, Gadêlha and Casademunt2017; Ling, Guo & Kanso Reference Ling, Guo and Kanso2018). Cilia can also beat collectively in dense arrays to form metachronal waves via hydrodynamic interactions (Meng et al. Reference Meng, Bennett, Uchida and Golestanian2021; Chakrabarti, Fürthauer & Shelley Reference Chakrabarti, Fürthauer and Shelley2022; Kanale et al. Reference Kanale, Ling, Guo, Fürthauer and Kanso2022), and such collective motions play crucial roles in the locomotion and material transport of many microorganisms and tissues (Lauga & Powers Reference Lauga and Powers2009; Faubel et al. Reference Faubel, Westendorf, Bodenschatz and Eichele2016; Juan et al. Reference Juan, Guillermina, Mathijssen, He, Jan, Marshall and Prakash2020).
Inspired by natural cilia, the design and fabrication of artificial cilia have attracted growing interest and are important for a wide range of applications, such as propelling microrobotics (Ye, Régnier & Sitti Reference Ye, Régnier and Sitti2013; Lum et al. Reference Lum, Ye, Dong, Marvi, Erin, Hu and Sitti2016; Liu et al. Reference Liu, Qin, Zhu, Yang and Luo2020; Hu, Zhang & Shelley Reference Hu, Zhang and Shelley2022), and pumping and mixing fluid in fluidics. In artificial systems at microscales, the beating mechanism of natural cilia seems impractical to realize. To mimic the ciliary beating patterns and generate non-reciprocal motions, various external actuation mechanisms have been explored, including light (Van Oosten, Bastiaansen & Broer Reference Van Oosten, Bastiaansen and Broer2009), pneumatic (Milana et al. Reference Milana, Gorissen, Peerlinck, De Volder and Reynaerts2019), electric fields (den Toonder et al. Reference den Toonder2008) and especially magnetic fields (Khaderi et al. Reference Khaderi, Baltussen, Anderson, Ioan, Den Toonder and Onck2009; Shields et al. Reference Shields, Fiser, Evans, Falvo, Washburn and Superfine2010; Khaderi et al. Reference Khaderi, Craus, Hussong, Schorr, Belardi, Westerweel, Prucker, Rühe, Den Toonder and Onck2011; Lum et al. Reference Lum, Ye, Dong, Marvi, Erin, Hu and Sitti2016; Hanasoge et al. Reference Hanasoge, Ballard, Hesketh and Alexeev2017; Meng et al. Reference Meng, Matsunaga, Yeomans and Golestanian2019; Dong et al. Reference Dong, Lum, Hu, Zhang, Ren, Onck and Sitti2020; Gu et al. Reference Gu2020).
Instead of applying distributed motor forces, De Canio, Lauga & Goldstein (Reference De Canio, Lauga and Goldstein2017) showed that a single tangential follower force acting on the tip of a clamped elastic filament in a viscous fluid can also induce buckling and spontaneous oscillations through a Hopf bifurcation. Using the same phenomenological model, Man & Kanso (Reference Man and Kanso2020) demonstrated that multiple active filaments can display different synchronization states through hydrodynamic interactions. Although this driving mechanism does not require an external time periodicity, it is considered theoretically and not as a practical fluid pumping mechanism since the follower force is required to always remain tangential to the filament tip. To efficiently pump fluid at low 
$Re$, non-reciprocal trajectories with large swept areas are needed (Osterman & Vilfan Reference Osterman and Vilfan2011), which may be achieved through the buckling of an elastic filament under compression.
In this work, we propose a fluid pumping mechanism that can both exploit the buckling instability and is also more practical experimentally. We consider a composite cilium consisting of an elastic filament and a spherical particle attached at the filament tip, moving in a three-dimensional Stokesian fluid. The filament base is clamped to a no-slip wall and tilted from the normal direction of the wall. An external time-periodic force always parallel to the wall acts on the particle. Similar to the follower force model, the component of the driving force tangential to the filament tip in our model can induce filament buckling. The oscillation of the cilium is sustained by the component normal to the filament tip.
The benefit of introducing a spherical particle is twofold. First, the particle may be charged or carry a magnetic moment, allowing easier experimental realizations of a driving force applied at the filament tip using electric or magnetic fields. Second, the drag force of a spherical particle scales linearly as the particle radius 
$b$, and the characteristic filament force upon the fluid scales approximately as 
$L/\ln (L/a)$ (Cox Reference Cox1970), where 
$L$ is the filament length and 
$a$ is the filament radius. Therefore, it is possible that the flux generated by the particle is comparable to or larger than the flux due to the filament as long as 
$b/L \gtrsim 1/\ln (L/a)$, which is a small value for a slender filament. Previous studies on artificial cilia only focused on elastic filaments or films, the effect of an attached particle has not been considered.
We first model the composite cilium numerically using a slender body theory. We observe that the cilium generates a large net flux at large tilt angles, accompanied by a buckling instability. The flux generated by the particle is indeed larger than the flux generated by the filament. The trajectories that the particle traces out depend sensitively on the buckling direction of the filament. As the buckling direction reverses, the particle trajectory and the generated net flux show abrupt changes. We also develop a reduced segmental model with a finite number of degrees of freedom that can reproduce similar dynamic behaviours. Finally, we discuss a possible experimental realization and estimate the magnitude of relevant parameters.
2. Model for an artificial cilium
We assume that the driving force acting on the particle is along the 
$y$ direction with a simple sinusoidal form, 
$\boldsymbol {F}_0(t) = F_0(t)\hat {\boldsymbol {y}} = A\sin (2{\rm \pi} t/\tau )\hat {\boldsymbol {y}}$, with 
$A$ the amplitude and 
$\tau$ the period. We first model an elastic filament and then derive a segmental model using rigid segments.
2.1. Dynamics of the composite cilium
Consider a slender, inextensible and elastic filament (with the aspect ratio 
$\epsilon = a/L \ll 1$) and denote the filament position by 
$\boldsymbol {r}(s, t)$ with the arc length 
$s\in [-L/2, L/2]$. The unit tangent vector 
$\boldsymbol {p} = (\cos \theta, \sin \theta )$, with 
$\theta$ the tangent angle between 
$\boldsymbol {p}$ and 
$\hat {\boldsymbol {x}}$ (figure 1) and the unit normal vector 
$\boldsymbol {p}^{\perp } = (-\sin \theta, \cos \theta )$. The filament force per unit length is given by the Euler–Bernoulli elasticity 
$\boldsymbol {f} = -B\boldsymbol {r}_{ssss}+(T\boldsymbol {r}_s)_s$, where 
$B$ is the bending rigidity and 
$T$ is the filament tension. The filament force may be derived from the bending energy formulation with 
$T$ acting as a Lagrange multiplier to enforce the filament inextensibility. From a non-local slender body theory (Johnson Reference Johnson1980), the velocity of the filament is given by
\begin{equation} 8{\rm \pi} \mu (\boldsymbol{r}_t-\boldsymbol{U}_{{p}\to {f}}) = \varLambda[\,\boldsymbol{f}] + \mathcal{K}[\,\boldsymbol{f}], \end{equation}
where 
$\mu$ is the fluid viscosity and 
$\boldsymbol {U}_{{p}\to {f}}$ is the velocity generated by the particle at the filament. The local operator 
$\varLambda [\boldsymbol {f}] = [c(\boldsymbol {I}+\boldsymbol {r}_s\boldsymbol {r}_s)+2(\boldsymbol {I}-\boldsymbol {r}_s\boldsymbol {r}_s)] \boldsymbol{\cdot} \boldsymbol {f}$, where 
$c = |\ln (\epsilon ^2 e)|$. The local operator captures the local drag anisotropy. From its inversion, we get the perpendicular and parallel anisotropic friction coefficients, 
$\xi _{\perp } = 8{\rm \pi} \mu /(c+2)$ and 
$\xi _{\parallel } = 4{\rm \pi} \mu /c$. In the limit of infinitely slender filaments one recovers the resistive force theory with 
$\xi _{\perp }/\xi _{\parallel } \approx 2$. The integral operator 
$\mathcal {K}[\boldsymbol {f}]$ captures the non-local interaction within the filament, which is interpreted as a disturbance velocity by the filament in the presence of a no-slip wall in our numerical computation. The filament velocity is then determined by a balance of viscous drag and the filament force 
$\boldsymbol {f}$. We write 
$8{\rm \pi} \mu (\boldsymbol {r}_t - \boldsymbol {U}) = \varLambda [\boldsymbol {f}]$, where 
$\boldsymbol {U}$ includes 
$\boldsymbol {U}_{{p}\to {f}}$ and the contribution from the filament motion.
Figure 1. Schematic of a composite artificial cilium consisting of an elastic filament clamped at the base and a spherical particle attached at the filament tip. The cilium is driven by an external time-periodic force 
$\boldsymbol {F}_0(t)$ acting on the particle. (a) Elastic model, (b) segmental model.
The filament tip is clamped to the particle surface, and the motion of the particle is fully determined by the translation and rotation of the filament at 
$s=L/2$. The particle position 
$\boldsymbol {r}_{{p}} = (\boldsymbol {r}+b\boldsymbol {p})_{s=L/2}$ with 
$b$ the particle radius, and its velocity 
$\boldsymbol {v}_{{p}} = (\boldsymbol {r}_t+b\theta _t\boldsymbol {p}^{\perp })_{s=L/2}$. We compute the disturbance velocity 
$\boldsymbol {U}$ as
\begin{equation} \boldsymbol{U}(s) = \int_{{-}L/2}^{L/2} \boldsymbol{\mathsf{G}}_{\delta}[\boldsymbol{r}(s), \boldsymbol{r}(s')]\boldsymbol{\cdot} \boldsymbol{f}(s')\,{\rm d}s' + 6{\rm \pi}\mu\chi b \boldsymbol{\mathsf{G}}_{\delta}[\boldsymbol{r}(s), \boldsymbol{r}_{{p}}] \boldsymbol{\cdot} \boldsymbol{v}_{{p}}. \end{equation}
Here, to account for the effect of the no-slip boundary, we use the regularized Blake tensor for a three-dimensional flow 
$\boldsymbol{\mathsf{G}}_{\delta}$ with 
$\delta$ the regularization parameter (Blake Reference Blake1971; Ainley et al. Reference Ainley, Durkin, Embid, Boindala and Cortez2008); 
$\chi$ is a correction to the free-space Stokes drag (Happel & Brenner Reference Happel and Brenner2012), and 
$\chi$ increases as the particle approaches the no-slip boundary.
To ensure the filament inextensibility, we require 
$\boldsymbol {r}_s\boldsymbol{\cdot} \boldsymbol {r}_{st} = 0$. Differentiating 
$\boldsymbol {r}_t$ with respect to the arc length and taking the tangent component, we obtain the tension equation
\begin{equation} 2cT_{ss}-(c+2)\theta_{s}^{2}T ={-}8{\rm \pi}\mu\boldsymbol{U}_s\boldsymbol{\cdot} \boldsymbol{p}-6cB\theta_{ss}^{2}-(7c+2)B\theta_{s}\theta_{sss}+(c+2)B\theta_{s}^{4}. \end{equation}
The normal component of 
$\boldsymbol {r}_{st}$ gives the equation of 
$\theta$, 
$\theta _t = \boldsymbol {r}_{st}\boldsymbol{\cdot}\boldsymbol {p}^{\perp }$
\begin{align} 8{\rm \pi}\mu\theta_{t}+(c+2)B\theta_{ssss} = 8{\rm \pi}\mu \boldsymbol{U}_{s} \boldsymbol{\cdot} \boldsymbol{p}^{{\perp}} + (9c+6)B\theta_{s}^{2}\theta_{ss}+(3c+2)T_{s}\theta_{s}+(c+2)T\theta_{ss}. \end{align} We scale length on 
$L$, force on 
$B/L^2$ and time on the period of the driving force 
$\tau$. The resulting elastoviscous number is 
$\eta = L / (B\tau /8{\rm \pi} \mu )^{1/4}$, which compares the viscous force with the elastic force. The other two control parameters include the tilt angle 
$\theta _0$ and the ratio of the particle radius to the filament length, 
$\beta = b/L$.
2.2. Boundary conditions and numerical methods
The orientation of the filament at 
$s=-L/2$ is fixed with 
$\theta = \theta _0$. The filament also has a zero velocity at 
$s=-L/2$, i.e. 
$\boldsymbol {r}_t = 0$. Separating the normal and tangent components of the filament velocity, this condition leads to 
$(T_s + 3B\theta _s\theta _{ss})_{s=-L/2} = 0$ and 
$(B\theta _{sss} - B\theta _s^3 - \theta _s T)_{s=-L/2} = 0$. The force and torque balance equations on the particle are
$$\begin{gather} ({-}B\boldsymbol{r}_{sss}+T\boldsymbol{r}_s)_{s={L}/{2}} ={-}6{\rm \pi}\chi\mu b(\boldsymbol{v}_{{p}}-\boldsymbol{U}_{{f}\to{p}})+\boldsymbol{F}_0(t), \end{gather}$$
$$\begin{gather}(B\boldsymbol{r}_{ss}\times \boldsymbol{r}_s)_{s={L}/{2}} = 8{\rm \pi}\mu b^3 \omega_{{p}}\hat{\boldsymbol{z}}, \end{gather}$$
where 
$\boldsymbol {U}_{{f}\to {p}}$ is the velocity induced by the filament at the particle centre and captures the effect of filament motion on the particle motion, 
$8{\rm \pi} \mu b^3$ is the rotational drag coefficient of a spherical particle and the particle angular velocity 
$\omega _{{p}} = (\theta _t)_{s=L/2}$. Since both 
$\boldsymbol {v}_{{p}}$ and 
$\omega _{{p}}$ are determined by the motion of the filament tip, (2.5) and (2.6) are translated into boundary conditions for 
$\theta$ and 
$T$. We solve (2.3) and (2.4) numerically using a second-order finite difference scheme (Tornberg & Shelley Reference Tornberg and Shelley2004). The coupled tension and 
$\theta$ equations with the nonlinear boundary conditions are solved using Newton's method. More details on the numerical methods can be found in Appendix A.
2.3. Segmental model
We replace an elastic filament with three rigid segments linked by torsional springs at the joints (figure 1b). The torsional spring exerts a torque proportional to the relative angle deflection between neighbouring segments. The total length of the segments is fixed, 
$\sum _{j=1}^{3} L_j = L$, and the length ratios 
$\gamma _j = L_j/L$. The centreline of each segment 
$\boldsymbol {r}_j = \boldsymbol {r}_j^c + s_j\boldsymbol {p}_j$, where 
$\boldsymbol {r}_j^c$ is the centre-of-mass position and the unit tangent vector 
$\boldsymbol {p}_j = (\cos \theta _j, \sin \theta _j)$. We keep the clamped conditions by requiring that 
$\theta _1 = \theta _0$ and 
$\boldsymbol {p}_3$ passes through the particle centre. Integrating (2.1) for each segment along the arc length and ignoring the non-local interaction
\begin{equation} 8{\rm \pi}\mu L_j\dot{\boldsymbol{r}}_j^c = [c(\boldsymbol{I}+\boldsymbol{p}_j\boldsymbol{p}_j) +2(\boldsymbol{I}-\boldsymbol{p}_j\boldsymbol{p}_j)]\boldsymbol{\cdot}\boldsymbol{F}_j, \end{equation}
where 
$\boldsymbol {F}_j$ is the total filament force upon the fluid. The torque-free condition of segment 3 with respect to 
$J_2$ is
\begin{equation} K(\theta_2-\theta_3)\hat{\boldsymbol{z}} + \boldsymbol{\sigma}^{h}_3|_{J_2} + (L_3+b)\boldsymbol{p}_3\times ({-}6{\rm \pi}\chi\mu b\boldsymbol{v}_{{p}}+\boldsymbol{F}_0) = 0, \end{equation}
where 
$K$ is the elastic modulus of the spring, the particle velocity 
$\boldsymbol {v}_{{p}} = \dot {\boldsymbol {r}}^c_3+(L_3/2+b)\dot {\boldsymbol {p}}_3$ and the hydrodynamic torque acting on segment 3 about 
$J_2$ is computed as
\begin{equation} \boldsymbol{\sigma}_3^{{h}}|_{J_2} ={-} \int_{{-}L_3/2}^{L_3/2}(s_3+L_3/2)\boldsymbol{p}_3\times \boldsymbol{f}_3 \,{\rm d}s_3. \end{equation}
The torque balance equation of segments 2 and 3 about joint 
$J_1$ is
\begin{align}&K(\theta_1-\theta_2)\hat{\boldsymbol{z}} + \boldsymbol{\sigma}_2^{{h}}|_{J_1} + \boldsymbol{\sigma}_3^{{h}}|_{J_2} - L_2\boldsymbol{p}_2\times\boldsymbol{F}_3 + [L_2\boldsymbol{p}_2+(L_3+b)\boldsymbol{p}_3]\nonumber\\ &\quad \times({-}6{\rm \pi}\chi\mu b\boldsymbol{v}_{{p}}+\boldsymbol{F}_0) = 0,\end{align}
where 
$\boldsymbol {\sigma }_2^{{h}}|_{J_1}$ is the hydrodynamic torque acting on segment 2 about 
$J_1$, and the last term computes the torque of the driving force and the viscous drag of the particle. The system is closed by the constraints that the velocities of neighbouring segments at the joints are the same. We scale length on 
$L$, time on 
$\tau$, force on 
$K/L$ and torque on 
$K$. The resulting elastoviscous number 
$\eta = L/(K\tau /8{\rm \pi} \mu )^{1/3}$. Hereafter, we use dimensionless variables with the same notation for both the elastic model and the segmental model.
3. Results
3.1. Fluid pumping
The motion of the artificial cilium reaches a steady state after a few periods. When 
$\theta _0 = 0$, the cilium beating patterns are periodic and symmetric over one period (figure 2a). The particle follows a symmetric ‘figure-8’ trajectory. As a result, the time-averaged disturbance flow field shows no net fluid pumping along the 
$y$ direction (figure 2b). However, when 
$\theta _0 < 0$, we observe a distinct power (
$F_0 < 0$) and recovery stroke (
$F_0 > 0$). A typical example is shown in figure 2(c). The filament is stretched out during the power stroke (red). As 
$\boldsymbol {F}_0$ reverses, the filament is bent with large deformation and moves towards the 
$+y$ direction during the recovery stroke (blue). The particle traces out an asymmetric trajectory and the disturbance flow field shows a clear net flux along the 
$-y$ direction (figure 2d). The centre of mass of the filament also traces out an asymmetric trajectory, similar to the natural cilium (Brumley et al. Reference Brumley, Wan, Polin and Goldstein2014).
Figure 2. Time lapse of the artificial cilium over one forcing period in the steady state and the time-averaged disturbance flow field with 
$\eta =3.5$ and 
$\beta =0.08$ for (a,b) 
$\theta _0=0$ (see supplementary movie 1 available at https://doi.org/10.1017/jfm.2023.436), and (c,d) 
$\theta _0=-1.16$ (supplementary movie 2). The force amplitude 
$A=65$. Time runs from red to blue. The cyan lines in (a,c) indicate the initial orientations. The dark dashed lines show the particle trajectories and the green dashed lines show the centre-of-mass position of the filament. The colour bars in (c,d) show the magnitude of the flow speed and arrows indicate flow directions.
The pumping performance can be characterized by the flux of the disturbance flow field obtained by integrating the Blake tensor over the 
$x$–
$z$ plane perpendicular to the pumping direction (Liron Reference Liron1978). The resulting instantaneous flux due to a point force of unit strength along the 
$+y$ direction located at a distance 
$h$ from the no-slip wall is given by 
$h/{\rm \pi} \mu$. The net flux 
$Q$ of the composite cilium consists of the flux generated by the filament 
$Q_{{f}}$ and the flux generated by the particle 
$Q_{{p}}$
\begin{equation} Q = Q_{{f}} + Q_{{p}} = \frac{1}{{\rm \pi} \mu} \left\langle\int_{{-}L/2}^{L/2}x(s)f_y(s)\,{\rm d}s\right\rangle + 6b \langle \chi x_{{p}}\boldsymbol{v}_{{p}}\boldsymbol{\cdot} \hat{\boldsymbol{y}}\rangle, \end{equation}
where 
$\langle \rangle$ denotes time average over a period in steady state and the distance of the particle from the wall 
$x_{{p}} = (x+b\cos \theta )|_{s={L}/{2}}$. The particle flux 
$Q_{{p}} \approx 6b\langle x_p \boldsymbol {v}_p\boldsymbol{\cdot} \hat {\boldsymbol {y}} \rangle = 6b\langle x_p \dot {y}_p \rangle = 6bS/\tau$, where 
$S$ is the area swept by the particle (Osterman & Vilfan Reference Osterman and Vilfan2011). This states that 
$Q_{{p}}$ is proportional to the area enclosed by the non-reciprocal trajectory of the particle.
Figure 3(a) shows the effect of 
$\theta _0$ with fixed force amplitude 
$A=65$. As 
$|\theta _0|$ increases (tilting towards the 
$-y$ direction), both 
$|Q_{{f}}|$ and 
$|Q_{{p}}|$ increase. Since the area swept by the particle is much larger than the area swept by the centre-of-mass position of the filament (figure 2c), 
$|Q_{{p}}|$ is always larger than 
$|Q_{{f}}|$. Surprisingly, for sufficiently large 
$|\theta _0|$, 
$|Q_{{f}}|$ and 
$|Q_{{p}}|$ drop abruptly to smaller values. The time lapse of the filament deformation shows that the abrupt change of 
$Q_{{p}}$ is accompanied by a reversal of filament bending direction. As shown by the two insets in figure 3(a), the filament is bent downwards during the recovery stroke for 
$\theta _0 = -1.12$ and upwards for 
$\theta _0 = -0.96$. This difference in the filament deformation leads to a difference in the particle trajectories, and therefore an abrupt change in 
$Q_{{p}}$. A similar discontinuity is observed when varying the force amplitude 
$A$ with fixed 
$\theta _0=-0.9$ (figure 3b). As 
$A$ is increased, the flux increases. When 
$A$ becomes sufficiently large, both 
$|Q_{{f}}|$ and 
$|Q_{{p}}|$ jump to larger values with a reversal of the bending direction.
Figure 3. (a,b) Time-averaged flux 
$|Q_{{f}}|$ (left axis, dark triangles), 
$|Q_{{p}}|$ (left axis, dark squares) and the total flux 
$|Q|$ (left axis, dark circles) as a function of (a) the initial tilt angle 
$\theta _0$ and (b) the force amplitude 
$A$ with 
$\beta =0.08$ and 
$\eta =3.5$. The flux is scaled by 
$B/\mu L$. In (a) 
$A = 65$, and in (b) 
$\theta _0 = -0.9$. The deviations 
$|\Delta \theta |=|\bar {\theta }_{{m}}-\theta _0|$ are shown on the right axes by the blue diamonds. The insets show the time lapse of the artificial cilium for different values of (a) 
$\theta _0$ (supplementary movies 3 and 4) and (b) 
$A$ corresponding to the filled red and green squares. The initial orientations are marked by the cyan lines. Panel (c) shows 
$|Q_{{p}}|$ as a function of 
$\theta _0$ and 
$A$.
The filament bending direction is likely determined by deviation of cilium orientation at the start of the recovery stroke from the initial orientation, which is also the natural orientation at rest with no external force. We compute the deviation as 
$\Delta \theta = \bar {\theta }_{{m}}(t=n)-\theta _0$, where the average orientation in steady state 
$\bar {\theta }_{{m}} = \int _{-1/2}^{1/2}\theta (s)\,{\rm d}s$ and 
$n$ is an integer. The deviation 
$|\Delta \theta |$ is largest for the symmetric case with 
$\theta _0 = 0$ (figure 2a). Figure 3(a,b) shows that 
$|\Delta \theta |$ decreases as the cilium is more tilted (with fixed 
$A$) or as the driving force 
$A$ decreases (with fixed 
$\theta _0$). The filament changes from downward bending to upward bending as 
$|\Delta \theta |$ exceeds a critical value around 0.15 in figure 3(a) and 0.18 in figure 3(b). As an example, the case of 
$\theta _0 = -1.12$ (downward bending) in figure 3(a) shows smaller deviation than the case of 
$\theta _0 = -0.96$ (upward bending). A similar observation can be made by comparing the two insets in figure 3(b) with 
$A = 50.9$ and 
$A = 69.8$.
We then compute the density plot of 
$Q_{{p}}$ as functions of 
$\theta _0$ and 
$A$. The values of 
$A$ and 
$\theta _0$ are limited to avoid the cilium touching the no-slip wall. As shown in figure 3(c), two regions with large negative flux are identified at large 
$A$ and 
$\theta _0$ corresponding to upward and downward bending. The sharp boundary separating these two regions spans a wide range for 
$\theta _0 \lessapprox -0.8$ and 
$A \gtrapprox 55$.
3.2. Linear stability analysis
When a filament is under compression at its tip along the tangential direction, the filament buckles as the compression exceeds a critical value (Landau & Lifshitz Reference Landau and Lifshitz1986). To investigate whether the tangent component of 
$\boldsymbol {F}_0$ is large enough to induce buckling in our system, especially at large values of 
$|\theta _0|$, we perform a linear stability analysis on the composite cilium with a tangential compression force 
$\varGamma$ acting on the particle.
Consider a small deformation from an initially straight filament with 
$\theta _0 = 0$, then 
$x\sim s$, 
$\boldsymbol {p} \sim (1, y_x)$, and 
$\boldsymbol {p}^{\perp } \sim (-y_x, 1)$. We use non-dimensional equations and ignore the non-local interaction. The linearized tension equation is 
$T_{ss} = 0$ with 
$T_{s} = 0$ at 
$s=-1/2$ and 
$T+3/2c\beta \chi T_{s} = -\varGamma$ at 
$s=1/2$. This leads to 
$T = -\varGamma$. We linearize (2.1) as
\begin{equation} -\alpha \varGamma y_{xx}-\alpha y_{xxxx} = y_{t}, \end{equation}
where 
$\alpha = \eta ^{-4}(c+2)$. The boundary conditions are
\begin{equation} y = 0,\quad y_{x} = 0, \end{equation}
at 
$s=-1/2$, and
\begin{equation} y_{xx} + \beta^3\eta^{4} y_{xt} = 0,\quad y_{xxx} - 3/4\eta^{4}\beta\chi (y_{t}+\beta y_{xt}) = 0, \end{equation}
at 
$s=1/2$. We consider perturbations of the form 
$y(x,t) = \phi (x)e^{\lambda t}$, and solve the resulting eigenvalue problem numerically using centred finite differences in the bulk and sided differences at the boundaries. Figure 4 shows that the real part of the largest eigenvalue 
${Re}(\lambda )$ decreases first and then increases as 
$\varGamma$ is increased. The system becomes unstable if 
${Re}(\lambda ) > 0$. For 
$\beta = 0$, the critical value 
$\varGamma ^{\ast } \approx 37.6$, which agrees with the result given in De Canio et al. (Reference De Canio, Lauga and Goldstein2017). The critical value becomes smaller as 
$\beta$ is increased. For 
$\beta = 0.08$, 
$\varGamma ^{\ast } \approx 22.6$.
Figure 4. Linear stability analysis. (a) The real part of the largest eigenvalue 
${Re}(\gamma )$ as a function of 
$\varGamma$ for different values of 
$\beta$. (b) In steady state, the 
$y$-component position of the particle (dark curve, left axis), actuation force 
$F_0$ (red solid curve, right axis) and filament tangent force 
$F_{{tang}}$ at 
$s=1/2$ (red dashed curve, right axis) as a function of time with 
$A = 65$, 
$\eta = 3.5$, 
$\beta = 0.08$ and 
$\theta _0 = -1.12$. The green dashed lines indicate 
$F_0^{\ast }$, the value of 
$F_0$ at the maximum of 
$|F_{{tang}}|$ around 
$t \approx 7.1$.
We observe a signature of buckling instability in our simulations. As shown in figure 5(b), computed with 
$\theta _0=-1.12$ (corresponding to the left inset in figure 3a), although 
$F_0$ increases rapidly first during the recovery stroke starting at 
$t=7.0$, the particle remains almost fixed for a period of time with little change in its 
$y$-component position. Meanwhile, the magnitude of the filament tangent force at 
$s=1/2$, expressed as 
$F_{{tang}}|_{s=1/2} = [(-\boldsymbol {r}_{sss}+T\boldsymbol {r}_s)\boldsymbol{\cdot} \boldsymbol {p}]_{s=1/2}$, increases and reaches a maximum at 
$t\approx 7.1$. The filament then buckles and the particle moves towards the 
$+y$ direction with 
$F_{{tang}}$ released. Therefore, we use 
$F^{\ast } = F_0(t\approx 7.1)$ to estimate the compression force acting on the composite cilia as 
$\varGamma = |F^{\ast }\sin (\theta _0)| \approx 37.9$. We apply the same approximation to the case of 
$\theta _0 = -0.96$ (the right inset in figure 3a) and obtain 
$\varGamma \approx 41.0$. Both estimates are larger than 
$\varGamma ^{\ast }$, indicating that buckling instability indeed occurs at large 
$|\theta _0|$ and the abrupt change in the net flux is caused by a reversal of the filament buckling direction.
Figure 5. Results of the segmental model with 
$\beta =0.08$, 
$\eta = 3.5$, and 
$\gamma _j=1/3$. (a,b) The time-averaged flux generated by the particle 
$|Q_{p}|$ (dark squares, left axis) as a function of (a) 
$\theta _0$ and (b) 
$A$. In (a) 
$A = 18.0$, and in (b) 
$\theta _0 = -1.0$. The flux is scale by 
$K/\mu$. Insets show the time lapse of the artificial cilia for different values of (a) 
$\theta _0$ (supplementary movies 5 and 6) and (b) 
$A$ corresponding to the red and green squares. Cyan lines indicate initial orientations. The deviations 
$|\Delta \theta |=|\bar {\theta }_{{m}}-\theta _0|$ are shown on the right axes by the blue circles. (c) Relative deflection 
$\theta _3-\theta _2$ under a constant driving force 
$F_0 = 18.0$ as a function of time 
$t$ for different values of 
$\delta \theta$ with 
$\theta _0 = -1.0$.
3.3. Segmental model
In the segmental model, the generated flux is also dominated by the contribution from the particle, especially at large 
$|\theta _0|$. Similar to the elastic model, abrupt change in 
$Q_{{p}}$ is also observed when varying 
$\theta _0$. Figure 5(a) shows that 
$|Q_{{p}}|$ first increases monotonically as 
$|\theta _0|$ increases from 0. For sufficiently large 
$|\theta _0|$, 
$|Q_{{p}}|$ suddenly jumps to a smaller value, along with a reversal of the buckling direction during the recovery stroke: segment 2 and segment 3 buckle upwards (
$\theta _2 > \theta _3$, see inset with 
$\theta _0 = -1.0$) when 
$\theta _0 \gtrsim -1.0$ and downwards (
$\theta _2 < \theta _3$, see inset with 
$\theta _0 = -1.2$) when 
$\theta _0 \lesssim -1.0$. An abrupt change in 
$Q_{{p}}$ is also found when the force amplitude 
$A$ is changed (figure 5b). As 
$A$ exceeds a critical value around 22.5, segments 2 and 3 switch from a downward buckling to an upward buckling with a significant increase in 
$|Q_{{p}}|$. Comparing the two cases shown by the insets of figure 5(b) with 
$A=13.6$ and 
$20.6$, 
$|Q_{{p}}|$ is tripled with an apparent increase in the enclosed area by the particle trajectory. The average orientation of the cilium is 
$\bar {\theta }_{{m}} = (\theta _1+\theta _2+\theta _3)/3$, and the deviation from the natural orientation 
$\Delta \theta = \bar {\theta }_{{m}}-\theta _0$. The transition from downward buckling to upward buckling is accompanied with a pronounced increase in 
$|\Delta \theta |$.
Finally, to verify that the reversal of the buckling direction is indeed caused by the deviation from the natural orientation, we perform the following numerical experiment. We apply a constant driving force 
$F_0=18.0$ along the 
$+y$ direction and evolve the system for a time duration of 
$0.5$. The deviation from a straight line is varied: we set 
$\theta _1(t=0) = \theta _0$, 
$\theta _2(t=0) = \theta _0-\delta \theta$ and 
$\theta _3 (t=0)= \theta _0-2\delta \theta$, where 
$\delta \theta$ is the magnitude of the deviation. The buckling direction is characterized by the relative deflection between segments 2 and 3. As shown in figure 5(c), segments 2 and 3 buckle downwards with 
$\theta _3-\theta _2 > 0$ for 
$\delta \theta \lessapprox 0.13$. An abrupt change occurs when 
$\delta \theta \gtrapprox 0.13$, segments 2 and 3 buckle upwards with 
$\theta _3-\theta _2 < 0$ for most of the time and reach equilibrium positions as 
$t\to 0.5$.
4. Conclusions and discussion
In this work we have studied the dynamics of a spherical particle constrained by an elastic filament as a simple model of an artificial cilium at low 
$Re$. We constructed an elastic model using slender body theory and derived a reduced segmental model with linked rigid segments. We found that the particle trajectory is strongly non-reciprocal at large tilt angle due to the buckling of the filament and a net fluid pumping parallel to the no-slip wall is generated. The particle trajectory and the induced flux depend sensitively on the buckling direction of the filament. Using the segmental model, we demonstrated that as the deviation of the cilium orientation at the start of the recovery stroke from the natural orientation exceeds a threshold, a reversal of the buckling direction occurs, leading to abrupt changes in the particle trajectories and the net flux.
The composite cilium we proposed may be fabricated as a whole experimentally at millimetre scale or larger by moulding silicone elastomers (Hu et al. Reference Hu, Lum, Mastrangeli and Sitti2018; Lu et al. Reference Lu, Zhang, Yang, Huang, Fukuda, Wang and Shen2018; Gu et al. Reference Gu2020), such as Ecoflex and Sylgard 184. Magnetic microparticles, like NdFeB, may be embedded within the spherical particle to provide a net moment after pre-magnetization. The cilium can then be driven by an external oscillating non-uniform magnetic field. To check if the parameter ranges for the observed pumping behaviours are realistic, we perform order-of-magnitude estimates of relevant parameters. We assume the filament length 
$L = 2$ mm, the radius 
$a = 0.1$ mm and the particle radius 
$b=0.2$ mm. The bending rigidity 
$B = {\rm \pi}Y a^4/2 \sim 10^{-11}\,{\rm N}\,{\rm m}^2$ for Ecoflex, where 
$Y$ is the Young's modulus (Vaicekauskaite et al. Reference Vaicekauskaite, Mazurek, Vudayagiri and Skov2020). Using numbers reported in previous experiments (Hu et al. Reference Hu, Lum, Mastrangeli and Sitti2018), the magnetization, which depends on the mass ratio of the magnetic particles and the magnitude of the magnetic field, may reach 
$M \sim 5\times 10^5\,{\rm A}\,{\rm m}^{-1}$ for a field strength around 
$1$ T. The resulting magnetic moment of the particle 
$m = 4{\rm \pi} b^3 M/3 \sim 10^{-5}\,{\rm A}\,{\rm m}^2$. To be comparable to the buckling threshold, 
$m\delta \sim 10 B/L^2$, the field gradient 
$\delta \sim 1$ T/m, which is approximately one or two orders of magnitude smaller than the gradient around common rare-earth magnets and easily achievable.
The hydrodynamic interactions with the no-slip wall have a small effect on the overall pump performance but affect the transition points between the upward and downward buckling. We also performed limited simulations with different actuation profiles. Including a weak second harmonic generates a faster increase of the actuation force and a larger bending deformation of the filament during the recovery stroke, shifting the particle closer towards the wall. This leads to a larger swept area and improves the pump performance. For elastic filaments free to undergo three-dimensional motions, tangential compression along the filament can induce three-dimensional spinning (Ling et al. Reference Ling, Guo and Kanso2018). The filament in our model is confined to planar motion in the 
$x$–
$y$ plane and the stability against perturbations in the 
$z$ direction is not analysed, but we speculate that the component of the driving force normal to the filament tip may favour motions in the 
$x$–
$y$ plane. We only considered a single cilium, and the results may be quite different in a cilium array due to hydrodynamic interactions. With a phase-dependent driving force, different synchronization states may arise by varying the ratio of particle radius to the filament length, leading to a different pump performance (Kotar et al. Reference Kotar, Leoni, Bassetti, Lagomarsino and Cicuta2010, Reference Kotar, Debono, Bruot, Box, Phillips, Simpson, Hanna and Cicuta2013; Chakrabarti & Saintillan Reference Chakrabarti and Saintillan2019; Man & Kanso Reference Man and Kanso2020; Chakrabarti et al. Reference Chakrabarti, Fürthauer and Shelley2022; Kanale et al. Reference Kanale, Ling, Guo, Fürthauer and Kanso2022).
Supplementary movies
Supplementary movies are available at https://doi.org/10.1017/jfm.2023.436.
Acknowledgements
We thank J. Zhang for helpful comments and suggestions. We also thank the anonymous reviewers for their constructive comments and suggestions.
Funding
This work is supported by National Natural Science Foundation of China (grant nos 12275332, 12047503, and 12247130) and Chinese Academy of Sciences and Wenzhou Institute (grant no. WIUCASQD2023009). The computations of this work were conducted on the HPC cluster of ITP-CAS.
Declaration of interests
The authors report no conflict of interest.
Appendix A. Numerical methods for the elastic model
We solve the elastic model using a second-order finite difference method. We discretize the arc length with a uniform grid of size 
$\Delta s$, 
$s_j = j\Delta s-1/2$ with 
$j = 0, 1, \ldots, 1/\Delta s$, and denote the quantities at 
$s_j$ with a subscript 
$j$. We discretize time as 
$t_n = n\Delta t$ and denote with a superscript 
$n$ the quantities at the current time step 
$t_n$. Schematically, we write the non-dimensional 
$\theta$ equation (2.4) as
\begin{equation} \frac{3\theta^{n+1}-4\theta^n+\theta^{n-1}}{2\Delta t} + (\alpha\theta_{ssss} - 3\zeta \theta_s^2\theta_{ss} - \zeta T_s\theta_s - \alpha T\theta_{ss})^{n+1} = [\boldsymbol{U}_s\boldsymbol{\cdot} \boldsymbol{p}^{{\perp}}]^{n}, \end{equation}
where 
$\alpha = \eta ^{-4}(c+2)$, 
$\zeta = \eta ^{-4}(3c+2)$. Denoting the solutions at the 
$k$th Newton interaction with a superscript 
$k$, and linearizing (A1) around current guesses,
$$\begin{align} &\delta\theta -
\frac{2\Delta t \alpha}{3}\alpha\theta_{ss}^k\delta T +
\frac{2\Delta
t}{3}({-}6\zeta\theta_s^k\theta_{ss}^k-\zeta
T_s^k)\delta\theta_s + \frac{2\Delta
t}{3}[{-}3\zeta(\theta_s^k)^2-\alpha
T^k]\delta\theta_{ss}\nonumber\\ &\quad + \frac{2\Delta t
\alpha}{3}\delta\theta_{ssss} - \frac{2\Delta t
\zeta}{3}\theta_s^k\delta T_s = \frac{4}{3}\theta^n -
\frac{1}{3}\theta^{n-1} + \frac{2\Delta
t}{3}[\boldsymbol{U}_s\boldsymbol{\cdot}
\boldsymbol{p}^{{\perp}}]^{n} + G[\theta^k, T^k],
\end{align}$$
where 
$G[\theta ^k, T^k]$ collects terms evaluated at iteration 
$k$. The tension equation (2.3) is linearized as
$$\begin{gather} 2c\delta T_{ss}-(c+2)(\theta^k_s)^2\delta T + [{-}2(c+2)\theta_s^k T^k +(7c+2)\theta^k_{sss}-4(c+2)(\theta^k_s)^3]\delta\theta_s\nonumber\\ +\, 12c\theta^k_{ss}\delta\theta_{ss} + (7c+2)\theta^k_s\delta\theta_{sss} ={-}\eta^4 [\boldsymbol{U}_s\boldsymbol{\cdot} \boldsymbol{p}]^n + M[\theta^k, T^k], \end{gather}$$
where 
$M[\theta ^k, T^k]$ collects terms evaluated at iteration 
$k$. The boundary conditions are linearized in a similar way. Results from previous time steps are used as the initial guesses. Solving the resulting linear system for 
$\delta \theta$ and 
$\delta T$ and iterating until converge, we obtain 
$\theta ^{n+1}$ and 
$T^{n+1}$. For most of our simulations, 
$\Delta s = 10^{-2}$, 
$\Delta t = 5\times 10^{-4}$, and the regularization parameter of the Blake tensor 
$\delta = 0.03$.
