text:
# Charge pumping in a quantum wire driven by a series of local time-periodic potentials

A parametric electron pump has attracted considerable attention in recent years.  It is a device that generates a dc current at zero bias by cyclic deformations of system parameters. The quantum pumping mechanism was originally proposed by Thouless, who studied the integrated particle current produced by a slow periodic variation of the potential, and showed that in a finite torus the integral of the current over a period can vary continuously, but it must have an integer value in an infinite periodic system with full bands. Such quantized charge transport was proposed to become an electric current standard.

Quite recently, the charge pumping was observed experimentally.  For technical reasons, instead of measuring charge currents, the pumped dc voltage $`V_{dot}`$ is measured in a quantum dot where two gates with oscillating voltages control the deformation of the shape of the dot. For weak pumping, the observed charge pumping has a sinusoidal dependence on the phase difference $`\varphi `$ between the two shape-distorting ac voltages applied to the gates, and is proportional to the square of pumping strength $`V`$. For strong pumping, the pumped current deviates from the square dependence on $`V`$ and becomes nonsinusoidal, being always antisymmetric about $`\varphi =\pi `$. The charge pumping may have a close relation to the adiabatic Berry’s phase since the evolution of the system is cyclic and is controlled by several system parameters, referred to as the parametric pumping. Based on this understanding, the total charge pumped per cycle is proportional to the area enclosed by the path in the parameter space, and nonzero pumping current requires at least two parameters.  The pumped charge drived by two parameters should be zero if two parameters are in phase ($`\varphi =0`$)) since the area enclosed by the path is zero. However, it is in contradiction with the observed current $`I(\varphi =0)0`$. One of possible mechanisms of nonzero currents for $`\varphi =0`$ is photovoltaic effects introduced in Ref., where a surprising result, nonzero dc current generated by a single pumping gate voltage, is also reported. The general physics of a quantum pump has been the subject of several theoretical analyses. Zhou et al demonstrated that at low temperatures both the magnitude and the sign of the pumped charges are sample specific quantities, and the typical value in disordered (chaotic) systems turns out to be determined by quantum interference effects. Another general expression for the average transmitted charge current was derived by Brouwer  under the adiabatic condition and based on the time-dependent S-matrix method, which appears to be quite successful for (adiabatic) weak pumping. Adiabaticity here means that the oscillating period $`\tau `$ of the system is much larger than the Wigner delay time $`\tau _w`$. Note that the adiabatic condition does not simply imply that the pumping strength $`V`$ should be very small. In fact, the adiabatic condition requires that $`\tau `$ must be larger as $`V`$ increases. On the other hand, the pumping was not weak in the experiments.  The main purpose of the paper is to develop a theory, which is also applicable in the case of strong pumping. By using the Floquet theorem, the photon-assisted transport has been taken into account. We calculate the pumped current through a mesoscopic region in the presence of time-periodic potentials. Our main results in the weak pumping regime, as well as those in the strong pumping regime are consistent with the experiment reported in Ref..

Consider electrons transmitting through a one-dimensional scattering region ranging from $`x_0`$ to $`x_0+\delta `$. The potential is given by

$$V(x,t)=\{\begin{array}{cc}0,& x_0<0,x>x_0+\delta \\ V_s(x,t),& x_0xx_0+\delta \end{array}$$
(1)
with $`V_s(x,t)=V_0+V_scos(\omega t+\varphi _s)`$. The Schrödinger equation can be written as

$$i\mathrm{}\frac{\mathrm{\Psi }(x,t)}{t}=\frac{\mathrm{}^2}{2m^{}}\frac{^2\mathrm{\Psi }(x,t)}{x^2}+V(x,t)\mathrm{\Psi }(x,t),$$
(2)
with $`m^{}`$ as the electron effective mass. Equation 2 can be solved by using the Floquet theorem. By setting $`\mathrm{\Psi }_{Fl}(x,t)=e^{iE_{Fl}t/\mathrm{}}\psi (x,t),`$ where $`E_{Fl}`$ is the Floquet eigenenergy and $`\psi (x,t)`$ is a periodic function $`\psi (x,t)=\psi (x,t+\tau )`$ with period $`\tau =2\pi /\omega `$, the Schrödinger equation takes the form

$$E_{Fl}\psi (x,t)=\frac{\mathrm{}^2}{2m^{}}\frac{^2\psi (x,t)}{x^2}i\mathrm{}\frac{\psi (x,t)}{t}+V(x,t)\psi (x,t).$$

Substituting $`\psi (x,t)=g(x)f(t)`$, we have two separated equations with an introduced constant $`E`$,

$`{\displaystyle \frac{\mathrm{}^2}{2m^{}}}{\displaystyle \frac{^2}{x^2}}g(x)+V_0g(x)`$ $`=`$ $`Eg(x),`$ (3)
$`i\mathrm{}{\displaystyle \frac{f(t)}{t}}V_scos(\omega t+\varphi _s)f(t)`$ $`=`$ $`(EE_{Fl})f(t).`$ (4)
Integrating Eq.(4) gives

$$f(t)=e^{\frac{iV_ssin\varphi _s}{\mathrm{}\omega }\frac{i(EE_{Fl})t}{\mathrm{}}}\underset{\nu =\mathrm{}}{\overset{\mathrm{}}{}}e^{\nu \varphi _s}J_\nu (\frac{V_s}{\mathrm{}\omega })e^{i\nu \omega t},$$
(5)
where $`J_\nu (x)`$ is the Bessel function of the first kind of order $`\nu `$. Since $`f(t)`$ is periodic in time with period $`\tau `$, it follows from Eq.(5) that $`EE_{Fl}=m\omega `$ with $`m`$ as an integer. The equation for $`g(x)`$ has a solution

$$g(x)=e^{\pm ik_m^sx},(k_m^s)^2=2m^{}(E_{Fl}+m\mathrm{}\omega V_0)/\mathrm{}^2.$$
(6)
Thus $`\psi (x,t)`$ becomes

$$\psi _m(x,t)=e^{i(V_s/\mathrm{}\omega )sin\varphi _s\pm ik_m^sx}\underset{n}{}F_{nm}e^{in\omega t},$$
(7)
with $`F_{nm}=exp[i(nm)\varphi _s]J_{nm}(V_s/\mathrm{}\omega )`$.

We consider an incoming wave from the left with the energy $`E_0=\mathrm{}^2k_0^2/2m^{}`$, then the outgoing waves should be divided into different modes $`E_n`$, which satisfies $`E_n=E_0+n\mathrm{}\omega `$ with $`n=0,\pm 1,\pm 2,\mathrm{}`$. The propagating modes mean that $`E_n>0`$, while the evanescent modes mean that $`E_n0`$. The latter exists only in the neighborhood of the oscillating barrier and do not propagate. Denote $`k_n=\sqrt{2m^{}E_n}/\mathrm{}`$, the solutions of the Schrödinger equation can be written as

$`\mathrm{\Psi }_l`$ $`=`$ $`{\displaystyle \underset{n=\mathrm{}}{\overset{\mathrm{}}{}}}(A_n^ie^{ik_nx}+A_n^oe^{ik_nx})e^{iE_nt/\mathrm{}},(x<x_0)`$
$`\mathrm{\Psi }_s`$ $`=`$ $`e^{\frac{iE_{Fl}t}{\mathrm{}}}{\displaystyle \underset{m,n=\mathrm{}}{\overset{\mathrm{}}{}}}(a_me^{ik_m^sx}+b_me^{ik_m^sx})F_{nm}e^{in\omega t},`$
$`(`$ $`x_0xx_0+\delta )`$
$`\mathrm{\Psi }_r`$ $`=`$ $`{\displaystyle \underset{n=\mathrm{}}{\overset{\mathrm{}}{}}}(B_n^ie^{ik_nx}+B_n^oe^{ik_nx})e^{iE_nt/\mathrm{}},(x>x_0+\delta )`$
where $`A_n^i`$ and $`B_n^i`$ are the probability amplitudes of the incoming waves from the left and right, respectively, while $`A_n^o`$ and $`B_n^o`$ are those of the outgoing waves. We can characterize the barrier by a scattering matrix $`S`$ which is a matrix connecting the incoming and outgoing channels

$$\left(\begin{array}{c}𝐀^o\\ 𝐁^o\end{array}\right)=S\left(\begin{array}{c}𝐀^i\\ 𝐁^i\end{array}\right),$$

where the $`S`$ matrix can be derived by the matching conditions for the wave function $`\mathrm{\Psi }(x,t)`$ and its derivative $`_x\mathrm{\Psi }(x,t)`$ at $`x=x_0`$ and $`x=x_0+\delta `$. After eliminating $`a_m`$ and $`b_m`$, we have

$$S=\left(\begin{array}{cc}R_{}& T_{}\\ T_{}& R_{}\end{array}\right),$$
(8)
where $`T_{}=L_LTL_R^1,`$ $`R_{}=L_R^1RL_R^1,`$ $`T_{}=L_R^1TL_L,`$ and $`R_{}=L_LRL_L`$. Here the left (right) arrow indicates incoming waves from right (left), the matrices $`L_L`$ and $`L_R`$ are defined as $`(L_L)_{mn}=exp[ik_nx_0]\delta _{mn}`$ and $`(L_R)_{mn}=exp[ik_n(x_0+\delta )]\delta _{mn}`$. $`T`$ and $`R`$ are given by

$`T`$ $`=`$ $`(C_1^1D_1+C_2^1D_2)/2,`$ (9)
$`R`$ $`=`$ $`(C_1^1D_1C_2^1D_2)/2,`$ (10)
where

$`C_1`$ $`=`$ $`(L_s\stackrel{~}{I})K_sF^{}(L_s+\stackrel{~}{I})F^{}K,`$
$`D_1`$ $`=`$ $`(L_s\stackrel{~}{I})K_sF^{}(L_s+\stackrel{~}{I})F^{}K,`$
$`C_2`$ $`=`$ $`(L_s+\stackrel{~}{I})K_sF^{}(L_s\stackrel{~}{I})F^{}K,`$
$`D_2`$ $`=`$ $`(L_s+\stackrel{~}{I})K_sF^{}+(L_s\stackrel{~}{I})F^{}K,`$
with the matrices $`(L_s)_{mn}=exp[ik_n^s\delta ]\delta _{mn}`$, $`(K_s)_{mn}=k_n^s\delta _{mn}`$, $`K_{mn}=k_n\delta _{mn}`$, $`\stackrel{~}{I}`$ as the unit matrix and $`F^{}`$ as the Hermitian conjugate of $`F`$. The electronic transport properties of the scattering region may be obtained straightforward from Eq.(8).

The above method may be generalized to $`l`$ time-periodic barriers described by

$$V(x,t)=\{\begin{array}{cc}0,& x<0,x>a_l\\ V_1(x,t),& 0x<a_1\\ V_2(x,t),& a_1x<a_2\\ \mathrm{}& \mathrm{}\\ V_l(x,t),& a_{l1}xa_l\end{array}$$
(11)
where $`V_1(x,t)=V_1^0+V_1cos(\omega _1t+\varphi _1)`$, $`V_2(x,t)=V_2^0+V_2cos(\omega _2t+\varphi _2),\mathrm{}`$ and $`V_l(x,t)=V_l^0+V_lcos(\omega _lt+\varphi _l)`$. This potential may be more a realistic model for experiments. Obviously the transport properties for each barrier can be characterized by an $`S`$ matrix given by

$$S^\alpha =\left(\begin{array}{cc}R_{}^\alpha & T_{}^\alpha \\ T_{}^\alpha & R_{}^\alpha \end{array}\right),$$

where $`\alpha =1,2,\mathrm{},l`$, $`T_{}^\alpha `$, $`T_{}^\alpha `$, $`R_{}^\alpha `$ and $`R_{}^\alpha `$ can be derived by the same method presented above. Now the propagating mode $`E_n`$ should be replaced by

$$E(n_j)=E_0+\underset{n_j=\mathrm{}}{\overset{\mathrm{}}{}}n_j\mathrm{}\omega _j.(j=1,\mathrm{},l)$$
(12)
The associated transfer-matrix $`M^\alpha `$ for the $`\alpha `$th barrier may be derived directly from the $`S^\alpha `$ matrix

$$M^\alpha =\left(\begin{array}{cc}(T_{}^\alpha )^1& (T_{}^\alpha )^1R_{}^\alpha \\ R_{}^\alpha (T_{}^\alpha )^1& T_{}^\alpha R_{}^\alpha (T_{}^\alpha )^1R_{}^\alpha \end{array}\right).$$

The total transfer-matrix $`M^t`$ for all those barriers is determined by

$$M^t=\left(\begin{array}{cc}M_{11}^t& M_{12}^t\\ M_{21}^t& M_{22}^t\end{array}\right)=M^lM^{nl}\mathrm{}M^1,$$

where $`M_{ij}^t(i,j=1,2)`$ are the partitioned matrices with the same size as $`T_{}^\alpha `$. The total scattering matrix $`S^t`$ can be derived from $`M^t`$ as

$$S^t=\left(\begin{array}{cc}(M_{11}^t)^1M_{12}^t& (M_{11}^t)^1\\ M_{22}^tM_{21}^t(M_{11}^t)^1M_{12}^t& M_{21}^t(M_{11}^t)^1\end{array}\right).$$

In each cycle a net charge current may pass through the scattering region in the direction determined from the detailed form of $`S^t`$ matrix. We define a net transmission coefficient (for an incoming wave in mode $`E_0=E`$) by

$$T_{net}=\underset{E(n_j)>0}{}\sqrt{\frac{2E(n_j)}{m^{}}}(|T_{,n_j0}^t(E_0)|^2|T_{,n_j0}^t(E_0)|^2).$$

The average net current per period $`\tau `$ (for $`E_0`$) through those barriers is $`j(E_0)=T_{net}(E_0)`$. If the system is connected through two ideal leads to two electron reservoirs with the same chemical potential $`\mu `$, the average pumped current per period $`\tau `$ is given by

$$I(\mu )=e_0^{\mathrm{}}𝑑Eg(E)f(E\mu )T_{net}(E),$$
(13)
where $`g(E)=\sqrt{2m^{}/E}/h`$ is the density of electrons contributing to the current in one direction, and $`f(E\mu )`$ is the Fermi-Dirac distribution. At zero temperature, it becomes

$$I(\mu )=\frac{2e}{h}_0^\mu 𝑑E\sqrt{\frac{m^{}}{2E}}T_{net}(E).$$
(14)
Another important quantity is the Wigner delay time which gives the time delay of the scattered electron due to its interaction with the scattering field ( here the oscillating potential). It relates to the S-matrix by

$$\tau _w(E)=\frac{i\mathrm{}}{N_c}Tr[(S^t)^{}\frac{dS^t}{dE}]=\frac{i\mathrm{}}{N_c}\frac{d}{dE}ln(detS^t),$$
(15)
where $`N_c`$ is the number of open channels. Physically, the Wigner time represents the time spent by a wave packet passing through the scattering region. The charge pumping is supposed to be adiabatic when $`\tau `$ is much greater than the Wigner delay time $`\tau _w`$.

It is obvious that the net charge transfer in one cycle is zero for a single time-periodic barrier since $`T_{}=T_{}`$. Then the simplest system which may induce the nontrivial charge pumping should include at least two barriers. As an example we consider a mesoscopic system with two time-periodic barriers connected through ideal leads to two electron reservoirs with the same chemical potentials $`\mu `$. The potentials are described by $`V_1(x,t)=V_1^0+V_1cos(\omega t)`$, $`V_2(x,t)=V_2^0`$, and $`V_3(x,t)=V_3^0+V_3cos(\omega t+\varphi )`$. This appears to be a simplified model for the Switkes et al’s experiment, nevertheless it turns out that some essential characteristics can be exhibited, as we will address below. In the following numerical calculations, $`\mu =75mev`$, $`m^{}=0.067m_e`$ (with $`m_e`$ as the mass of the free electron), $`V_2^0=30mev`$, $`V_2=0`$, and $`N_c`$ is determined by a natural condition: $`|T_{}|^2+|R_{}|^21.0c_e`$ with $`c_e`$ ($`=1.0\times 10^4`$ in this paper) as a defined error.

The general characteristics of quantum pumping are shown in Figs.1 and 2. The parameters in Figs.1, 2, and 3 are chosen as $`V_1^0=V_3^0=50mev`$, and $`V_1=V_3=V`$. Figure 1 shows that the pumped current $`I(V)`$ is proportional to $`V^2`$ for small pumping amplitude $`V`$, with the proportional factor depending on the driving frequency ($`\mathrm{}=1`$). But it deviates from $`V^2`$-dependence for the strong pumping case. On the other hand, the pumped current is sinusoidal dependence on $`\varphi `$ for weak pumping, and becomes nonsinusoidal dependence on $`\varphi `$ when $`V`$ increases, as seen in Fig.2. Another important characteristic shown in Fig.2 is that $`I(\pi +\varphi ,V)=I(\varphi ,V)`$ for all amplitude strengths, and $`|I(\varphi ,V)|`$ is maximum at $`\varphi =\pi /2`$ or $`\varphi =3\pi /2`$ for weak pumping. Remarkably, these results for $`I(\varphi ,V)`$ are in agreement with the experimental observation in Ref..

Figure 3 shows that sharp peaks in the Wigner time occur at the resonance insert energies $`E_0=n\mathrm{}\omega `$. Besides, jumps in the pumped current as a function of $`E_0`$ appear at the peak positions of the Wigner time. The direction of the current depends crucially on the insert energy. It is interesting to note that the adiabatic condition is not necessary in our calculations. Figure 3 indicates that the maximum value of $`\tau _w`$ is about $`5.5ns`$ for $`\mathrm{}\omega =6.0mev`$ (corresponding to $`\tau 0.7ps`$ ), which is much greater than the pumping cyclic time $`\tau `$. Then we may say that the method described here is beyond adiabaticity. Actually, the nonadiabatic effects are only important for the strongly photon-assisted transport since $`\tau _w`$ is greater than $`\tau `$ only if the energy of the incoming wave is approximately equal to the resonance energy for photon-assisted tunneling. Physically, by emitting or absorbing photons, the outgoing waves may be at the quantum states different from that of the incoming wave. Consequently, the adiabatic condition, which requires that the quantum state is at the same instant state in the whole evolution, is not satisfied. Note that the formula derived by Brouwer may be valid merely under the adiabatic condition, and thus the method developed here may be quite useful.

It is quilt intriguing to note from Fig.4 that $`I(\varphi =0)`$ is nonzero for $`V_1V_3`$, while the corresponding areas enclosed by the path in the parameter space $`\{V_1(x,t),V_3(x,t)\}`$ are zero. Although the pumped currents in the above case were predicated to be zero under the adiabatic approximation, the deviation from zero is reported experimentally at strong pumping, just as we observed here in terms of a rigorous theoretical analysis which is also valid for strong pumping. Moreover, from comparison with that $`I(\varphi =0)=0`$ for $`V_1=V_3`$, it is now clear that the present nonzero pumped currents stem from the spatial asymmetry of potentials $`V_1V_3`$, which is coincident with the result obtained by Wagner in Ref.: the nonzero currents may be observed in a single osscillating potential but with asymmetric static potential. Actually, to observe a pump current at zero applied bias, it seems that the inversion symmetry should be broken, either in real or in $`k`$ space.

The fact that $`I(\varphi =0)`$ is nonzero at strong pumping may be understood based on a scenario of the nonadiabatic geometric phase. Pumped currents are determined by geometric phase accumulated in the evolution.  Under the adiabatic approximation, $`I(\varphi =0)=0`$ is predicted theoretically because the corresponding adiabatic geometric phase is zero. While it is now clear that the nonadiabatic geometric phase may be nonzero even in the case where the area enclosed by the path in the parameter space is zero (thus the adiabatic phase is zero). Therefore, the nonadiabatic correction to the currents should be taken into account for strong pumping whenever the adiabatic condition is not well satisfied. Physically, it is reasonable to believe that the observed nonzero pumping at phase $`\varphi =0`$ for the strong pumping stems from the nonadiabatic correction when the inversion symmetry is broken. Practically, the asymmetric spatial potential might be present in the experiment, which may originate from either the shape-distorting ac voltages, or from the internal potential established during transport. Since the current calculated in this approach is conserved since $`|T_{}|^2+|R_{}|^2=1.0`$, no internal potential appears explicitly in the present formulism. It is worth pointing out that nonzero pumped currents for $`\varphi =\pi `$ are also seen in Fig.4, which seems to contradict with that in Ref . Also note that a nonzero $`I(\varphi =\pi )`$ was also predicted by another totally different theoretical study, so this contradication is still an interesting open question at present.

In summary, we developed a method to calculate the pumped current and Wigner delay time in a mesoscopic system with a series of time-periodic barriers connected to two electron reservoirs, which appears to be applicable for strong pumping cases.

We thank the support from a RGC grant of Hong Kong(Grant No. HKU7118/00P) and a CRCG grant at the University of Hong Kong.


format: text

text:
# A simulation study of the dynamics of a driven filament in an Aristotelian fluid

## 1 Introduction

For a class of biologically important polymeric materials elasticity is crucial. Their typical lengths (microns or less) are comparable with the scale on which rigidity prevents them from collapsing. The cytoskeletal filaments actin and microtubules  fall in this category, as do cilia and flagella. The latter are motile assemblies of microtubules and other proteins. Because of their size and typical velocities, the motion of these filaments is nearly always in the low Reynolds number regime. This is an inertialess, Aristotelian, world where the dynamics of a surrounding fluid become time-reversible. As a notable consequence, it is difficult to generate any propulsion on this scale . Nonetheless, cytoskeletal filaments are involved in cellular and microorganism motility. Perhaps the most widely known example is that of the flagellum of a sperm cell, that enables it to swim along the ovaric tubes. The internal drive of a flagellum, however, is rather complicated . It involves many internal degrees of freedom and active components. On the other hand, modern micromanipulation techniques, such as optical and magnetic trapping, open up the possibility of perturbing otherwise passive filaments with a simplified and controlled drive. This provides a potentially useful model system for which one may study the fundamentals of motility.

It is this problem we concern ourselves with here. Specifically, we consider the flexible one-armed swimmer. That is, an elastic filament that is wiggled at one end. If the filament were rigid, the reversible motion of the surrounding fluid would ensure that this mechanism generates no propulsion (the “scallop theorem” as Purcell termed it). However, the flexibility of the arm breaks the time reversal symmetry for the motion of the assemblage. This makes propulsion, in principle, possible. For any microscopic filament the factors that determine its dynamic behavior are the same. Namely, the equations of motion will be essentially inertialess. The motion itself will be determined by a balance between forces driving the filament, friction forces exerted as the surrounding fluid opposes any motion, and bending forces that try to restore the (straight) equilibrium state. For relatively simple model systems, there has recently been theoretical progress in solving analytically the “hyperdiffusion” equation that, in the limit of small amplitude motion, describes the movement of such a filament. Wiggins and Goldstein  considered the motion of a single filament driven at one end by an external perturbation. Their analysis emphasized that there are two very different regimes; one where bending forces dominate and the filament behaves like a rigid rod, and a second where the viscous damping of the fluid has the effect of suppressing the propagation of elastic waves. For the one armed swimmer, this leads to an optimal set of parameters that maximize either the swimming speed or swimming efficiency. The same analysis gives predictions for the shape of such a wiggled filament that can be compared with the response observed in a micro-manipulation experiment . By comparing experimental results with theory, structural properties of the filament were inferred.

Wiggins and Goldstein consider the flexible one armed swimmer in the limit of small amplitude motion. With the assumption of inertialess dynamics, one can derive the following equation    for the function $`h(s,t)`$, describing the displacement of the filament from the horizontal axis as a function of time $`t`$ and arclength $`s`$.

$$\xi _{}_th=k\mathrm{\Delta }^2h$$
(1)
Here $`\xi _{}`$ determines the viscous force, treated as simply a transverse viscous drag and $`k`$ the stiffness of the filament. This “hyperdiffusion” equation has to be solved subject to appropriate boundary conditions (corresponding to different forms of external driving). Simple active driving mechanisms could be an oscillating constraint on the end position, or an oscillatory torque applied at one extremity. The former can be regarded as the simplest example of what has been called *elastohydrodynamics* as it involves the balance of viscous and elastic forces. The latter is a more plausible biological mechanism as it involves no net external force. Both these mechanisms are considered by Wiggins and Goldstein and we also consider both here.

To summarize the predictions of the theory, for a given amplitude of driving the remaining parameters can be grouped together to define a dimensionless “Sperm number”,

$$Sp=\left(\frac{l^4\omega \xi ^{}}{k}\right)^{\frac{1}{4}}$$
(2)
where $`l`$ is the length of the filament and $`\omega `$ the wiggling frequency.

This characterizes the relative magnitudes of the viscous and bending forces. A low value implies that bending forces dominate, a high value viscous forces. As a function of the Sperm number, the theory predicts

* The swimming speed and efficiency (defined as the amount of energy consumed, relative to the amount of energy required to simply drag a passive filament through the fluid at the same velocity), go to zero as Sp goes to zero. This is the stiff limit where the motion is reversible and the scallop theorem applies
* At a sperm number Sp$`4`$ there is a maximum in the both the swimming speed and efficiency (although not at exactly the same value)
* At high sperm numbers a plateau region where the speed and efficiency become independent of Sp, albeit at values lower than the peak.

In this paper we describe a numerical model that allows us to simulate such a driven filament. With the model we can calculate the dynamics, free from restrictions such as small amplitude motion, and with greater scope to specify the type of active forces driving the motion and the boundary conditions applicable for a given physical situation. With such a model, we can test theoretical predictions and also study more complex problems where no analytic solution is available. Here we do both. Looking at small amplitude motion we compare with the theory. Moving on to large amplitude motion we establish to what extent the small amplitude approximation limits the theory.

## 2 Model and Simulation

Our model solves the equations of motion of a discretized elastic filament immersed in a low Reynolds number fluid. Any form of internal and external forcing can be imposed but we restrict ourselves here to an active force, acting on one extremity, that is periodic in time. The hydrodynamics is kept to the approximation of slender body flow , where the local velocity-force relation is reduced to a simple expression in terms of friction coefficients that are shape and position independent. They do nonetheless reflect the difference between friction transverse and longitudinal to the filament. For the problems we are concerned with here the planar driving forces produce planar motions. The model would apply equally well were this not to be the case.

Considering a continuous description of the filament in space and time, one can specify this at any given instant $`t`$ by a curve $`𝐫(s,t)`$, giving a point in space for any value of the arclength parameter $`s`$ (figure 2).

To describe the dynamics we need the local forces acting on the filament. The latter are related to the energy of the model system. Specifically, we have

* A bending elasticity, described by the Hamiltonian

$$H=\frac{1}{2}k_0^LC(s)^2𝑑s$$
(3)
where $`C(s)=\left(\frac{^2𝐫}{s^2}\right)^2`$ is the local curvature and $`k`$ the stiffness.
* A constraint of inextensibility, which can be expressed in terms of the tangent vector as

$$\left|\frac{𝐫}{s}\right|=1$$

and imposes the condition that the filament is, to a first approximation, inextensible.
* An over-damped (negligible mass) equation of motion, which can be written as

$$_t𝐫(s,t)=\left(\frac{1}{\xi _{}}\widehat{𝐭}\widehat{𝐭}+\frac{1}{\xi _{}}\widehat{𝐧}\widehat{𝐧}\right)\frac{\delta H}{\delta 𝐫(s)}$$
(4)
Here, following slender-body theory, the effect of the surrounding fluid is taken as a drag force that is proportional and opposite to the local velocity. This is anisotropic due to the elongated shape of the filament. This requires the presence of a longitudinal drag coefficient $`\xi _{}`$ associated with the projector $`\widehat{𝐭}\widehat{𝐭}`$ along the tangent vector $`𝐭`$, together with a transverse coefficient $`\xi _{}`$ acting along the normal vector $`𝐧`$.

Accordingly, one obtains two equations of motion, one for the evolution of the filament shape, and the other for the tension force $`\tau (s,t)`$, that enforces locally the inextensibility constraint. Expressing the curve shape as the angle $`\psi (s,t)`$ that its local tangent forms with a fixed $`\widehat{𝐱}`$ direction, one can write these equations as (see ):

$$_t\psi =\frac{1}{\xi _{}}\left(k_s^4\psi +\tau _s^2\psi +_s\psi _s\tau \right)+\frac{1}{\xi _{}}_s\psi \left(k_s^2\psi _s\psi +_s\tau \right)$$
(5)
and

$$_s^2\tau \frac{\xi _{}}{\xi _{}}(_s\psi )^2\tau =k_s(_s\psi _s^2\psi )+\frac{\xi _{}}{\xi _{}}_s\psi (k_s^3\psi )$$
(6)
The two nonlinear equations above have then to be solved subject to appropriate boundary conditions. For example, no external forces and torques for a free tail. For the wiggling problems we examine here, the non-equilibrium drive (oscillating end position or torque) is, in these terms, simply a time-dependent boundary condition. Through a functional expansion about the obvious solution for zero drive $`\psi _0(s,t)=0`$, $`\tau _0(s,t)=0`$,

$$\begin{array}{cc}\psi & =ϵ\psi _1+ϵ^2\psi _2+\mathrm{}\\ \tau & =ϵ\tau _1+ϵ^2\tau _2+\mathrm{}\end{array}$$
(7)
it is straightforward to obtain, to second order in $`ϵ`$ the (decoupled) equations

$$_t\psi _1=\frac{k}{\xi _{}}_s^4\psi _1$$

for $`\psi _1`$ and

$$_t^2\tau _2=_s(k_s\psi _1_s^2\psi _1)\frac{\xi _{}}{\xi _{}}(k_s\psi _1_s^3\psi _1)$$

for the tension. Furthermore, expressing the shape of the filament in terms of the transverse and longitudinal “absolute” displacements $`u(s,t)`$ and $`h(s,t)`$ from the direction $`\widehat{𝐱}`$ of the filament’s resting position, one gets to equation 1 for the time evolution of $`h`$ to second order in $`ϵ`$.

In the simulations we use a particle model to solve equations 5 and 6 numerically using an approach similar to molecular dynamics. Time is discretized and the filament is described as a set of $`n`$ point particles rigidly connected by $`n1`$ “links”. The interaction between the particles is constructed so as to reproduce the appropriate collective behavior. For convenience in implementing the algorithm, we do not simulate the over-damped motion, given by equation 4 directly. This would correspond to the zero mass case. Rather, we solve the damped Newton equation for an object with “small” total mass $`m`$. By making the mass small enough we can reproduce the required inertia-less mass independent behavior .

The bending forces acting on the individual particles are defined as follows. If we consider three consecutive discretization points, their positions will lie on one unique circle of radius, $`R_i`$,

$$C_i^2=\left(\frac{1}{R_i}\right)^2=\frac{2}{\mathrm{\Delta }l^2}\left(1\mathrm{cos}(\theta _i)\right)$$

where $`\mathrm{\Delta }l=L/(n1)`$ is the link length and $`\theta _i`$ the angle between two links at the position of bead $`i`$. We introduce a bending potential $`U_i`$ of the form

$$U_i=A\left(1\mathrm{cos}(\theta _i)\right)$$

so that the total bending energy will be

$$H=\underset{i=2}{\overset{n}{}}U_i=A\underset{i=2}{\overset{n}{}}\left(1\mathrm{cos}(\theta _i)\right)$$

which we can compare with a discretization of the integral in 3

$$H=\frac{k}{2}\mathrm{\Delta }l\underset{i=2}{\overset{n}{}}\left(C_i^2\right)$$

This leads to the identification $`A=\frac{k}{\mathrm{\Delta }l}`$. A more sophisticated approach , where the problem is mapped onto the worm-like chain model of Kratky and Porod , leads to a slightly different expression, $`A=\frac{n1/2}{L}k`$. The two expressions are equivalent in the limit $`n\mathrm{}`$, where they reproduce the bending energy of the continuous filament, but they differ for the finite number of beads used in the model. The latter leads to faster convergence in the results as the number of discretization point particles is increased. We therefore chose to adopt it.

The inextensibility constraint is implemented by introducing equal and opposite forces along the links between particles. The magnitude of the forces is computed by imposing a fixed distance $`\mathrm{\Delta }l`$ between consecutive beads at each time step. This is a straightforward matter from the computational point of view, as it involves only the inversion of a tridiagonal matrix .

The viscous drag forces acting on the particles of the model filament are taken as $`F_{ij}=\delta _{ij}(\xi _{}\widehat{𝐭}\widehat{𝐭}+\xi _{}\widehat{𝐧}\widehat{𝐧})𝐯`$, where $`\widehat{𝐭}`$ and $`\widehat{𝐧}`$ are respectively unit vectors parallel and normal to the filament, $`\xi _{},\xi _{}`$ are the longitudinal and transverse friction coefficients, and $`𝐯`$ is the local velocity. This means that hydrodynamics is approximated as a local effect on the filament so the hydrodynamic interaction between different points along the curve does not vary. The global shape of the curve enters only through the anisotropy of the viscous drag coefficients acting on individual points. The ratio of the two coefficients depends on the geometric details of the filament analyzed. For cilia, flagella, or cytoskeletal filaments, its value is typically taken between $`1.4`$ and $`2`$ . We chose to adopt an arbitrary $`1.5`$ in most of our simulations, but we also explored different values, including the cases where the two drags are equal or their ratio is lower than one. The time evolution is evaluated in a molecular dynamics-like fashion, with the only slight subtlety that the Verlet algorithm has to be modified to allow for the velocity dependent anisotropic viscous force . Finally, the active drive at the head is simply implemented as a constraint on the first or first two particles. That is, $`y_1=h_0\mathrm{cos}(\omega t)`$ for the oscillating constraint or a periodic torque, $`T_x=B\mathrm{sin}\left(\omega t\right)`$ realized as a couple of forces applied to the first two beads. Here $`\omega `$ is the driving frequency.

## 3 Results for small deviations

### 3.1 Wave Patterns

Using the ”Sperm Number” Sp defined in section 1, we can characterize the relative magnitude of the viscous and bending forces. To recapitulate, a low value of Sp indicates that bending forces dominate, whereas for low values the dominant forces are viscous. One reason for defining this number comes from the solution of equation 1  . In fact, Sp can be interpreted as a rescaled filament length, where the rescaling factor is a characteristic length $`l_{}=\left(\frac{k}{\omega \xi _{}}\right)^{1/4}`$ that can be used to non-dimensionalize the equation. Both for the oscillating constraint and oscillating torque we recover the fact that the dynamic response, for a fixed driving amplitude, is solely dependent on Sp.

In figure 4 we have plotted the wave patterns for the filament at different values of Sp. These results were obtained using the oscillating constraint. That is, the transverse position at the wiggled end is forced to be sinusoidal in time. The amplitude of the motion is small, the maximum displacement being $`1\%`$ of the filament length. The pictures can be interpreted as “stroboscopic snapshots” of the filament’s motion. For small Sp, bending forces dominate and the stiff filament pivots around a fixed point. This motion is virtually symmetric with respect to time inversions (“reciprocal”). As Sp increases a (damped) wave travels along the filament and time reciprocity is broken. For increasing values of Sp, viscous forces overcome elastic forces and the characteristic length scale of damping of the traveling wave becomes smaller. This requires that the spacing between the beads in our discrete model must also be reduced to give a fixed degree of accuracy. The number of beads in the model (or equivalently the inverse bead spacing) were thus increased with increasing Sp to ensure that the results are within a percent of the true, continuum, values. The oscillating torque gives qualitatively similar results.

All these results are in agreement with the analytical findings of Wiggins and Goldstein  in the small amplitude approximation. The agreement is also quantitative.

### 3.2 Swimming

From the simulation we are also able to compute the velocity and efficiency of the movement generated transverse to the wiggling direction (due to the propulsive force generated by the presence of the active force) as a function of Sp. We define swimming of the immersed object as the generation of motion, through modifications of shape, in the direction along which no external force acts. Both the speed and efficiency, as Wiggins and Goldstein predict, display an optimum value at intermediate (but different) values of Sp. Subsequently they reach a plateau as viscous forces begin to dominate (Sp increases).

According to the “scallop theorem” of low Reynolds number hydrodynamics, reciprocal (time reversion invariant) motion generates no swimming . This is a consequence of the time-reversibility of Stokes flow and sets an important condition for the ability of microorganisms to swim. In our case, this implies that we expect no swimming as Sp approaches zero and the motion approaches reciprocity. This is confirmed by the result in figure 5. The optimum of the velocity is thus the result of a trade-off between non-reciprocity of the motion and damping of the traveling wave.

At this point we should also be able to compare our results quantitatively with those obtained analytically using the approximation of small deviations. However, in this respect the theoretical analysis is somewhat misleading. Computing the time average of the force, as in  and , yields the expression

$$\overline{F}=\frac{h_0^2\xi _{}\omega }{4\sqrt{2}}\text{Y(Sp)}$$

where Y(Sp) is a scaling function that can be computed exactly (figure 6). This expression depends only on the transverse friction coefficient and does not reduce to zero when $`\xi _{}=\xi _{}`$. As such, it is impossible to relate this to the swimming speed. This follows from the fact that if the condition $`\xi _{}=\xi _{}`$ is satisfied there can be no swimming. It is easy to show this must be the case (basically as a consequence of Newton’s third law). The main reason is that, if one considers one particle (i.e. a short piece of filament), the effective viscous drag that it experiences at any moment in time is decoupled from the local configuration of the filament if there is no anisotropy in the friction coefficients. Averaged over one cycle, this always leads, effectively, to reciprocal motion. All the forces sum to zero so there can be no displacement. This is shown more formally in the Appendix. Furthermore, our simulations do indeed yield no average velocity if the two friction coefficients are equal (we use this to check that there is no “numerical” swimming, due to the accumulated errors in the simulation). Thus the result given in equation 3.2, whilst analytically exact, is misleading (probably due to subtleties in the formalism of the over-damped equation of motion). To correct for this anomaly we used the theory and computed instead, following the procedure outlined in , the time average of the swimming velocity given the analytical solution for the shape . This yields (see Appendix)

$$v=h_0^2\left(1\frac{\xi _{}}{\xi _{}}\right)\frac{\omega }{4\sqrt{2}L}\text{Y(Sp)}$$
(8)
where Y(Sp) is again the scaling function specified by Wiggins and Goldstein in computing the average force (figure 6). Note that this expression (equation 8) predicts no swimming when

* Sp $`=0`$ and the motion is reciprocal in time (see fig. 2)
* When the two drag coefficients $`\xi _{}`$ and $`\xi _{}`$ are equal.

consistent with both the scallop theorem and Newton’s third Law.

It also predicts a change in the swimming direction if the friction coefficients are interchanged. Curiously, this reversal of direction has a biological analogue in the organism *Ochromonas* which has a flagellum decorated by lateral projections (mastigonemes) and swims in the same direction as that of the propagating wave. The body follows the flagellum, instead of preceding it as in sperm cells ( , p.11).

Comparing the modified analytical expression for the swimming speed with the simulations, the essential features predicted are obviously present. Both approach zero as $`\text{Sp}^4`$ for small Sperm numbers, but with increasing Sperm number display a maximum. In fact, a careful analysis shows that the agreement, in the small amplitude limit, is exact. The presence of a plateau at high Sp is hard to interpret, in the sense that it predicts velocities for even the “infinitely floppy” filament, where the wave pattern is completely damped in an infinitely small region close to the driven extremity. However, in our simulations we see the velocity dropping only when the size of this damping region is comparable to the distance between two subsequent discretization points, so we have to confirm the analytical result and explain this oddity, as we will see, as a feature of the small deviation approximation.

## 4 Large Angular Deviations

Our simulation contains the full nonlinear model for the dynamics of the filament, its only limitation being the discretization of space and time. Therefore, it is interesting to use it to investigate the limitations of the analytical model when the motion involves shapes that deviate significantly from straight. This is also closer to a real experimental (or biological) situation. The shapes we found often cannot be described by a function, as the displacement from the horizontal axis is not single-valued. This can be observed in figure 7, where we show an example of wave pattern for the case of oscillating large amplitude constraint. In this case, the maximum transverse displacement is $`60\%`$ of the tail length. Looking at this figure, it is obvious that the behavior predicted by equation 1 will be substantially modified.

The first notable area of disagreement is at high values of the Sperm number ($`\text{Sp}>>1`$) where we no longer find a plateau but a slow and steady drop in both speed and efficiency (figure 8). This effect is clearly a consequence of the non-negligible amplitude of the motion because for smaller amplitudes a plateau is indeed reached. This is a limitation of the theory, one respect in which large amplitude motions differ from the small amplitude limit. Further the results for a dimensionless amplitude of $`0.25`$ display a transient plateau that subsequently decays to zero. This implies that for any finite amplitude the dimensionless swimming speed always goes zero for large enough Sp. The smaller the amplitude the longer the plataeu persists, but only for negligible amplitude, is it the asymptotic behaviour. It should be noted that figures 5 and 8 should be interpreted with care. The swimming velocity is plotted in units of the fraction of the length per cylce. To obtain absolute swimming speeds, for a tail of given length and stiffness, we would need to multiply this dimensioless swimming speed by the frequency. The frequency itself is proportional to Sp<sup>4</sup> so a plateau in these plots still implies a swimming speed increasing proportionally with $`\omega `$. The drop from the plateau means that the actual swimming speed will increase with frequency, but a a slower rate. Thus, in practice the one armed swimmer can go as fast as he or she likes by wiggling fast enough.

Secondly, we find that the dependence of the optimum swimming speed, equation 8 predicts as the square of the amplitude of the oscillating constraint, becomes linear for higher amplitude oscillations (figure 9).

Thus far, we have not been able to show why this is the case, but we believe it is related to the following. For small amplitude motion the elastic wave simply propagates along an essentially straight filament. As the amplitude increases, this is no longer true because the filament itself is significantly bent and, so far as the damping is concerned, it is the distance along the filament that is relevant. This is no longer the same quantity as the absolute distance. This seems to lead to an increase in the effective length of the filament.

From these results it is clear that the amplitude of the drive for which the small deviation approximation breaks down depends on the value of Sp, being greater for smaller Sperm Numbers. At the optimal value for the speed, Sp$`=4`$, the approximation holds for maximum transverse displacements of up to $`20\%`$ of the tail length, which is well beyond the point one would expect the assumptions to be valid. However, for a realistic experiment with actin or a microtubule, the values of Sp are much higher than $`4`$, and the value for the theshold is much lower. For example, for an actin filament of $`50\mu m`$, driven at 1 cycle/second at an amplitude of $`25\%`$ of its length, we estimated a speed of about $`1,7\mu m/sec`$ with the small deviation model, whereas our simulation predicts a reduction of this value by a factor $`1/4`$.

For external driving in the form of a torque applied at one end, we have only considered large amplitude motion. Specifically, the pre-factor $`B`$ was adjusted to produce a maximum angle at the driven end of $`60^o`$. This clearly violates the small angle approximation of Wiggins and Goldstein  but is more consistent with the head deflections found in practice for swimming organisms. In figure 10 we plot the efficiency and the mean velocity as functions of Sp. Once again, the two curves agree qualitatively with those found analytically by Wiggins and Goldstein in that there is a peak speed and efficiency. The values are at slightly different values of Sp and, because the small angle approximation is violated, not quantitatively predicted by equation 8.

We should add here a few comments. Notably, the peak efficiency of less than $`1\%`$ seems very low. However, this depends strongly on the amplitude of the motion. If we go above the $`60^o`$ limit we have imposed here for the torque, or to driving amplitudes of greater than half the length of the filament, it is possible to reach values of $`2\%`$ before the motion becomes unstable. This is similar to the efficiency typical for both the helical screw mechanism used by bacteria and the sperm motion . Thus, the one-armed swimmer operating at peak efficiency is a plausible and not especially inefficient entity. Note also that the efficiency (which is dimensionless), as well as the swimming speed, decays to zero rather than reaching a plateau value. This means that in absolute terms the one-armed swimmer can carry on with increasing speed by increasing its wiggling frequency but only at the price of decreasing efficiency.

## 5 Conclusions

We have described a simulation method that can be used to study the motion of driven elastic filaments in a low Reynolds number flow. Here, the hydrodynamic friction is treated quite simply, consistent with comparing with analytically tractable theories. A more complete calculation of the hydrodynamic effects would no doubt be instructive. In particular, the friction coefficients are not, as we assume here, independent of distance along the filament. At the expense of a little more computational complexity, such effects could be incorporated into our model in a straightforward manner. We showed that, within this approximation, the picture suggested by Wiggins and Goldstein for the linear regime of small angular deviations from the straight position is essentially correct. Our results for the motion of the filament, show good agreement with their analytical calculations. There is an optimal balance between bending forces and viscous forces that leads to a maximum propulsive speed and efficiency. However, in one quantitative respect our results suggest that their analysis is limited. We could not relate their expression for the average force exerted by the filament to the swimming speed. Instead, we used their model to compute an expression for the average swimming velocity that is physically more plausible and agrees with the simulation results.

For large amplitude motion, we found that the dependence of the swimming speed on both sperm number and amplitude was significantly modified relative to the small amplitude case. Further, we postulate that this is due to the fact that in a highly distorted filament the wave travels along a notably different path than is the case for small amplitude motion. A quantitative understanding of this effect is still, however, lacking. Nonetheless, the general picture derived from the linear theory, of an optimal compromise between the bending required to break time reversibility and excessive damping suppressing motion along the filament, remains valid. The most significant difference we found was that there will come a point beyond which increasing the wiggling frequency leads to a drop in efficiency. The theory, on the other hand, predicts that the efficiency remains constant.

For realistic amplitudes of oscillation, we found that one-armed swimming is, speed and efficiency-wise, a plausible strategy a microorganism might use to get around. It is also a sight simpler than the helical screw mechanism used by most bacteria. This requires a rotary joint . Nonetheless, while we stand open to correction, we have not been able to identify a single organism that actually adopts this strategy. Perhaps the most interesting question surrounding the one armed swimmer is thus: why doesn’t it exist? Based on our results, we suggest two hypotheses. First, localized bending of the tail requires implausibly high energy densities. Second, the existence of an evolutionary barrier. It is useless trying to swim with a short or slow moving tail. Note that at small Sp (that is, low frequency and or a short appendage) there is nothing to be gained in terms of motility. This is not the case for either the helical screw mechanism, commonly used by bacteria, or the traveling wave, used by spermatozoa. Both of these give a maximum swimming speed and efficiency at low sperm numbers.

Regarding experimental studies of in vitro motility, one main problem so far is that the force involved has been too small to be detectable with an optical trapping experiment. This limitation could be resolved simply by time, as it is reasonable to expect that the resolution of experiments will improve. On the other hand, by means of the model one could try to find the region of the parameter space where this force is expected to be highest, and try to design an “optimal experiment” where the motility could actually be quantified.

We would like to thank Catalin Tanase, Marileen Dogterom and Daan Frenkel for discussion and help.

## Appendix

### 5.1 Unequal friction coefficients as a condition for motility

It is possible to show that there can be no movement if the viscous drag coefficients are equal. Most conveniently, we work with the the discrete model. Since the discrete model produces the continuum result in the limit that the number of beads $`N`$ goes to infinity, there is no loss of generality (so long as the answer does not depend on $`N`$). The equation of motion for the center of mass is

$$m\frac{d\overline{𝐯}}{dt}=\underset{i=1}{\overset{N}{}}𝐟_i$$
(9)
where $`\overline{𝐯}=\frac{1}{N}_{i=1}^N𝐯_i`$ is the center of mass velocity, $`m`$ is the total mass. The total force $`𝐟_i`$ on each bead consists of a bending force $`𝐟_{iB}`$ a tension force $`𝐟_{iT}`$ an hydrodynamic force $`𝐟_{iH}`$ and an “external” force $`𝐟_{iX}`$ which accounts for the external drive. We know that by definition

$$\underset{i=1}{\overset{N}{}}𝐟_{iB}=\underset{i=1}{\overset{N}{}}𝐟_{iT}=0$$

Now, the external periodic force is applied only at one extremity, so that $`_{i=1}^n𝐟_{iX}=𝐟_1(t)`$ and equation 9 can be written as

$$m\frac{d\overline{𝐯}}{dt}=\underset{i=1}{\overset{N}{}}𝐟_{iH}+𝐟_1(t)$$

Integrating on a cycle we get

$$m\mathrm{\Delta }\overline{𝐯}_{\text{cycle}}=_{\text{cycle}}\text{dt }\underset{i=1}{\overset{N}{}}𝐟_{iH}$$

The hydrodynamic force on bead $`i`$ is written in the form $`𝐟_{iH}=\xi _i𝐯_i\xi _i\overline{𝐯}_i`$ (with $`\xi _{i\times }=\xi _\times /N`$). Thus, the effective drag on one particle depends on the local configuration of the filament shape.

If the two friction coefficients are the same $`\xi _i=\xi _i=\xi _i`$ then

$$m\mathrm{\Delta }\overline{𝐯}_{\text{cycle}}=\xi _{\text{cycle}}\text{dt }\overline{𝐯}$$
(10)
which necessarily leads to zero (or decaying to zero) global velocity. On the other hand, if the two drags are different, the right hand side integral in equation 10 can be written as

$$_{\text{cycle}}\text{dt }\xi _{}(t)\overline{𝐯}(t)$$

where $`\xi _{}`$ is an effective drag which depends on time through the configuration of the filament. This integral in general gives a number once, in the spirit of resistive force theory, the configuration is plugged in, and swimming is not, therefore, precluded.

### 5.2 Analytical computation of the mean swimming velocity in the the small angular deviation approximation

Her we outline the procedure adopted to calculate analytically the average of the swimming speed using the small deviation approximation. This calculation largely follows the methodology used in  on a different model.

We can define the time average of the swimming speed (projected along its only nonzero component along the $`\widehat{𝐱}`$ direction) as

$$v=\underset{t\mathrm{}}{lim}\frac{1}{t}_0^L\text{dt }_t𝐫\widehat{𝐱}$$

The expression for $`_t𝐫`$ can be obtained from equation 4 in terms of the local angle $`\psi `$ as

$$_t𝐫=\frac{1}{\xi _{}}𝐧(k_s^3\psi +\tau _s\psi )+\frac{1}{\xi _{}}𝐭(k_s^2\psi _s\psi +_s\tau )$$

Fixing a reference frame one can consider the “comoving” frame with respect to the filament, and expand $`\psi `$ and $`\tau `$, together with the absolute displacements $`h`$ and $`u`$, and the swimming speed $`v`$ as in formula 7. Following this reasoning we can rewrite in vector notation the formula above for $`s=0`$ as

$$\begin{array}{c}(ϵv_1+ϵ^2v_2,0)+_t(ϵu_1(0)+ϵ^2u_2(0),ϵh_1(0)+ϵ^2h_2(0))=\\ \\ =\left[\frac{1}{\xi _{}}ϵ^2k\psi _1_s^3\psi _1+\frac{1}{\xi _{}}(kϵ^2_s^2\psi _1_s\psi _1+ϵ^2_s\tau _2)\right]_0\widehat{𝐱}+\left[\frac{1}{\xi _{}}(kϵ_s^3\psi _1)\right]_0\widehat{𝐲}\end{array}$$
(11)
where we stop the expansion to second order in $`ϵ`$. Expressing the equality for the different powers of $`ϵ`$ one gets $`v_1=0`$ and

$$v_2+_tu_2(0)=\psi _1(0)_th_1(0)+\frac{1}{\xi _{}}(k_s\psi _1(0)_s^2\psi _10+_s\tau _2(0))$$

$`\tau _2(0)`$ is obtained integrating equation 11. Taking into account the boundary conditions for $`h`$, $`u`$ and $`\tau `$ one gets to the expression

$$v_2+_tu_2(0)=\frac{\xi _{}\xi _{}}{\xi _{}L}_0^L\text{ds }_sh_1_th_1+\frac{1}{L}_0^L\text{ds }_0^s\text{ds’}\frac{1}{2}_t(_sh_1)^2$$

which, plugging in the analytical solution of (1) and time averaging, gives (8).