Nonlinear parametric excitation effect induces stability transitions in swimming direction of flexible superparamagnetic microswimmers
NNonlinear parametric excitation effect induces stability transitions in swimming direction of flexible superparamagnetic microswimmers
Yuval Harduf , Dongdong Jin , Yizhar Or and Li Zhang Abstract:
Microscopic artificial swimmers have recently become highly attractive due to their promising potential for biomedical applications. The pioneering work of Dreyfus et al (2005) have demonstrated the motion of a microswimmer with an undulating chain of superparamagnetic beads, which is actuated by applying an oscillating external magnetic field. Interestingly, they have shown that the swimmer's orientation undergoes a 90 o -transition when the magnetic field's oscillations amplitude is increased above a critical value. In this work, we further investigate this transitions both theoretically and experimentally by using numerical simulations and presenting a novel flexible microswimmer with a superparamagnetic head. We prove that this effect depends on both frequency and amplitude of the oscillating magnetic field, and demonstrate existence of an optimal amplitude achieving maximal swimming speed. Asymptotic analysis of a minimal two-link model reveals that the changes in the swimmer's direction represent stability transitions which are induced by a nonlinear parametric excitation. Faculty of Mechanical Engineering, Technion – Israel Institute of Technology, Haifa 3200003, Israel. Department of Mechanical and Automation Engineering, The Chinese University of Hong Kong, Shatin NT, Hong Kong SAR, China. Chow Yuk Ho Technology Center for Innovative Medicine, The Chinese University of Hong Kong, Shatin NT, Hong Kong SAR, China. † These authors contributed equally to this work. * Corresponding author, email: [email protected] and [email protected] ntroduction
Remotely controlled micro-robotic swimmers have recently become highly attractive in the engineering science community, due to their promising potential for biomedicine . Drawing biological inspiration from swimming microorganisms , micro-robotic swimmers are planned to be utilized for in-body biomedical tasks such as targeted drug delivery , therapeutic diagnosis, and even assisting sperm's motility . The motion of microswimmers is dominated by viscous drag forces while inertial effects are negligible . A common actuation method of microswimmers is applying time-varying external magnetic fields . Several works presented rigid helical microswimmers actuated by a rotating magnetic field which induces corkscrew-like swimming. Dreyfus et al presented a microswimmer composed of a flexible chain of superparamagnetic beads connected to a red blood cell, where a spatially uniform magnetic field B varies in time as B ( t ) =B x x ̂ +B y sin ( πft ) y ̂ = B x ( x ̂ + sin( ωt ) y ̂ ). (1) In a later work , it has been shown that when the amplitude ratio = B y / B x is increased beyond a critical value of √2 , microswimmers composed of a chain of beads reached a sharp 90 o -transition in their mean orientation. The "ineffective" swimmers studied in did not display any net swimming, since they were not attached to a load that is needed for breaking the swimmer's front-back symmetry. The 90 o -transition has also been confirmed by numerical simulations for a microswimmer with a load, which also displayed a change in the swimming direction from x to y . Nevertheless, no physical explanation has been provided to this effect. In this work, we further study this transition both theoretically and experimentally. Unlike previous works , we find that the conditions for transition can depend on both the ratio and the frequency f . For a fixed frequency, we also find an optimal value of which maximizes the swimming speed in y direction. These results are confirmed by using numerical simulations of a multi-link superparamagnetic microswimmer, and by presenting experimental results of a novel elastic microswimmer with a superparamagnetic head that breaks front-back symmetry. Next, we provide an analytic explanation to this transition by presenting a minimal theoretical model of a two-link swimmer with an elastic joint. We explicitly formulate the swimmer's dynamics as a low-dimensional nonlinear system with parametric excitation, which induces transitions in the stability of different periodic solutions, corresponding to oscillations about and swimming along x or y directions. In the limit of fast oscillations and low stiffness, asymptotic analysis using the method of multiple scales reveals the previously observed transition at β =√ . Approximate expressions for the mean swimming speed are obtained, which clearly indicate existence of an optimal value of β that attains maximal speed. In the limit of fast oscillations and high stiffness, asymptotic approximation of the system yields a scalar second-order differential equation with parametric excitation , which resembles the well-known Kapitza pendulum . Analysis of stability transitions in this system gives conditions that depend on both the ratio β and the frequency f , in agreement with our experiments and simulation results. umerical simulations We now present our numerical simulations of a multi-link microswimmer. A schematic description of the swimmer appears in Fig. 1a. It consists of a chain of
N=10 slender rigid links connected by rotary joints, where each link carries a superparamagnetic bead. The chain, whose total length is L , and is attached to a spherical head which represents the load. The external magnetic field undergoes planar oscillations as described in (1). A detailed derivation of the microswimmer's equations of motion appears in the Supplementary Information (SI), and key points are briefly summarized here. Elasticity of the chain is represented by torsion springs, such that the elastic bending torque acting on each joint is given by τ i =-kϕ i , where ϕ i is the relative angle at joint i . The magnetic field induces an external torque L i on each link, which is dictated by the interaction between dipoles generated by the superparamagnetic beads. For calculating viscous drag forces acting on the swimmer, we use Stokes drag for a sphere and resistive force theory for the slender links, while hydrodynamic interaction is neglected. All this gives a nonlinear system of first-order differential equations which are integrated numerically (see SI). We calibrate our system parameters using data from the experiments of Dreyfus et al , and then simulate the microswimmer's motion under a moderate amplitude ratio of β≈1 for different values of actuation frequency f of the magnetic field as given in (1). We find that the swimmer's net motion is along x direction, and that a maximal speed is obtained at an optimal value of the frequency. The results of these simulations show good agreement with those reported by Dreyfus et al , as shown in Fig. 1b. Next, we verify the results given by Roper et al , by considering an elastic chain without a load, composing an "ineffective swimmer". For a fixed actuation frequency f of the magnetic field in (1) and a fixed magnitude of the constant component of magnetic field B x , we gradually increase the amplitude of the oscillating component B y , until reaching a critical value where the swimmer's orientation converges in steady state to oscillations about y rather than x direction. For a frequency of f=50Hz , the critical value of B y as a function of B x is shown in Fig. 1c, showing a nearly constant critical ratio of β cr =√2 , in agreement with Roper et al . However, when the actuation frequency is decreased to f=5Hz , a noticeable change occurs. As shown in Fig. 1c, the value of β cr now increases significantly above √2 for large values of B x . This proves that the conditions for transition in the swimmer's direction depend on frequency f and the constant magnetic field B x , in addition to the amplitude ratio 𝛽 . In order to show a change in the swimming direction and not just mean orientation of the swimmer's oscillations, we conducted additional simulations for Dreyfus' microswimmer with a spherical load. Fig. 1d shows the mean swimming speed V as a function of the actuation frequency f for a constant amplitude ratio of β=1.6, clearly showing a transition from swimming in x direction to y direction at f ≈1.8Hz . This proves again that this transition strongly depends on the actuation frequency f , in addition to β . Finally, we simulate the microswimmer's motion under a constant frequency f=20Hz , and plot the speed V as a function of β , in Fig. 1e. The results show that the speed depends non-monotonically on β , and that an optimal value of β exists, which maximizes the swimming speed. Remarkably, this maximal speed is obtained beyond the transition to swimming in y direction. -2 -1 f [ Hz ] s p ee d V [ m = s e c ] Figure 1: (a) The planar multi-link microswimmer model with a non-magnetic spherical load. The rigid slender links carry small superparamagnetic beads (dashed circles), and the connecting joints contain torsion springs. (b) Normalized swimming speed V x / L as a function of actuation frequency (logarithmic scale). Comparison of our numerical simulations (solid lines) with Dreyfus' experimental data (discrete markers). Parameter values are given in the SI. (c) Critical value of the magnetic field's oscillating component B y = B x for transition from oscillations about x direction to y direction, as a function of the constant component B x . The 'x' markers are from the experiment of Roper et al under frequency f =50Hz. The 'o' markers are our simulations of a multi-link swimmer without a spherical load for f =50Hz. The '+' markers are our simulations of a multi-link swimmer without a spherical load for f =5Hz. The dashed line is of slope √2 . It can be seen that the critical value of B y / B x depends also on f and B x , and differs from √2 . (d) Numerical simulations of the mean speed V of Dreyfus' microswimmer with a spherical load as a function of actuation frequency f (logarithmic scale) for B y =1.6B x =13.92mT . The dashed vertical line denotes critical frequency where the steady-state solution changes from swimming in x direction to y direction. The labels V x and V y denote swimming speed in x and y directions, respectively. (e) Numerical simulations of the mean speed V of Dreyfus' microswimmer with a spherical load as a function of amplitude ratio (logarithmic scale) for B x =8.7mT and f=10Hz. The dashed vertical line denotes critical value of ≈ where the steady-state solution changes from swimming in x direction to y direction. It is noticeably larger than 𝛽 = √2 (dash-dotted line). The labels V x and V y denote swimming speed in x and y directions, respectively. It can be seen that there exist an optimal value of that maximizes swimming speed V y , which is greater than V x . a ! [ rad=sec ] V x = ! L x V y V B x ( mT ) c r i t i c a l B y ( m T ) -6 ratio - s p ee d V [ m = s e c ] V x V y - opt : : - transition : : - = p b c d e xperiment We now present experimental results with a newly fabricated microswimmer in order to verify our findings on the transition in the swimming direction. Our experimental microswimmer prototype consists of two parts: a superparamagnetic head actuated by an external oscillating magnetic field, and a flexible main body which enables breaking the time-reversibility of motion during body undulations. The main body of the microswimmer is made of polypyrrole (Ppy), a kind of conductive and elastic polymer. The superparamagnetic head was made by embedding Fe O nanoparticles (Fe O NPs) at the end of the Ppy filament. The fabrication process, including preparation of superparamagnetic Fe O NPs, electrophoretic codeposition of Fe O NPs and Ppy followed by electrochemical deposition of Ppy in an anodic aluminum oxide (AAO) template, is shown in Fig. 2 (Detailed information regarding all the experimental procedures are available in the SI). It is worth mentioning that the average diameter and zeta potential of Fe O NPs prepared by hydrothermal method in electrophoresis bath was approximately 80 nm and -50 mV, respectively (shown in Fig. S1(b) and (d)), demonstrating that they could be assembled into AAO template with 200 nm-diameter pores by electrophoretic deposition. Typical SEM images (see Fig. S1(e) and (f) in the SI) depict the morphology of as-prepared microswimmers, and it is found that the typical length and diameter of the microswimmers are 25-30 μm and approximately 200 nm, respectively. The microswimmer is actuated by a planar oscillating magnetic field generated using a 3-axis Helmholtz electromagnetic coil setup. The generated planar oscillating field B ( t ) was a spatially uniform time-varying external magnetic field, having a constant component of B x =5mT in x direction and a sinusoidally oscillating component of 𝑏 sin( ft ) in y direction, as described in Eq. (1). Under small-amplitude oscillation of the magnetic field where the ratio 𝛽 is as low as 0.5, the swimming behavior of the microswimmer is very close to previously reported flexible magnetic Figure 2:
Schematic of the fabrication process of the microswimmer and the input planar oscillating magnetic field icroswimmers . For better understanding, Fig. 2(a) shows a sequence of pictures illustrating the dynamic motion of microswimmer during a 20-second duration under an oscillating magnetic field with the input parameters of 𝛽 =0.5 and f =10 Hz. As expected, the superparamagnetic head of the microswimmer continuously tends to align with the direction of the magnetic field, resulting in undulations of the Ppy flexible main body. As a consequence, the microswimmer oscillates about and swims along the mean direction of the magnetic field, x axis, with a translational speed of V x =0.69 μm/s. Repeating this process for varied frequencies f , our experimental results indicate a resonance-like dependence of swimming velocity on f as shown in Fig. 2(b). That is, the maximal speed V x is attained at an optimal input frequency. This resonance-like effect is in agreement with our numerical simulations and the results of Dreyfus et al . Next, when the amplitude ratio 𝛽 is increased to 1, the dynamic motion and average speed of the microswimmer as a function of frequency are shown in Fig. 2(c) and (d). Figure 3:
The microswimmer's motion status under an oscillating magnetic field with different values of β and f. (a) Motion snapshots of the microswimmer for β =0.5, f =10 Hz. (b) Average speed of the microswimmer along x direction as a function of frequency f for β =0.5. (c) Motion snapshots for β =1, f =4 Hz. (d) Average speed of the microswimmer along x direction as a function of frequency f for β =1. (e) (i)-(v) display the 90°-transition in oscillating direction of microswimmer for β =2 while f is increased from 2 to 2.5 Hz. (vi)-(x) show the dynamic motion of microswimmer in y direction when β =2 and f =3 Hz. (f) Average speed in x and y directions as a function of frequency for β =2. In (a)(c)(e), the black arrows represent the oscillating direction of the microswimmer and all the scale bars are 10 μm. Full movies appear in the Supplementary Information (Movie S1, 2, 3). In (b)(d)(f), “ V = V x ” or “ V = V y ” represent that the microswimmer swims along the x or y direction, respectively. The average speeds and standard deviation values are obtained from at least three speeds obtained from an image sequence.
No significant difference compared to the condition when 𝛽 =0.5 is found except for the maximal speed and optimal frequency values. At this time, a maximal speed of 0.75-0.96 μm/s is obtained when the frequency of the magnetic field is 4 Hz. Remarkably, when we further increase 𝛽 to 2, a difference in the motion of the microswimmer is observed. As depicted in the images of (i) and (ii) of Fig. 2(e), the microswimmer still oscillates about the x axis when f is lower than 2 Hz. After a time of 20 seconds, when f is slightly increased to 2.5 Hz, the microswimmer almost immediately changes its oscillating direction from the x axis to the y axis, which is shown in (iii)-(v) of Fig. 2(e). This new motion condition (oscillations about the y axis) is very stable and further increasing f does not have any influence on it. Moreover, except for oscillating direction, 90° transition in the swimming direction of microswimmer also occurred. Images of (vi)-(x) of Fig. 2(e) show that the microswimmer starts swimming along the y direction with a velocity of ~1.34 μm/s (at f =3 Hz), which is comparable (and even larger) to the translational speed in x direction shown in Fig. 2(a) and (c). Finally, Fig. 2(f) shows a plot of the mean swimming speed as a function of frequency. We conclude that the 90°-transition in both oscillating orientation and swimming direction of the microswimmer can be obtained for a fixed value of 𝛽 by increasing the frequency f above a critical value. This finding is in agreement with our numerical simulations above. Next, we conduct a series of experiments in order to find the dependence of the transition conditions on both 𝛽 and f . Under different values of 𝛽 (from 0.5 to 3.5), the critical values of frequency, above which the microswimmer exhibits transition of motion to y direction, are shown in Fig. 3(a). In our experiments, when 𝛽 is lower than 2 the transition does not occur for any frequency, as indicated by the empty left region in Fig. 3(a). It is found that the critical frequency decreases with increasing 𝛽 . For better illustration, the transition of oscillating orientation under varying frequency for 𝛽 =3 is shown in Fig. S4 and Movie S4 in the SI. When the frequency f is held constant while 𝛽 is changed incrementally, we observe that an optimal value of 𝛽 exists, for which maximal swimming speed is achieved. This optimal speed is found to be within the region of “ V = V y ”, where the microswimmer swims along the y direction. For example, the swimming speed as a function of 𝛽 for f =3 Hz is shown in Fig. 3(b). For 0< 𝛽 <2, increasing 𝛽 leads to the increase of velocity in the x direction. For 𝛽 >2, the microswimmer begins to move in the y direction, and it is found that microswimmer exhibits a maximal speed of 1.52-1.69 μm/s when β =2.5. The dynamic motion with those actuation parameters are also shown in Fig. S5 and Movie S5 in the Supplementary Information. The existence of this optimal value of 𝛽 , which has not been observed in previous works , is also in agreement with our numerical simulations. Figure 4: (a) The critical transition frequency f as a function of the amplitude ratio β . Standard deviation values are obtained by repeating the 90°-transition in oscillating orientation three times. (b) The swimming speed of the microswimmer in x and y directions as a function of β for f =3 Hz. The average speeds and standard deviation values are obtained from at least three speeds calculated from an image sequence.
Theoretical two-link model
In order to provide a theoretical explanation to the new findings observed in this microswimmer, we study the simplest possible model that describes our swimmer - two slender rigid links connected by an elastic joint, which move in ( x,y ) plane only (Fig. 5a). The length of each link is l , and their cross-section radius is a . A similar model has been studied previously for a microswimmer with ferromagnetic links , and it has been shown that the external magnetic actuation enables this swimmer to overcome the well-known scallop theorem and generate non-reversible undulations that result in net swimming motion. In our work, the rightmost link (head) is made of a superparamgnetic material with anisotropy of and director t ̂ which is aligned with the link's longitudinal axis. The magnetic field B ( t ) polarizes the link and induces a magnetic moment which, in turn, interacts with the field and generates a torque L given by ˆ ˆ( )( ) v L t B t B , where v is the link's volume. Since the magnetic field B also lies in ( x,y ) plane, the torque it generates on the link simplifies to (see SI): (2) vB L z where B is the magnitude of the field B and is the relative angle between its direction and the link's axis t ̂ (Fig. 5a). This implies that the magnetic torque L vanishes when the directions of the magnetic field and the link are either aligned ( ) or perpendicular (γ=±π/2). It can readily be seen that for a constant magnetic field B , the aligned orientations are stable while the perpendicular ones are unstable. However, when B ( t ) is oscillating as in (1), interaction with effects of swimmer's elasticity and fluid's viscosity changes the stability characteristic in a complicated way, as explained below. The microswimmer's dynamic equations of motion can be explicitly obtained as a system of nonlinear first-order differential equations (see SI for detailed derivation). The system's dynamics can be analyzed asymptotically in different physical limits, as follows. In the limit of small amplitude ratio 𝛽 of the magnetic field in (1), the system's solution can be expanded as power series in 𝛽 using perturbation expansion (see SI). This solution involves stable oscillations about θ=ϕ=0, and swimming in x direction, whereas oscillations about =90 o are unstable. Leading-order expression for the mean swimming speed in x direction is obtained as:
22 42 2 2 2 2 t )64 25 (64 64 1 kkx k m m m ltV t t Ot t (3) where t m , t k are characteristic time scales of the system, given as:
2, where 12 log( )6 tt tm kx l viscous c l viscoust tv magnetic k elastc lic l /B ac (4) and is the fluid's dynamic viscosity. Eq. (3) captures the resonance-like frequency dependence of the swimmer for low values of 𝛽 , similar to that found in our simulations and experiments, and in agreement with Dreyfus et al . Next, we study the transitions that occur when 𝛽 >1, where this small-amplitude approximation is no longer valid. We first consider the limit of very fast oscillations of the magnetic field, i.e. ω -1 ≪ t m ,t k . Using the method of multiple scales , we define fast and slow nondimensional time scales t f = t and t s = t/t m and separate the system's dynamics into fast oscillations superimposed on a slowly evolving average solution. It can then be shown (see SI) that the slow-time average dynamics of the two angles θ ̅ ( t s ) ,ϕ ̅ ( t s ) can be obtained asymptotically as: ss ddtddt where = t m / t k . (5) The system of slow dynamics in (5) has only two equilibrium states, {θ ̅ e =ϕ ̅ e =0} and {θ ̅ e =90 ∘ , ϕ ̅ e =0}, which correspond to oscillations about x and y directions, respectively. Moreover, linearization of (5) about the equilibrium points reveals that θ ̅ e =0 is a stable equilibrium for β< √ , whereas θ ̅ e =90 ∘ becomes stable for β> √ (see SI). This is in agreement with the results obtained in Roper et al . Furthermore, averaging over the fast oscillations gives an approximation of the mean swimming speed of the swimmer in x and y directions as (see SI): x e ey e e lfVV lf (6) where ϵ= ( ωt m ) -1 . This confirms that oscillations about θ ̅ e =0 correspond to motion in x direction, whereas oscillations about θ ̅ e =90 ∘ correspond to motion in y direction. Importantly, (6) also indicates that there exist an optimal value of 𝛽 that maximizes the mean speed, as corroborated in our experimental observations. However, this limit of fast oscillations does not capture the dependence of the critical value of value of 𝛽 on the frequency f . In case where the oscillation frequency f is not so large, the system's equations can be simplified by assuming a limit of stiff spring, t k ≪ t m , which implies small joint angle ( t ). Linearizing by , one obtains a scalar second-order differential equation for ( t ) as (see SI):
32 32 10 cos(2 ) cos(2 ) 2 40 sin(2 ) sin( ) 2 cos(2 ) 10 (2 ) sin(2 ) =4 5 cos( ) sin( ) sicos( ) n2 t tt tt t (7) Equation (7) is a second-order system with parametric excitation , whose structure is remarkably similar to the well-known equation of Kapitza pendulum with tilted base oscillations , see Figure 6a. Using numerical integration of (7), periodic solutions can be computed, and their orbital stability properties can be obtained by evaluating the Floquet multipliers of the linearized variational equation for small perturbations about each periodic solution (see SI). This enables plotting stability regions of the two different solutions in the plane of ( , f ), as shown in the stability transition curves in Figure 6b. These curves show that the stability transitions depend on both amplitude 𝛽 and frequency f of the oscillations (in analogy to Kaptiza pendulum). It can be seen that for large frequency f , the transitions value approaches β→ √ in agreement with the results of Roper and our analysis. Another interesting observation is the existence of a significant region of bistability , where the solutions of swimming in x and y directions are both locally stable, with different regions of attractions in terms of initial conditions, as demonstrated in Figure 5d. Nevertheless, it is often practically difficult to observe co-existence of solutions in a real microswimmer. This is because one cannot precisely enforce any desired initial conditions, while some periodic solutions may have very narrow regions of attraction. Calibrating parameters of the theoretical models by fitting to the experiments gives t m = t k =0.1 sec , and l =25 m . Figure 6c plots a comparison of the stability transition curves -0.4 0 0.4 0.8 1.2 1.6-0.100.1 [rad] ? [ r a d ] Figure 5: (a) The planar two-link microswimmer model. The head link is superparamagnetic. (b-d) Solution trajectories of the two-link microswimmer in the plane of angles (𝜃, 𝜙) , for model parameters t k = t m = , and f =0.16. (b) Solutions under small oscillations = (t) and (t) oscillate about 0 o and the swimmer moves in x direction. (c) Solutions under large oscillations =
10 - convergence to a periodic trajectory where (t) oscillates about 90 o while (t) oscillates about 0 o , and the swimmer moves in y direction. (d) Solutions under intermediate oscillations = =0 o and =90 o are both stable, and convergence is determined by initial conditions. -0.01 0 0.01 0.02 0.03-10-505 x 10 -3 [rad] ? [ r a d ] [rad] ? [ r a d ] a b d [rad] c btained experimentally (squares) with our theoretical analysis, where the solid curves denote the theoretical transition curves of local stability, and the region enclosed between them is a region of theoretical bistability. In addition, the dashed line plotted in Fig. 6c is the curve of transition between swimming in x and y directions, under specific initial conditions of (0)= (0)=0. That is, the swimmer is initially straightened and aligned with x axis. This initial condition is practically realized in our experiments by initializing the microswimmer under a constant magnetic field in x direction. Therefore, this curve gives a meaningful comparison to our experiment, showing a remarkable agreement with the theoretical two-link model. Finally, Fig. 6d shows a comparison of the net swimming speeds V x and V y obtained from the two-link model (solid lines) with our experimental results, showing an excellent agreement. Figure 6: (a) The example of Kapitza pendulum , an inverted pendulum with an oscillating pivot. The pendulum has mass m and length l. Its tilt angle is denoted by 𝜓 and gravity acceleration is g . The pivot of the pendulum is oscillating along a direction with inclination angle 𝜑 , and displacement of 𝑎 ⋅ cos(𝜔𝑡) . The equation of motion of this pendulum is given by ψ ሷ - ቀ gl + aω l cosφ cos ( ωt )ቁ sinψ=- aω l sin φ cos ( ωt ) cosψ . The structure of this equation is remarkably similar to eq. (7), where the pendulum angle ψ is replaced by . Nevertheless, equation (7) also contains an additional "parametric damping term", i.e. proportional to 𝜃ሶ . (b) Plot of stability transition curves in the plane of parameters 𝛽 and f , the amplitude and frequency of the magnetic field's oscillations. The lower curve is stability limit of swimmer's oscillations about =0 o , and the upper curve is stability limit of swimmer's oscillations about =90 o . This plot shows that stability transitions depend on both 𝛽 and f , as observed in our experiments. The region between the two curves is bistability regime, where both solutions are stable under different initial conditions. The dashed curve represents the transition curevs under practical initial conditions of (c) Zoom into the plot of stability transition curves, compared to our experimental results. (d) Plot of net swimming speed as a function of for f =3Hz, compared to our experimental results. The model correctly captures the maximum in V y for optimal value of . The dashed line is the bistable solution of V x which has not been reached under practical initial conditions of p f t r a n s iti on [ H z ] =90°stable =0°stable a b c d n summary, we have studied stability transitions in microswimmers with superparamagnetic links under planar oscillations of an external magnetic field. Unlike previous studies, we have shown that the conditions for transition depend on various parameters, including amplitude and frequency of the oscillations. We used numerical simulations and experimental measurements of a microswimmer with superparamagnetic head. A simple two-link model enabled asymptotic analysis of the swimmer's nonlinear dynamics, and provided an explanation to the stability transitions which are induced by a nonlinear parametric excitation. References: Peyer, K. E., Zhang, L., & Nelson, B. J. (2013). Bio-inspired magnetic swimming microrobots for biomedical applications. Nanoscale, 5(4), 1259-1272. 2.
Sitti, M., Ceylan, H., Hu, W., Giltinan, J., Turan, M., Yim, S., & Diller, E. (2015). Biomedical applications of untethered mobile milli/microrobots. Proceedings of the IEEE, 103(2), 205-224. 3.
Krishnamurthy, D., Katsikis, G., Bhargava, A., & Prakash, M. (2016). Schistosoma mansoni cercariae swim efficiently by exploiting an elastohydrodynamic coupling. Nature Physics. 4.
Lauga, E., & Powers, T. R. (2009). The hydrodynamics of swimming microorganisms. Reports on Progress in Physics, 72(9), 096601. 5.
Gao, W., Kagan, D., Pak, O.S., Clawson, C., Campuzano, S., Chuluun‐Erdene, E., Shipton, E., Fullerton, E.E., Zhang, L., Lauga, E. and Wang, J. (2012). Cargo‐towing fuel‐free magnetic nanoswimmers for targeted drug delivery. small, 8(3), 460-467. 6.
Mhanna, R., Qiu, F., Zhang, L., Ding, Y., Sugihara, K., Zenobi‐Wong, M., & Nelson, B. J. (2014). Artificial bacterial flagella for remote‐controlled targeted single‐cell drug delivery. Small, 10(10), 1953-1957. 7.
Magdanz, V., & Schmidt, O. G. (2014). Spermbots: potential impact for drug delivery and assisted reproductive technologies. Expert Opinion On Drug Delivery 11(8),1125-1129. 8.
J. Happel and H. Brenner, Low Reynolds Number Hydrodynamics. Prentice-Hall, 1965. 9.
Purcell, E. M. (1977). Life at low Reynolds number. American journal of physics, 45(1), 3-11. 10.
Dreyfus, R., Baudry, J., Roper, M. L., Fermigier, M., Stone, H. A., & Bibette, J. (2005). Microscopic artificial swimmers. Nature, 437(7060), 862-865. 11.
Jang, B., Gutman, E., Stucki, N., Seitz, B.F., Wendel-García, P.D., Newton, T., Pokki, J., Ergeneman, O., Pané, S., Or, Y. and Nelson, B.J. (2015). Undulatory locomotion of magnetic multilink nanoswimmers. Nano letters, 15(7), 4829-4833. 12.
Tottori, S., Zhang, L., Qiu, F., Krawczyk, K. K., Franco‐Obregón, A., & Nelson, B. J. (2012). Magnetic helical micromachines: fabrication, controlled swimming, and cargo transport. Advanced materials, 24(6), 811-816. 13.
Walker, D., Kubler, M., Morozov, K. I., Fischer, P., & Leshansky, A. M. (2015). Optimal length of low Reynolds number nanopropellers. Nano letters, 15(7), 4412-4416. 14.
Morozov, K. I., & Leshansky, A. M. (2014). Dynamics and polarization of superparamagnetic chiral nanomotors in a rotating magnetic field. Nanoscale, 6(20), 12142-12150. 15.
Suter, M., Zhang, L., Siringil, E.C., Peters, C., Luehmann, T., Ergeneman, O., Peyer, K.E., Nelson, B.J. and Hierold, C. (2013). Superparamagnetic microrobots: fabrication by two-photon polymerization and biocompatibility. Biomedical microdevices, 15(6), 997-1003. 16.
Roper, M., Dreyfus, R., Baudry, J., Fermigier, M., Bibette, J., & Stone, H. A. (2008, April). Do magnetic micro-swimmers move like eukaryotic cells?. Proceedings of the Royal Society of London A: Mathematical, Physical and Engineering Sciences (Vol. 464, No. 2092, pp. 877-904). The Royal Society. 17.
Gauger, E., & Stark, H. (2006). Numerical study of a microscopic artificial swimmer. Physical Review E, 74(2), 021907. 18.
Nayfeh, A. H., & Mook, D. T. (2008). Nonlinear oscillations. John Wiley & Sons. Chicago 19.
Kapitza, P. L. (1951). Dynamic stability of a pendulum with an oscillating point of suspension. Zh. Eksp. Teor. Fiz, 21(5), 588-597. 20.
Cox, R. G. (1970). The motion of long slender bodies in a viscous fluid Part 1. General theory. Journal of Fluid mechanics, 44(04), 791-810. 21.
Gray, J., & Hancock, G. J. (1955). The propulsion of sea-urchin spermatozoa. Journal of Experimental Biology, 32(4), 802-814. 22.
Gutman, E., & Or, Y. (2014). Simple model of a planar undulating magnetic microswimmer. Physical Review E, 90(1), 013012. 23.
Nayfeh, A. H. (2008). Perturbation methods. John Wiley & Sons. 24.
Sudor, D. J., & Bishop, S. R. (1999). Inverted dynamics of a tilted parametric pendulum. European Journal of Mechanics-A/Solids, 18(3), 517-526. cknowledgment
The work of D.J. and L.Z. was supported by the Early Career Scheme (ECS) grant with Project No. 439113 and the General Research Fund (GRF) with Project No. 14209514 and 14203715 from the Research Grants Council (RGC) of Hong Kong SAR. The work of Y.H. and Y.O. was supported by the Israeli Science Foundation (ISF) under grant no. 567/14.
Author contributions
Y.O. and L.Z. initiated the research. Y.O. and Y.H. performed the numerical simulations. L.Z. and D.J. designed and fabricated the microswimmer. L.Z. and D.J. designed and provided the experimental magnetic setup for manipulation. D.J. conducted the magnetic actuation experiments and performed the analysis of the dynamic motion of microswimmers. Y.O. and Y.H. proposed the two-link model and achieved the theoretical analysis. D.J., Y.H., L.Z. and Y.O. cowrote the manuscript. Y.O. and L.Z. supervised the work and gave critical input. All authors contributed to discussions.
Additional information
Competing financial interests
The authors declare no competing financial interests. upplementary Information for manuscript: "Nonlinear parametric excitation effect induces stability transitions in swimming direction of flexible superparamagnetic microswimmers"
Yuval Harduf, Dongdong Jin, Yizhar Or and Li Zhang
Part 1 - Experiments
This part contains detailed information on the experiments. In particular, it presents the fabrication process and characterization of the flexible superparamagnetic microswimmers, as well as the magnetic actuation experiments.
The chemicals of iron(III) chloride hexahydrate (99%), iron(II) chloride tetrahydrate (99%), sodium citrate dihydrate (99%), ammonia solution (25-28%), citric acid monohydrate (99.5%), pyrrole (99%) and sodium hydroxide (97%) were purchased from Aladdin Chemicals. Anodic aluminum oxide (AAO) template with average pore diameter of ~200 nm and thickness of 60 μm was purchased from Shanghai Shangmu Technology. All the chemicals were used without further purification.
Synthesis and characterizations of Fe O nanoparticles The Fe O nanoparticles were prepared by hydrothermal method. In a typical procedure, 30 mL of ultrapure water (18.2 MΩ*cm) was first purged with dry argon gas in a round-bottom flask while stirring for 30 min to remove the oxygen in water. Subsequently, 0.5 ml of concentrated hydrochloric acid, 30 mmol of FeCl ·6H O, 15 mmol of FeCl ·4H O and 30 mmol of sodium citrate dihydrate were added into the flask. After the chemicals were completely dissolved, the pH value of obtained solution was adjusted as 9 by gradually adding ammonia, affording a black precipitate. The black product was poured into a teflon-lined stainless steel autoclave which was then sealed and heated at 180 ℃ for 2 h. Finally, the anoparticle product was collected by centrifugation at 10000 rpm for 5 min and washed with ultrapure water for three times, followed by redispersing in ultrapure water prior to use. Before fabricating microswimmers, the morphology, size distribution, magnetic property and surface charge distribution of the as-prepared Fe O nanoparticles were investigated (shown in Supplementary Fig. 1(a-d)) using transmission electron microscope (Tecnai Spirit 12, FEI), Zeta Sizer (Nano ZS, Malvern), and vibrating sample magnetometer (VSM 7410, Lake Shore). Typical TEM micrograph of the Fe O nanoparticles shown in Supplementary Fig. 1a displayed that Fe O nanoparticles had an average diameter of 10-20 nm. When we dispersed Fe O nanoparticles in water, however, aggregation of Fe O nanoparticles would occur, leading to the increasing size of Fe O particles. It could be observed from Supplementary Fig. 1b that the actual diameter of most Fe O nanoparticles (~68.4%) in electrophoresis bath was between 40 and 100 nm, and the average size of overall Fe O nanoparticles was calculated as ~80 nm. Therefore, given that the average pore size of AAO template was ~200 nm, the possibility of assembling the as-prepared Fe O particles into AAO template might be realized. Supplementary Fig. 1c showed the magnetic hysteresis loop of the Fe O nanoparticles at room temperature. The Fe O nanoparticles were superparamagnetic with a saturation value of 65.4 emu/g. Moreover, the average zeta potential of Fe O nanoparticles was calculated as ~-50 mV from Supplementary Fig. 1d, demonstrating that Fe O nanoparticles could be electrophoretic deposited in working electrode by applying a positive voltage between working electrode and counter electrode. Synthesis of microswimmers
Microswimmer with a superparamagnetic head was synthesized by electrophoretic codeposition of Fe O nanoparticles and polypyrrole followed by electrochemical deposition of polypyrrole. A thin gold film (~100 nm) was first coated on one side of AAO templete using magnetron sputtering and the AAO template was then assembled in a homemade teflon plating cell (shown in Supplementary Fig. 2a) to serve as a working electrode. Subsequently, he plating cell was put in a modified vacuum drier as shown in Supplementary Fig. 2b. For the electrophoretic deposition of Fe O nanoparticles based superparamagnetic head, 10 mL of aqueous solution containing ~1 mg of Fe O nanoparticles and 1 mmol of pyrrole acting as electrophoresis bath was added into the plating cell under vacuum to facilitate the nanoparticles to enter into the pores of AAO template. A thin stainless steel sheet was employed as a counter electrode, and the distance between working electrode and counter electrode was 10 mm. Then the electrophoretic deposition was performed at a constant electric potential of +50 V for 3 h. After the reaction was completed, the plating solution was poured and the assembled AAO template was washed by being immersed in ultrapure water for 30 min. Next, for the electrochemical deposition of polypyrrole nanowire based flexible tail, a Pt wire and an Ag/AgCl electrode were used as counter and reference electrodes, respectively. Then polypyrrole nanowires were electropolymerized at +0.80 V for 3 C from a plating solution containing 0.1 M citric acid and 0.1 M pyrrole. Finally, the AAO template was taken out from plating cell and dissolved in 3M NaOH for 30 min to obtain the microswimmers, which was then collected by centrifugation at 10000 rpm for 3 min and washed with ultrapure water for three times. All microswimmers were stored in ultrapure water at room temperature when not in use. Typical scaning electron microscopy (SM-7800F, JEOL) images of as-prepared microswimmers were shown in Supplementary Fig. 1e, f. Magnetic actuation experiment
The magnetic control experiments were conducted using a 3-axis Helmholtz electromagnetic coil setup as shown in Supplementary Fig. 3. The control signals were generated by a Sensoray 826 card and then the current was amplified to generate magnetic fields in the coils. There were three pairs of coils in the setup and in our experiments, the current through the middle coils was constant to generate a constant x vector, the current through the big coils was alternating to generate an oscillating y vector, while the current through the small coils was 0. In this way, a planar os cillating magnetic field could be generated. efore the swimming experiments, microswimmer was put in an open tank filled with ultrapure water by a pipette which was then put in the centre of the coil pairs. During the actuation experiments, a constant magnetic field of 5 mT in x direction (
𝐵⃗ = 𝐵 𝑥 ⃗⃗⃗⃗ ) was first set for 5s acting as the initial condition. Then an oscillating y vector with a certain value of β and 0.1 Hz of f was generated together with the constant x vector. We gradually increased the frequency and tried to look for the 90° transition in the oscillating orientation and swimming direction of microswimmer and calculate the swimming velocity. Subsequently, we turned off the magnetic field, set a constant x vector for 5s renewedly, and then changed the value of β , and conducted the magnetic actuatuion experiments in the same way. Each status with a certain value of β and f was kept for at least 5s. art 2 - Theory This part contains detailed information on the theoretical analysis. In particular, it presents the multi-link model of the microswimmer and its numerical simulations, as well as the simplified two-link model and the asymptotic analysis of its dynamics.
Our model of the swimmer as presented in Supplementary Fig. 6 consists of a chain of n slender links ("tail"), attached to a non-magnetic spherical "head", which represents the load. The links are connected by rotary joints with linear torsion springs, and carry superparamagnetic beads. The swimmer is submerged in a viscous fluid and is actuated by an external magnetic field. While Dreyfus's model [1] considered the tail as a continuous filament of length L , we treat it here as a chain of discrete, rigid links, with equal lengths Ll n and a circular cross-section of radius a . In both cases the red blood cell is modeled as a rigid sphere of radius h r . We define a world-fixed reference frame whose axes are x̂, ŷ, ẑ . Only planar motion of the swimmer in x̂ − ŷ plane is considered, and all rotations are about ẑ axis. We attach a reference frame to each link, whose axes are ˆ ˆ ˆ, , i i t n z , where the direction ˆ i t is aligned with the longitudinal axis of the link and ˆ i n is perpendicular to it. The location of each link is defined as its geometrical center and its orientation is defined as the angle between its axis ˆ i t and the ˆ x axis. We denote T x y as the location of the head link in the world-fixed frame, is the angle between the head's axis ˆ t and the ˆ x axis, and i is the relative angle between the axis of the th i link ˆ i t and the axis of the th i link ˆ i 1 t . .1.1 The magnetic torques acting on the links The swimmer is actuated by a uniform, planar magnetic field denoted by t B , which is assumed to vary within the x̂ − ŷ plane and has the form: si 1n x t B B (S1) In Dreyfus et al [1] the paramagnetic beads have a preferred magnetization direction. In the assembly process, it is assumed that all of the beads' directions align with the backbone of the continuous filament. It is also assumed that each bead experiences magnetic interactions only with the two beads adjacent to it. The magnetic torque per unit length in Ref [1] is of the form: (S2)
Where s is a coordinate along the filament's length, χ = , χ ⊥ are the susceptibility components of a single bead in the preferred magnetization direction and in the direction perpendicular to it, respectively, is the vacuum permeability, b r is the beads' radius and ˆ s t is the local unit vector tangent to the filament. In our model, where the continuous filament is regarded as a discrete chain of links, the torque acting upon each link is obtained by integrating (S2) over the length of a link, and has the form: (S3) Where: / 4(1 / 6)(1 / 12)1 (S4) We denote the torques acting at the joints as i . Since the links are connected by torsion springs with a linear spring constant k , the torques acting at the joints are given by: i i k (S5) e would like to determine the springs' constant in such a way that the deflection of our "discrete filament" will be comparable with the deflection of Dreyfus' filament in Ref [1]. We use Euler-Bernoulli beam theory [2] to model the deflection of Dreyfus' filament, and we assume that its properties are constant along its length. Looking at an infinitesimal cross-section of the beam we get: d s M sds (S6) Where is the filament's bending rigidity from Dreyfus's model. We assume that the bending torque is constant along the beam and integrate along the beam's length: l l d ds MlMdsds (S7) And so we obtain our spring constant to be:
Mk Lnl (S8)
Due to the extremely small Reynolds number of the problem, the fluid is governed by Stokes flow, which in turn dictates that the drag forces are linear in the velocity [3] . We model the drag forces using resistive force theory (RFT) [4], [5] which gives the following relation: (S9)
Where , i M i f is the force vector and torque acting upon the th i link, i u is the link's velocity vector, i is the link's angular velocity. , , t n m c c c are the tangential, normal and rotational drag coefficients accordingly. According to resistive force theory [4], [5] , for a slender link these drag coefficients are: n t tm lc c llc ca (S10) And for the spherical head these drag coefficients are [3] : n t m c r rc c (S11) here is the fluid's viscosity. In Dreyfus' model [1] the drag is modeled in the same manner and the drag coefficients are the same, only defined per unit length. It is assumed that the hydrodynamic interaction between links is negligible and so the hydrodynamic forces on a single link are only a result of that link's velocity. We define a body coordinates vector T x y b q and a shape coordinates vector Tn s q . Both vectors compose the full vector of generalized coordinates , . b s q q q We define the orientation of the th i link as the angle between the ˆ i t axis and the ˆ x axis as i and write them in terms of the body and shape coordinates: ii jj (S12) We denote the location of the swimmer as b r , the location of each link as i r and the location of each joint as j b . These vectors are given by:
11 1 221 cos cos2sin sin2cos cos , 3sin sin2 hh i ii i xy lr lrl for i n
1i b2 1i 1 rr rrr r (S13) cos2 sin jj l jj 1 b r (S14) Next, we write the linear and angular velocities of each link in terms of the body and shape velocities, expressed in the body-fixed reference frame: i ii i b i s uV T q E q (S15) here i u is the vector of velocities and i is the angular velocity of the th i link about the ˆ z axis. The matrices i T and i E are given as follows [6] : I I J r rT J r bE E E E E (S16)
Where for i jI for i j ij and J We can now write relation (S9) in matrix form, using "resistance matrices": , i M ihyd i i i fF V R (S17) Where i R are given explicitly as follows. For the spherical link (head): hh rr R (S18) And for the slender links: i i ii i i i l il la i R (S19) Since the swimmer is of micro-scale, the flow's Reynolds number is extremely small. Neglecting all inertial effects, the swimmer moves quasi-statically in equilibrium. Formulating a force and torque balance on the swimmer's body gives: , i Ti hyd i b
T F F (S20)
Substituting (S17) and (S15) into (S20) we get: bb bu i b iTs b q q T F R R (S21)
Formulating torque balance at each joint leads to (see [6]): , i Ti hyd i b s
E F F τ (S22)
Where s n satisfies s - k q s from (S5). Substituting (S17), (S15) and (S5) into (S22) gives : Tbu uu k i Ts i b sb q q E F q
R R (S23)
Where: , , . bb bu uu
T Ti i Ti i ii i ii
T T E E ET
R R R R R R (S24)
It can be shown (see [6]) that the matrices , i i T E in (S20) and (S22) are the same as the ones in equation (S15). Collecting the equations into matrix form, one obtains: ( ) ( ), t bA q q q (S25) where ,( ) , ( ) bb buTbu uu kt ii i b sTTi b A T Fq b q E F q
R RR R (S26)
Equation (S25) is a system of n first-order, time-dependent, nonlinear ordinary differential equations, which are linearly coupled. The solution of this system under given initial conditions can be obtained computationally via numerical integration. All of our numerical analyses were conducted using MATLAB's integration function ODE15s. In order to obtain the mean speed of the swimmer in a period we used the following scheme: 1.
Set physical parameters 2.
Integrate over p n periods of the magnetic field (typically p n is in the range of
100 ~ 1000 , in order for the system to reach steady state) 3.
The mean swimming speed is calculated using the following formula: p pp p x n nn ny V xV xy y (S27) n all of the simulations we used a swimmer whose tail consists of n links, and unless stated otherwise, the numerical values for the physical parameters used in the simulations are taken from the work of Dreyfus et al [1] . For the simulation results in Fig. 1b in the manuscript, these values are summarized in Supplementary Table 1. For the simulation results in Fig. 1c in the manuscript, we used the "green squares" data from Supplementary Table 2, without a spherical head, under varying values of the magnetic field. For the simulation results in Fig. 1d and 1e, we used the "green squares" data from Supplementary Table 3, including a spherical head, under varying values of the magnetic field. We now introduce the two-link model, which is the simplest possible model that captures the microswimmer's dynamics and stability transitions. It consists of two slender links connected by a single rotary joint with a torsion spring, as shown in Supplementary Fig.7. The microswimmer's "head" link is made of superparamagnetic material having uniaxial anisotropy with an easy axis along the longitudinal axis of the link 𝒕̂ . The susceptibility tensor of the uniaxial anisotropic link is: T I t Iχ t (S28)
The magnetic field B polarizes the head link and induces a magnetic moment: v χ B M (S29) Where v is the link's volume. The magnetic field generates a torque acting on the polarized link: L = M× B = ∆χ v ( t ̂ × B )( t ̂ ∙ B ) =𝐿 𝑚 𝐳̂= ∆χvB sin (2 γ ) 𝐳̂ (S30) Where 𝛾 is the angle between the magnetic field's vector B and the link's direction ˆ t (see Supplementary Fig.7), and B is the magnitude of the magnetic field B . The magnetic field ( ) t B undergoes planar oscillations as given in (S1). ormulation of the dynamic equations of the two-link microswimmer follows the same procedure as described above for the multi-link swimmer. Some differences are that there is no spherical link, and that the torque on the magnetic link is given by (S30) rather than (S3). Using this, one can obtain a system of differential equations of the form (S25), where the generalized coordinates are ( , , , ). x y q An additional simplification is obtained if one replaces body velocity , x y by the components , t n v v expressed in the body-fixed reference frame ˆ ˆ, t n . The kinematic relation between these velocity components is given by: cos sinsin cos tn vx vy (S31) The use of body-fixed velocity components , t n v v eliminates the dependence of the left-hand side of (S25) on the absolute orientation angle , resulting in considerable simplification, as: )( ) , ,( tn vv t A b (S32)
Where: t l ll ll l l lc l ll l A
2) 4) 300(cos( ) sin( )sin( ))( cos( )sin( ) sin( )) x lvB t tk b After inversion of the matrix A in (S32), the dynamic equations of the two-link microswimmer are obtained as: sin( ) cos( ) sin(2 ) sin (2 ) 1 2 cos(2 )sin(2 ) sin( )(cos( ) 3)2(cos(2 ) 17) 2(cos(2 ) 17)sin (cos( )(5 cos(2 )) sin(2 ) cos(4 ) 2 4 cos(2 )sin(2 ) 216(cos(2 ) 17) t m kn m l lv t tltft ftft ftv m km k ft ftft ltt ttft (cos( ) 3)(cos(2 ) 1 17) k t (S33) Where 𝑡 𝑚 , 𝑡 𝑘 are characteristic time scales given by: , 6 12 t tm kx l viscous c l viscoust tv magnetic k elast cc B i (S34) Note that (S33) is independent of the 𝑥, 𝑦 coordinates, since the fluid domain is unbounded. Therefore, the dynamic behavior of the swimmer is mostly encapsulated in the coupled nonlinear second-order system of the two angles , . We study it below using asymptotic analysis under three different physical limits. We now consider the dynamics in (S33) in the limit of small-amplitude oscillations of the magnetic field B ( t ) in (S1). By assuming <<1, we can expand the solution into a power series in as: ... t t t t q q q q (S35) We expand the expressions for , from (S33) into a Taylor series about
0, 0 and substitute (S35) into the expanded expressions. Next, we equate coefficients of same powers of in order to obtain a system of differential equations for each order of approximation. The first-order approximation yields the following system:
11 11 m k mm k m t t tt tt t (S36)
This is a linear system whose characteristic polynomial is:
1) 818( 5 m k m k tt t t
Since , m k t t are positive values, the system is stable. This implies that transient terms in the solution of (S36), which take the form e - t, decay to zero in time. The steady-state solution of the system in (S36) contains harmonic terms which depend on the frequency , and can be obtained as: k m mk k mk m k m m k m k m mk k mk k mk m k m m k m k m t t t tt t tt t t t t t t t t tt t t tt t tt t t t t t t t tt tt cos4 1 m tt (S37) Next, it is possible to obtain the first-order solution for the body-fixed velocities , t n v v in (S33) and expand the relation (S31) in in order to find the leading-order dynamics of forward motion x ( t ), which is of order 2, and given by: sin( t) sin( t)16 16 8 16 16( ) m m m mk l l l l lt t tx t tt (S38) Substituting the expressions for , from (S37) into (S38), one obtains: ( ) sin 2 cos x x x x Ot tt V DC (S39) Where
22 2 2 2 24 2 6 2 4 2 2 2 2 222 2 2 2 22 3 4 2 2 22
64 25 64 64 1832 325 1408 1152 2 29 64 32 14 64 25 64 64 14 (65 112 ) 45 120 64 3 4 kx k m k m mk m k k k m m k k m mx k m k m mk m k m k k k m m k mx k t lt t t t tt t t t t t t t t t tC lt t t t tt t t t t t t t t t tD tV
22 2 2 2
64 25 64 64 1 m k m m lt t t t
Note that the expression for the mean speed x V implies that it vanishes for extreme frequencies and maximized for some intermediate optimal frequency opt which can be obtained rather simply: opt m k lt t (S40) Importantly, the small- approximation gives zero net motion in y direction. Moreover, the leading-order system in (S36) is always stable, and expansion about the perpendicular orientation = /2, =0, yields a linear system which is always unstable . Thus, the small-amplitude asymptotic approximation cannot capture the stability transitions and changes in the swimming direction from x to y . We now analyze the two-link microswimmer's dynamics in the limit of fast oscillation frequency of the magnetic field, by using the method of multiple scales [7],[8] . We denote ϵ = (ω𝑡 𝑚 ) -1 , ϵ= (ω𝑡 𝑘 ) -1 and assume that 𝜖 <<1 and = O (1). We define the nondimensiona l time t t , and introduce fast and slow time scales 𝑇 , 𝑇 as: T = t ̃ = ωt , T = ϵt ̃ = t / t m We expect to obtain solutions of fast oscillations superposed on slow evolution of "average" values, and so we assume solutions of the form: , ( ), ,
T TT T , and expand them into power series in 𝜖 as: (S41) From the chain rule, time-differentiation is expanded as: (S42) and hence we define the time-derivative operators: n n
D T for n =0,1. Substituting (S41) and (S42) into the equations for in (S33), one obtains: (S43) We expand (S43) into a Taylor series in , about 0 and equate coefficients of powers of 𝜖 . Equating coefficients of order 𝜖 gives: D TD T (S44)
That is, 𝜃 and 𝜙 are functions of the slow time 𝑇 only. Equating coefficients of order 𝜖 in the expansion of (S43) gives: T TD D T TD D (S45)
Substituting the solutions for , from (S44) into (S45) gives: D D T TT TD D (S46)
Integrating (S46) with respect to the fast time T gives: sin(2 )(cos(2 ) 19) cos(2 )(cos(2 ) 19)sin(2 ) cos( )16(cos(2 ) 17) 2(cos(2 ) 17)2 2 sin(2 )(cos(2 ) 19) 8(cos( , ) ( ) 3)(cos( )(4 ) 3(4 )) ( )16(cos(2 ) 17) T T T TD D T T sin( ) cos( )(cos( ) 3) cos(2 )(cos( ) 3)sin(2 ) cos( )8(cos( ) 3) 2(cos( ) 3)4(3 cos( )( 2 ) 6 ) 2 2 sin( ) cos( )(cos( ) 3) (8(cos( ) 3)( , ) )
T T T TD D T T (S47) ince we assume an oscillating solution in the fast time T , we require that the diverging secular terms (multiplied by T ) in (S47) are eliminated. We obtain the following system for the zero-order approximation of the slow dynamics of , : ddTddT (S48) This is a second-order coupled nonlinear system of time-independent differential equations in the slow time T . This system has two equilibrium points: { 0, 0} e e and { , 0}2 e e (more precisely, e k for 𝑘 = 0,1,2 … but the all other value are ignored due to symmetry). Linearization about e yields:
21 0 021 0 0 ee DD (S49) The characteristic polynomial for this system is:
12 cos 2 2 cos 216516 e e
In order for a 2 nd order system to be locally asymptotically stable all the coefficients of its characteristic polynomial are required to be positive, hence for e the system is stable when , while for e the system is stable when . Therefore, steady-state solutions for , depend on the value of as:
20, 2 20 forfor (S50) ext, we look for approximate solutions of higher orders for , in steady state for oscillations about either e or e . Substituting the steady state value of , (S50) into (S47) yields: ee (S51) From (S33), the next-order approximation is: e e ee e e
D T T TD T T
Integrating in T gives e e ee ee e e T TT TT T
02 21 0 1 1 1 1 0 2 1 n(3 )1 1cos(2 ) sin(2 )+ + -2 cos(2 )-D4 2 e e
TT T T (S52)
Requiring elimination of (diverging) secular terms yields: ee DD (S53) This system (S53) is similar to the system in (S49) and has the same stability conditions. The steady-state solution of the system (S49) for a stable point 𝜃 𝑒 converges to zero. Substituting the steady state solution into (S52) yields: e e ee e e TT (S54) urther expansions of the solutions to higher orders is done iteratively in a similar manner, up to any desired order in 𝜖 . The approximate steady-state solution is then obtained by substituting the steady-state solution of each order into (S41): (S55) Next, the expressions in (S55) are substituted into (S33) in order to obtain leading-order expressions for the body velocity components , t n v v . These body-frame velocities are transformed into world-frame velocities , x y by using the relation (S31), expanding it as a power series in ( , ) about ( , 0) e , and substituting (S55). The leading-order expressions for , x y contain oscillating terms and constant terms V x, V y , which are precisely the mean speed components, given by: (S56) It can be seen from (S56) that in the case where oscillations about e are stable the swimmer moves in ˆ x direction, x y V V . On the other hand, in the case where oscillations about e are stable the swimmer moves in ˆ y direction, x y V V . Finally, the expressions in (S56) imply the existence of an optimal value of that maximizes the speed V y in the regime of . .2.3 Asymptotic analysis under a stiff spring We now consider the limit where the torsion spring is very stiff and oscillations frequency are very fast. This implies that ϵ = (ω𝑡 𝑚 ) -1 , where ϵ <<1, while ( ωt k ) -1 = O (1). We use the same nondimensional time t t , and expand the joint angle into a power series in ϵ as: (S57) Substituting (S57) into the system (S33), the equation for is obtained as: (S58) For the equation in from (S33), we equate expressions in different of powers of ϵ : (S59) It can easily be seen that the solution for t decays asymptotically to zero, since the function f in (S59) is bounded between two positive constants. Thus, in steady-state we can use the first-order approximation 𝜙 = ϵ 𝜙 . Substituting into (S33) and re-normalizin g time and frequency by m t , one obtains the system: t tt t (S60) The system (S60) is nonlinear in t while t appears linearly. Using the time-differentiation operator D , we substitute D into (S60) and eliminate in order to obtain a single equation for as: D t tD t D t
This equation can be rearranged as a second-order Kapitza-like differential equation: θ ̈ + ( α +10 cos ( θ ) ( β cos ( ωt ) - β +2 ) +40 β sin ( θ ) sin ( ωt )) θ ̇ + ( α ( β ) + αβ cos ( ωt ) -10 ωβ sin ( ωt )) sin ( θ ) =4(5 ωβ cos ( ωt ) + αβ sin ( ωt ) ) cos ( θ ) (S61) The nonlinear system in S(61) typically has a periodic solution of nonlinear oscillations , where 2 / . p p t t T T , We now study the stability of this periodic solution by considering a small perturbation about it: ( ) ( ) ( ) p t t t (S62) Substituting (S62) into S(61), expanding trigonometric expressions and rearranging gives:
222 2 22 22 p pp p pp pp p tt tt t
22 2 22 2
5t ) +2 sin(2 ) sin(t )41 cos(2 ) 2 sin(2 ) sin(2 )32 161 5sin( ) cos( ) cos(2 )8 85 sin(2 ) sin (t ) 1 +2 cos(2 ) sin(t )5 04 p pppp p t tt t (S63) he terms in the square brackets in (S63) vanish since p t is a solution of S(61). Assuming small perturbation t , the higher-order terms appearing in braces in (S63) also vanish, and we are left with the variational equation:
222 2 22 2 p pp p pp p tt tt t (S64)
This is a Hill-type equation [8] , a second-order linear homogenous differential equation with time-periodic coefficients [8] . According to Floquet theory [8] , this equation has two fundamental solutions , t t , which satisfy i i i t T t for
1, 2 i , where i are called "Floquet multipliers". The conditions for local orbital stability of the periodic solution p t is thus reduced to the requirement that i for
1, 2 i , which implies that any small perturbation t decays asymptotically to zero. The Floquet multipliers can be found numerically as follows. The linearity of (S64) implies that its solution satisfies , where . t tT z Pz z = (S65) Columns of the matrix P can be obtained from solutions of T z under initial conditions of T z and T z . The Floquet multipliers are simply the eigenvalues of 𝐏 . Thus, the stability of a periodic solution can be obtained by numerical integration of (S64) over a single time period T , once the periodic solution p t of S(61) is found numerically. References [1]
Dreyfus, R., Baudry, J., Roper, M. L., Fermigier, M., Stone, H. A., & Bibette, J. (2005). Microscopic artificial swimmers. Nature, 437(7060), 862-865 . [2] Gere, J. M., & Timoshenko, S. P. (1997). Mechanics of materials, 1997. PWS-KENT Publishing Company, ISBN 0, 534(92174), 4 . Chicago. 3]
J. Happel and H. Brenner, Low Reynolds Number Hydrodynamics. Prentice-Hall, 1965 . [4] Cox, R. G. (1970). The motion of long slender bodies in a viscous fluid Part 1. General theory. Journal of Fluid mechanics, 44(04), 791-810. [5]
Gray, J., & Hancock, G. J. (1955). The propulsion of sea-urchin spermatozoa. Journal of Experimental Biology, 32(4), 802-814. [6]
Wiezel, O., & Or, Y. (2016). Optimization and small-amplitude analysis of Purcell's three-link microswimmer model.
Royal Society, Proceedings A . 472: 20160425, 2016. [7]
Nayfeh, A. H. (2008). Perturbation methods. John Wiley & Sons . [8] Nayfeh, A. H., & Mook, D. T. (1995).
Nonlinear Oscillations.
Wilet-VCH. upplementary Figure Legends:
Supplementary Figure 1|Characterizations of Fe O nanoparticles and microswimmers. a , TEM micrograph of Fe O nanoparticles. b , Size distribution of Fe O nanoparticles. The standard deviation values were obtained by repeating the measurement for three times. c , Magnetic hysteresis loop of Fe O nanoparticles. d , Surface charge distribution of Fe O nanoparticles. e-f , SEM images of the soft microswimmers. Supplementary Figure 2|Illustration of the homemade teflon plating cell (a) and modified vacuum drier (b).
Supplementary Figure 3|Setup of the dynamic magnetic field generation and actuation system, which consists of three pairs of Helmholtz electromagnetic coils assembled in an orthogonal manner
Supplementary Figure 4|The 90°-transition in oscillating direction of microswimmer when β was 3 and f was increased from 0.3 to 0.5 Hz. Movie clips are available in the Supplementary Movie 4.
Supplementary Figure 5|The dynamic motion of microswimmer in y direction when β was 2.5 and f was 3 Hz. Full movies were enclosed in Supplementary Movie 5.
Supplementary Figure 6|Schematic model of the multi-link swimmer.
Supplementary Figure 7|Model of the two-link microswimmer. upplementary Table 4|List of parameters used for numerical analysis, taken from Dreyfus et al [1]
PARAMETER MEANING
GREEN SQUARES ■ BLUE DIAMONDS RED CIRCLES Flexural rigidity −22 [ 𝐽𝑚] −22 [ 𝐽𝑚] −22 [ 𝐽𝑚] L Filament length x B Constant field component y B Oscillating field component amplitude s a Blood cell radius 𝝌 = Susceptibility in the preferred direction Susceptibility in the perpendicular direction n c Normal drag coefficient a Bead radius0.5[𝜇𝑚]