Transportation efficiency of hydrodynamically coupled spherical oscillators in low Reynolds number fluids

The cilia and flagella are the fundamental organs for microorganisms as the instruments of swimming in low Reynolds number environments [1, 2] and transporting materials like protein [3]. The undergoing mechanics in their wave-like behaviors is also a perfect paradigm for the study of active matter [4]. The synchronization of the beating pattern of these filaments is common [5]. Usually, for example, bacteria perform a two-gait “run-and-tumble” motion in fluids [6]. From repeating such motion the bacteria gain the net momentum in order to move. Thus it is worth investigating the mechanism behind such behavior in an analytical interpretation.

Some experimental studies take efforts to discover the relation between viscous hydrodynamic interaction and synchronization, as Taylor originally proposed [7]. Experiments on carpets of bacteria with active flagella [8] and arrays of artificial magnetically actuated cilia [9, 10, 11] have found collective effects via hydrodynamic interactions such as complex flow patterns and collective phase shifts. Through the synchronization of filaments, it is well known that the energy dissipation of the model can be minimized [12], so one can expect that such collective effects provide sufficient influence over the performance of the oscillator model.

The dynamic of passively induced movement generated by the flagella is a highly complex elastohydrodynamic problem [13]. Therefore, in order to solve such a model, a computationally simple but still practically accurate model is required. It has been a popular idea of representing slender filaments by spheres for hydrodynamic modeling [14]. Uchida and Golestanian propose a generic criterion to reach stabilized states for two beads with arbitrary but periodic trajectory and force profile [15]. Recently, the details of hydrodynamic effects are categorized by direct and indirect effects from the slender structure of filaments, from approximating the filament as a point mass orbiting around the cell [16].

One representative of a minimal model can be found in Ref. [17]. In this literature, a pair of spherical beads orbiting in the same shape of periodic trajectory near a wall is investigated. At first, the linear stability of the synchronized state is analyzed and examined by, for instance, circular, linear, and elliptical trajectories. Furthermore, the analysis of nonlinear evolution for these trajectories is done in the far-field limit, as well as their near-field interaction between two beads due to the finite size of the trajectory.

For a better comprehension of the low Reynolds number hydrodynamics, there are literatures such as Kotar et al. [18], which focuses on the optimal condition of hydrodynamic synchronization for rotors in circular trajectories. This article numerically examines the model setups proposed by Niedermayer et al. [19] and Uchida & Golestanian [15], that a smaller spring stiffness or larger asymmetry leads to stronger coupling condition and synchronization. Another conclusion is, if some thermal noise is included in the original setup, the synchronization still holds for not very high temperatures. Liao & Lauga [20] have examined a model of two spherical bodies in circular orbits of late as well, and the condition of minimal energy dissipation rate is discovered.

In this article, we build a simple model, namely, a chain of spherical beads that oscillate in a one-dimensional orbit and are hydrodynamically coupled with one another, and focus on the transportation efficiency of the system, which is rarely focused on by previous literature. Our results are expected to contribute to describe the ciliary motion in mucociliary transport of the respiratory system and oviduct transport of the female reproductive system [21, 22, 23, 24]. Moreover, the metachronal pattern we discovered in the case of many beads provides theoretical implications for some motion generated from motile cilia in bacteria like Paramecium cells [25].

The flow of text is as the following: We will first introduce a minimal model, with only two beads inside. In the two-bead model, we provide the theoretical analysis of the dynamical system using a perturbation method, yielding expressions for net transportation, work, and efficiency. We then numerically analyze the case of many beads, indicating the positive cooperative effect via hydrodynamic interaction between the beads.

As an analytically tractable model of a flow generator, a one-dimensional setup is established and illustrated in Fig. 1. Specifically, consider $N$ spherical beads that can move only along the $x$-axis, obeying

$\displaystyle\gamma\left[\dot{x}_{i}-v(x_{i})\right]=F_{i},$

(1)

where $\gamma=6\pi\mu a$ is the drag coefficient for the fluid viscosity $\mu$, $a$ is the radius of the beads, $x_{i}$ ($i=1,2,\ldots,N)$ is the position of bead $i$, $v(x_{i})$ is the flow field at $x_{i}$ generated by the beads $j\neq i$, and $F_{i}$ is the active force applied to bead $i$. The flow field is given by

$\displaystyle v(x_{i})\coloneqq\sum_{j\neq i}H_{ij}F_{j},$

(2)

where $H_{ij}$ $(1\leq i,j\leq N)$ is the Oseen tensor for the Stokeslet in incompressible Newtonian fluid [26]. Because only one-dimensional motion in the $x$-direction is considered, the Oseen tensor is given as

$\displaystyle H_{ij}=\frac{1}{\gamma}\frac{3a}{2|x_{ij}|}+O\left[\left(\frac{a}{d}\right)^{3}\right],$

(3)

where $x_{ij}=x_{i}-x_{j}$ and $d$ is the distance between the equilibrium positions of the neighboring beads. Hereafter, all $O\left[\left(\frac{a}{d}\right)^{3}\right]$ terms are neglected. Also consider a test particle of radius $a$, obeying

$\displaystyle\gamma\left[\dot{x}_{0}-v(x_{0})\right]=0,$

(4)

where

$\displaystyle v(x_{0})\coloneqq\sum_{j=1}^{N}H_{0j}F_{j},$

(5)

$v(x_{0})$ is the velocity profile of the test particle. Denote the typical distance between the flow generator and the test particle by $l$. Here, we assume $l\gg d$ and employ the following approximation throughout this study:

$\displaystyle H_{0j}=\frac{1}{\gamma}\frac{3a}{2l}.$

(6)

With this approximation, $l$ only affects the magnitude of displacement of $x_{0}$ and can be chosen arbitrarily. Active force $F_{i}$ is modeled as follows. Assume that each bead is connected to a pinned point by an elastic beam, which exerts force by periodically changing its natural length. Specifically, consider

	$$\displaystyle F_{i}=k\left[x_{i}^{*}+L(\phi_{i})-x_{i}\right],$$		(7)
	$$\displaystyle L(\phi_{i})=\epsilon L_{0}\sin{\phi_{i}},$$		(8)
	$$\displaystyle\phi_{i}=\omega t+\psi_{i}$$		(9)

where $x_{i}^{*}+L(\phi_{i})$ is the equilibrium position, $L(\phi_{i})$ is the change of the natural length, $\phi_{i}$ is the phase of beating, $k$ is the spring constant, $\psi_{i}$ is the phase offset, and $\epsilon L_{0}$ is the oscillation amplitude of the natural length, $\epsilon$ is a nondimensional quantity introduced for later convenience. Without loss of generality, set $x_{i}^{*}=(i-1)d$ and $\psi_{1}=0$. One can interpret that beams repeat contraction and expansion with period $T=\frac{2\pi}{\omega}$ with some phase differences between the neighbors. In the following numerical simulations, without specific mentions we fix numerical values of parameters as: $a=0.1$, $\gamma=1,l=1000,L_{0}=1,\omega=1$ and vary $N,d,\epsilon$ and $\psi_{i}$ $(i=2,\ldots,N)$. Typical behavior for $N=4$ is displayed in Fig. 2.

To quantify the system performance, define the following efficiency for the transportation of the test particle:

$\displaystyle\hat{\eta}\coloneqq\frac{\gamma\big{<}v_{0}\big{>}^{2}}{P}=\frac{\gamma R^{2}}{WT},$

(10)

where $\big{<}\cdot\big{>}$ denotes the average over one cycle, $P$ and $W$ are the total power and work exerted by the flow generator for one cycle, respectively, and $R$ is the net transportation of the test particle for one period. More precisely, define

$\displaystyle R\coloneqq\lim_{n\to\infty}\int_{nT}^{(n+1)T}\dot{x}_{0}\ dt=\lim_{n\to\infty}\int_{nT}^{(n+1)T}v(x_{0})\ dt=\lim_{n\to\infty}\frac{3a}{2\gamma l}\int_{nT}^{(n+1)T}\sum_{j=1}^{N}F_{j}\ dt,$

(11)

where $n$ is the number of cycles and consider a large $n$ limit to eliminate the transient behavior. Also define

$\displaystyle W\coloneqq\sum_{i=1}^{N}W_{i},$

(12)

where

$\displaystyle W_{i}\coloneqq\lim_{n\to\infty}\int_{nT}^{(n+1)T}\big{(}\dot{x}_{i}-v(x_{i})\big{)}F_{i}\ dt=\lim_{n\to\infty}\frac{1}{\gamma}\int_{nT}^{(n+1)T}F_{i}^{2}\ dt.$

(13)

One may evaluate $R$, $W$, and thus $\hat{\eta}$ if the time-asymptotic behavior of $x_{i}(t)$ $(i=1,\ldots,N)$ is obtained.

We are interested in the cooperative effect on the efficiency and thus later compare the efficiency for different $N$ values. Thus, it is convenient to introduce the per-unit efficiency:

$\displaystyle\eta=\frac{\hat{\eta}}{N}=\frac{\gamma R^{2}}{WTN}.$

(14)

A positive cooperative effect is indicated if $\eta$ significantly increases with $N$.

A similar kind of efficiency could be found, such as the Lighthill hydrodynamic efficiency (cf. Ref. [27]) which happens to be one of the most popular among literatures (for instances, in Ref. [28, 29]). There, the particle is approximately spherical and deformable in order to move by morphing and the expression is

$\displaystyle\hat{\eta}^{\prime}\coloneqq\frac{\gamma\big{<}v_{0}\big{>}^{2}}{\big{<}P\big{>}}$

(15)

where $v_{0}$ is the velocity of the particle and $P$ is the work done from morphing. This is fundamentally the same definition as ours if one replaces the work done from morphing to the work done by beads.

Here, by considering the minimal model, the flow generator consisting of only two beads $(N=2)$, and perturbatively solving such system, we derive the expressions for $R,W$, and $\hat{\eta}$. For simplicity, $\psi_{2}$ is denoted as $\psi$ in this section; thus, $\phi_{1}=\omega t$ and $\phi_{2}=\omega t+\psi$. Then, our model for $N=2$ is

	$\displaystyle\dot{x}_{1}=\frac{3a}{2\gamma|x_{2}-x_{1}|}F_{2}+\frac{1}{\gamma}F_{1},$		(16a)
	$\displaystyle\dot{x}_{2}=\frac{3a}{2\gamma|x_{2}-x_{1}|}F_{1}+\frac{1}{\gamma}F_{2}.$		(16b)

The perturbation method might be applied by treating $\epsilon$ as a small parameter. Firstly introduce the dimensionless version of canonical basis as

	$\displaystyle X:=\frac{x_{1}+x_{2}}{d},$		(17)
	$\displaystyle Y:=\frac{x_{2}-x_{1}}{d}.$		(18)

For simplicity, always assume $Y>0$, which is a reasonable assumption physically, meaning that the beads never overlap with each other, so no collision occurs. Then, Eq. (16) becomes

	$\displaystyle\dot{X}$	$\displaystyle=\frac{k}{\gamma}\left(\frac{3a}{2dY}+1\right)\left[1+\epsilon\hat{L}_{0}(\sin{\phi_{1}}+\sin{\phi_{2}})-X\right],$		(19a)
	$\displaystyle\dot{Y}$	$\displaystyle=\frac{k}{\gamma}\left(-\frac{3a}{2dY}+1\right)\left[1+\epsilon\hat{L}_{0}(\sin{\phi_{2}}-\sin{\phi_{1}})-Y\right],$		(19b)

where $\hat{L}_{0}:=\frac{L_{0}}{d}$. Now, solve this using the perturbative forms:

	$\displaystyle X(t)$	$\displaystyle=1+\epsilon X_{1}(t)+\epsilon^{2}X_{2}(t)+O(\epsilon^{3}),$		(20)
	$\displaystyle Y(t)$	$\displaystyle=1+\epsilon Y_{1}(t)+\epsilon^{2}Y_{2}(t)+O(\epsilon^{3}).$		(21)

Note

$\displaystyle\frac{1}{Y}=1-\epsilon Y_{1}+\epsilon^{2}(-Y_{2}+Y_{1}^{2})+O(\epsilon^{3}).$

(22)

Therefore, by collecting $O(\epsilon)$ in Eq. (19),

	$\displaystyle\dot{X}_{1}={}$	$\displaystyle\alpha\left[m(t)-X_{1}\right],$		(23a)
	$\displaystyle\dot{Y}_{1}={}$	$\displaystyle\beta\left[h(t)-Y_{1}\right],$		(23b)

where $\alpha:=\frac{k}{\gamma}(\frac{3a}{2d}+1)$, $\beta:=\frac{k}{\gamma}(-\frac{3a}{2d}+1)$, $m(t):=\sin{\phi_{1}}+\sin{\phi_{2}}$, and $h(t):=\sin{\phi_{2}}-\sin{\phi_{1}}$. Similarly, for $O(\epsilon^{2})$,

	$\displaystyle\dot{X}_{2}$	$\displaystyle=-\frac{3akY_{1}}{2\gamma d}\left[m(t)-X_{1}\right]-\alpha X_{2},$		(24a)
	$\displaystyle\dot{Y}_{2}$	$\displaystyle=\frac{3akY_{1}}{2\gamma d}\left[h(t)-Y_{1}\right]-\beta Y_{2}.$		(24b)

Solving these equations, the time-asymptotic expressions are obtained as

	$\displaystyle{X}_{1}=$	$\displaystyle\frac{2\alpha\cos\left(\frac{\psi}{2}\right)\left[\alpha\sin\left(\omega t+\frac{\psi}{2}\right)-\omega\cos\left(\omega t+\frac{\psi}{2}\right)\right]}{\alpha^{2}+\omega^{2}},$		(25a)
	$\displaystyle{Y}_{1}=$	$\displaystyle\frac{2\beta\sin\left(\frac{\psi}{2}\right)\left[\beta\cos\left(\omega t+\frac{\psi}{2}\right)+\omega\sin\left(\omega t+\frac{\psi}{2}\right)\right]}{\beta^{2}+\omega^{2}},$		(25b)
	$\displaystyle{X}_{2}=$	$\displaystyle-\frac{3a\beta k\omega\sin(\psi)}{2\alpha\gamma d\left(\alpha^{2}+\omega^{2}\right)\left(\alpha^{2}+4\omega^{2}\right)\left(\beta^{2}+\omega^{2}\right)}$
		$\displaystyle\left[\left(\alpha^{2}+4\omega^{2}\right)\left(\alpha\beta+\omega^{2}\right)+\alpha\omega\left(\alpha^{2}+3\alpha\beta-2\omega^{2}\right)\sin(2\omega t+\psi)+\alpha\left(\alpha^{2}\beta-3\alpha\omega^{2}-2\beta\omega^{2}\right)\cos(2\omega t+\psi)\right],$		(25c)
	$\displaystyle{Y}_{2}=$	$\displaystyle-\frac{3a\beta k\omega\sin^{2}\left(\frac{\psi}{2}\right)\left[\left(2\omega^{3}-4\beta^{2}\omega\right)\cos(2\omega t+\psi)+\beta\left(\beta^{2}-5\omega^{2}\right)\sin(2\omega t+\psi)\right]}{\gamma d\left(\beta^{2}+\omega^{2}\right)^{2}\left(\beta^{2}+4\omega^{2}\right)}.$		(25d)

With these expressions, one may now calculate $R$ and $W$. Substituting Eq. (25) into Eq. (11),

$\displaystyle R=\frac{9\pi a^{2}k^{2}\sin(\psi)\beta\left(\alpha\beta+\omega^{2}\right)}{2l\gamma^{2}\alpha\left(\alpha^{2}+\omega^{2}\right)\left(\beta^{2}+\omega^{2}\right)}\epsilon^{2}\left[1+O(\epsilon)\right]$

(26)

which vanishes in the first order of $\epsilon$ while generally not in the second order of $\epsilon$ (cf. Eq. (10) and Eq. (25)). Substituting Eq. (25) into Eq. (12) and using $F_{1}^{2}+F_{2}^{2}=\frac{(F_{1}+F_{2})^{2}+(F_{1}-F_{2})^{2}}{2}$,

$\displaystyle W=$

$\displaystyle\frac{\pi k^{2}\omega d^{2}\left[\left(\beta^{2}-\alpha^{2}\right)\cos(\psi)+\alpha^{2}+\beta^{2}+2\omega^{2}\right]}{\gamma\left(\alpha^{2}+\omega^{2}\right)\left(\beta^{2}+\omega^{2}\right)}\epsilon^{2}\left[1+O(\epsilon)\right]$

(27)

which implies similar minimal energy dissipation versus phase difference as in Ref. [20]. Therefore,

$\displaystyle\eta$

$\displaystyle=\frac{81a^{4}\beta^{2}k^{2}\sin^{2}(\psi)\left(\alpha\beta+\omega^{2}\right)^{2}}{16\gamma^{2}\alpha^{2}d^{2}l^{2}\left(\alpha^{2}+\omega^{2}\right)\left(\beta^{2}+\omega^{2}\right)\left[\left(\beta^{2}-\alpha^{2}\right)\cos(\psi)+\alpha^{2}+\beta^{2}+2\omega^{2}\right]}\epsilon^{2}\left[1+O(\epsilon)\right].$

(28)

It should be noted that all quantities start from $O(\epsilon^{2})$ because the contributions of first-order terms are averaged out over one cycle. As shown in Fig. 3, these expressions are in excellent agreement with the simulation results.

Furthermore, by recalling that $\frac{a}{d}$ is assumed to be small and using $\alpha=\beta=\frac{k}{\gamma}+O(\frac{a}{d})$, the expressions can be much simplified as

	$\displaystyle R$	$\displaystyle=\frac{9\pi a^{2}L_{0}^{2}k^{2}\sin(\psi)}{2ld^{2}(k^{2}+\gamma^{2}\omega^{2})}\epsilon^{2}\left[1+O\left(\epsilon,\frac{a}{d}\right)\right],$		(29a)
	$\displaystyle W$	$\displaystyle=\frac{2\pi L_{0}^{2}k^{2}\omega\gamma}{\left(k^{2}+\gamma^{2}\omega^{2}\right)}\epsilon^{2}\left[1+O\left(\epsilon,\frac{a}{d}\right)\right],$		(29b)
	$\displaystyle\eta$	$\displaystyle=\frac{81a^{4}L_{0}^{2}k^{2}\sin^{2}(\psi)}{32l^{2}d^{4}\left(k^{2}+\gamma^{2}\omega^{2}\right)}\epsilon^{2}\left[1+O\left(\epsilon,\frac{a}{d}\right)\right].$		(29c)

From these expressions, observe that $R$ decreases with increasing $d$, proportionally to $1/d^{2}$. Thus, $R$ vanishes at a large $d$ limit, indicating the essential role of hydrodynamic interaction on the generation of the net flow. In contrast, $W$ is almost independent of $d$ and the leading term is just the sum of the works of two independent beads; the hydrodynamic interaction gives little contribution to the work for $d\gg a$. Therefore, $\eta$ is approximately proportional to $R^{2}$ and its maxima approximately coincide with the maxima of $R$.

In this section, a more general setup is considered, where the number of bead $N$ is no longer restricted by two. In the simulation, we set $\psi_{i}=i\Delta\psi$, where $\Delta\psi$ is a constant describing the phase difference. In this setup, one only needs to consider $0\leq\Delta\psi\leq\pi$ because of the symmetry of the system. Such a pattern mimics the metachronal wave ubiquitously found in the collective motion of flagella or cilia (cf. Ref. [30]). For $N=3$, as shown in Fig. 4, this setup is numerically proven to approximately produce the maximum efficiency output for the cases of $N=3$. The figure represents the efficiency as a function of $\psi_{2}$ and $\psi_{3}$, in which the optimality is found at $(\psi_{2},\psi_{3})\approx(1.3,2.6)$, thus $\psi_{3}\approx 2\psi_{2}$.

We first observe the dependence of per-unit efficiency $\eta$ on $\Delta\psi$, shown in Fig. 5(a). It is seen that the optimal phase difference, denoted as $\Delta\psi^{*}$, decreases as $N$ increases. In Fig. 5(b), we further plot the dependence of $\Delta\psi^{*}$ on the bead number $N$, suggesting that $\Delta\psi^{*}$ approaches a certain value as $N$ increases. There, we test two different $d$ values and find that the dependence on $d$ is rather weak.

We then observe the $N$-dependencies of $R,W$ and $\eta$ values for $\Delta\psi=\Delta\psi^{*}$, denoted as $R^{*},W^{*}$, and $\eta^{*}$, respectively. We plot these values in Fig. 6, where the vertical axes are scaled as $R^{*}\frac{d^{2}}{\epsilon^{2}},W^{*}\frac{1}{\epsilon^{2}}$, and $\eta^{*}\frac{d^{4}}{\epsilon^{2}}$. This scaling is motivated from the theoretical predicted dependence on $d$ and $\epsilon$ for $N=2$, given in Eq. (29); i.e., $R\propto\frac{\epsilon^{2}}{d^{2}},W\propto\epsilon^{2}$, and $\eta\propto\frac{\epsilon^{2}}{d^{4}}$. Figure 6(a) indicates that $R$ increases approximately linearly with $N$, particularly for large $d$. Moreover, as shown in Fig. 6(b), $W^{*}$ increases virtually linearly with $N$, indicating that the work was approximately equal to the summation of independent beads. Assuming $R,W\propto N$, $\eta=\frac{R^{2}}{WN}$ should be independent of $N$; however, as shown in Fig. 6(c), $\eta$ is actually strongly dependent on $N$, indicating that the growth of $R$ significantly deviates from a linear growth. Importantly, $\eta$ for $N>2$ is larger than $\eta$ for $N=2$, indicating the positive cooperative effect of hydrodynamic coupling. It should also be noted that $\eta$ saturates or even turns to decrease for large $N$, which could come from our simple setting for the phase offsets and can possibly be improved by considering a more complex pattern of $\psi_{i}$. It is also observed that all the quantities at given $N$ in Fig. 6 have similar magnitudes, suggesting that the dependence $R\propto\frac{\epsilon^{2}}{d^{2}},W\propto\epsilon^{2}$, and $\eta\propto\frac{\epsilon^{2}}{d^{4}}$ roughly holds true for $N\geq 3$.

We introduced a simple one-dimensional coupled sphere motor model, with the generated flow and transportation efficiency investigated both theoretically and numerically. For the minimal case with two beads, the proper theoretical expression for both the net transportation and the efficiency are derived by the perturbation method, and it is proven that these results are in excellent agreement with the numerical results. Furthermore, we numerically investigated the case of $N\geq 3$ and found that the efficiency increases with $N$.

For a more practical usage, the effect of a wall, to which the system is pinned, should also be examined; such a setup can be treated by using the Blake tensor [31, 17] instead of the Oseen tensor. It could be very helpful to derive theoretical expressions for the net transportation and efficiency for this case. In our model, oscillators were designed non-autonomous for simplicity; to understand the role of synchronization and self-organized metachronal wave pattern, it is essential to extend the model to the case of autonomous oscillators and compose a theoretical framework. Such an effort would substantially contribute to deepening the understanding of the hydrodynamic effect of interactive objects in low Reynolds numbers. Moreover, it is worth discovering whether the higher order interactions will eventually destroy the positive coupling effect completely for very large number of oscillators, and if this is also true with effect of a wall. In addition, a two-dimensional movement model could be considered, as many other publications have focused, while the complexity in finding a theoretical expression could be fairly challenging.