Modeling Propulsion of Soft Magnetic Nanowires

1. Introduction

Development of artificial nanomachines that can controllably propel through complex fluidic environments is one of the most exciting challenges of nanotechnology. The emergent approaches range from catalytically driven nanowires, Janus particles to thermally, light- and ultrasound-driven colloids (Wang, 2013). One of the promising methods is magnetic actuation. While the traditional techniques are based on strong gradient magnetic fields to generate a force for remote towing of magnetic nanoparticle, the alternative approach relies on a weak (milli Tesla) uniform rotating magnetic field that serves to apply a torque twirling the nanomotor. Given that the particle shape admits non-trivial rotation-translation coupling, such torque-driven twirling will result in net propulsion. Notice that while the field gradient can tow an isotropic (i.e., spherical) particle, it is not very efficient due to the need of large field strength (of the order of Tesla) required for creating appreciable variance of the field at the size of the nanoparticle. Although propulsion based on external magnetic torque requires more complex particle shape [e.g., helical (Ghosh and Fischer, 2009; Zhang et al., 2009)], it offers a remote, fuel-free and engineless propulsion in a variety of fluidic environments with typical velocities considerably exceeding the speed of gradient towing.

This technology has been extensively studied over the last decade by a number of groups. Various methods, such as “top-down” approach (Zhang et al., 2009), delamination of magnetic stripes (Smith et al., 2011), glancing angle deposition Ghosh and Fischer (2009), direct laser writing (Tottori et al., 2012), biotemplated synthesis using biological spiral organelles (Gao et al., 2014), two-photon polymerization of a curable magnetic polymer composite (Peters et al., 2013), spiraling microfluidic flow lithography (Yu et al., 2017), and other techniques have been developed for fabrication of μm-size and sub-μm-size (Schamel et al., 2014) helical motors. These bio-inspired helical motors propel unidirectionally (along the field rotation axis) when driven by a rotating magnetic field similar to a twirling bacterial flagellum. However, fabrication of three-dimensional (3D) chiral/helical nanostructures requires sophisticated procedures and the search for simpler alternatives is under way.

One alternative relies on the fact that in torque-driven propulsion the propeller's shape may not necessarily be helical or even chiral. It was recently demonstrated that geometrically achiral objects made of three interconnected magnetized microbeads can be steered quite efficiently by an in-plane rotating magnetic field (Cheang et al., 2014). These findings suggested that the two-dimensional (2D) ferromagnetic propellers could be of practical interest, as they can be mass-fabricated via standard photolithography methods (Tottori and Nelson, 2018). It was theoretically predicted (Morozov et al., 2017) and demonstrated experimentally [using up-scaled cm-size propeller (Sachs et al., 2018)] that unidirectional propulsion of 2D magnetic structures is feasible, however it requires non-trivial off-plane magnetization. Since 2D structures are prone to magnetize in-plane, uniform off-plane magnetization at nano/microscale cannot be easily achieved.

Another possibility to avoid sophisticated microfabrication relies on spontaneous aggregation of magnetic nanoparticles into random-shaped 3D microclusters (Vach et al., 2013, 2015). These random aggregates can also be steered through fluid by an external torque, however they appear to be significantly less efficient “swimmers” on average in comparison to nanomotors with preprogrammed (optimal) geometry and magnetization (Mirzae et al., 2018).

There is an additional possibility that uses even simpler one-dimensional (1D) soft magnetic nanowires. The first demonstration of soft artificial “swimmer” was provided by Dreyfus et al. (2005) whereas a linear chain of magnetic microbeads linked by DNA and attached to a red blood cell was actuated by a plane oscillatory magnetic field. This “swimmer” was undergoing in-plane undulations and propelled similar to a flagellum of an eukaryotic cell. A minimal design of the planar undulating magnetic microswimmer made of just two rigid links connected by a torsional spring was suggested by Gutman and Or (2015) and the corresponding nanowire-based analog was demonstrated by Jang et al. (2015). Highly efficient two-arm magnetic nanoswimmer exhibiting complex 3D (“freestyle”) undulations driven by an in-plane oscillatory magnetic field was reported in Li et al. (2017).

Flexible 1D nanowire-based propellers steered by a rotating field were reported by Gao et al. (2010). These nanomachines were fabricated by electrodeposition from a composite nanowire that had a rigid magnetic (Ni) head, flexible (porous Ag) middle segment and passive rigid (Au) tail. The flexible segment deforms and supposedly acquires helicity due to an interplay of viscous and elastic forces when actuated by rotating magnetic field and the nanowire propels similar to the 3D rigid helical motors. The nanowire propeller in (Gao et al., 2010) had a total length 6.5 μm and diameter 100 nm and exhibited propulsion with speed of ~5 μm/s when actuated by the magnetic field rotating in-plane with frequency 10–15 Hz. Pak et al. (2011) proposed a similar two-segment (Ni-Ag) design whereas nanowire propeller was driven by a conically rotating magnetic field (permanent magnetic field applied along the axis of the rotating field). This Ni-Ag propeller showed even faster propulsion in comparison to Gao et al. (2010). Interestingly, the three-segment nanowire of Gao et al. (2010) could propel either head-forward (fast) or tail-forward (slow), whereas the two-segment nanowire in Pak et al. (2011) could only propel head-forward.

Theoretical modeling of soft nanowire propellers is quite limited. Gauger and Stark (2006) put forward a discrete bead-spring model for simulating externally driven undulatory propulsion of a flexible filament powered by an oscillating magnetic field reported by Dreyfus et al. (2005). Only planar undulations (via bending) were considered and the hydrodynamics was modeled using approximate Rotne-Prager method. Pak et al. (2011) put forward an approximate elastohydrodynamic model based on the dynamic deformation of initially straight elastic filament that is steadily rotated at its one (clamped) end, while the other end is free. The filament is rotated in a way, that its long axis forms a constant conical angle with the rotation axis mimicking actuation by a conically rotating magnetic field. Although this model showed a reasonable agreement with the experiments, it only imitates the real problem in which both ends of the nanowire are free. Various analytical theories concerning two-, three-, or multi-link propellers take advantage of simplified (local) hydrodynamics that neglects hydrodynamic interaction between the rigid links (see e.g., Gutman and Or, 2015; Jang et al., 2015; Alouges et al., 2019).

In the present paper we develop a numerical approach for modeling of flexible, semi-flexible or multi-link externally driven nanopropellers. The numerical algorithm is based on a discrete bead-spring model of the filament. Multipole Expansion (ME) algorithm allows for rigorous account of non-local intra-filament hydrodynamic interaction mediated by a viscous fluid, while springs connecting the beads model the elastic deformation (due to bending and torsion). The ability to accurately model viscous hydrodynamic forces is important for accurate simulations of non-slender objects and large-amplitude deformations (see e.g., Berman et al., 2013 for the importance of non-local hydrodynamics in locomotion powered by large-amplitude undulations). This is particularly relevant for externally actuated nanowire propellers, whereas the deformation is not known in advance and determined by an interplay of elastic, magnetic and hydrodynamic forces. Below we demonstrate the applicability of the proposed approach toward simulations of soft nanowire-based propellers (Gao et al., 2010; Pak et al., 2011).

2. Results

2.1. Bead-Spring Model of the Flexible Filament

2.1.1. Elastic and Magnetic Energies

The bead-spring model of the flexible filament comprises of N identical spherical beads of radius a (beads of different radii a_i can also be used) connected by N − 1 springs. Each bead is characterized by its position r_i and a triad of orthogonal (binormal, normal, and tangent) unit vectors {b_i, n_i, t_i} (Maggs, 2000; Wada and Netz, 2012), see the schematic illustration in Figure 1. The ith bead is linked to its neighbors by springs connected to its surface at points r_i − a_it_i and r_i + a_it_i.

FIGURE 1

Figure 1. Illustration of the bead-spring model: magnetized beads connected by springs actuated by an in-plane rotating magnetic field. Triads of local orthogonal unit vectors {b_i, n_i, t_i} affixed to each bead together with magnetic moments m_i are showing.

There are three contributions into the elastic energy: stretching, bending and torsion, $E^{el}=E_{st}+E_b+E_t$. The stretching energy reads

$\begin{array}{ll}{E}_{st}=\frac{K}{2}{\displaystyle \sum _{i=1}^{N-1}}{\left(|r_i-a_i{t}_i-r_{i+1}-a_{i+1}{t}_{i+1}|-l_0\right)}^2,& (1)\end{array}$

where K is the spring constant and l₀ is the length of the spring at rest. The constant K is assumed to be sufficiently large to preserve the curvilinear length of inextensible or weakly extensible filament. We use non-zero values of the rest lengths l₀ to prevent the beads from overlapping and avoid the necessity to introduce steric repulsions between neighboring beads. In calculations we typically used the values l₀ ~ 0.25a. Notice that K should not necessarily be constant along the filament. When the filament is made of several segments (e.g., of different rigidity), the value of the spring constant can be prescribed separately to each segment.

The bending energy describes the preferable orientation of the tangential vector t_i with respect to the line of centers r_i,i+1 connecting the neighboring beads. The characteristic bending energy of a bead-spring contact is $\frac{A}{a}(1-t_i·{\hat{r}}_{i,i+1})$, where A is the bending modulus and ${\hat{r}}_{i,i+1}$ is the unit vector along the line of centers of i and i + 1 beads, ${\hat{r}}_{i,i+1}=\frac{r_i-r_{i+1}}{|r_i-r_{i+1}|}$. Summing up all bead-spring contacts and symmetrizing the bending energy relative to both filament ends, one obtains

$\begin{array}{ll}{E}_b=\frac{A}{2a}{\displaystyle \sum _{i=1}^{N-1}}\left[2-(t_i+t_{i+1})·{\hat{r}}_{i,i+1}\right],& (2)\end{array}$

where A is the bending modulus.

Notice that the bending modulus A was chosen in such a way that the relation (2) would reduce to the form $E_b=(A/2a)\sum _i(1-t_i·t_{i+1})$ commonly used for chains of beads in contact (i.e., for l₀ = 0) (Maggs, 2000; Wada and Netz, 2012), for small deviations of t_i and t_i+1 from the vector r_{i, i+1}. Let us choose the local spherical coordinate system with the polar axis aligned with the unit vector ${\hat{r}}_{i,i+1}$, i.e., ${\hat{r}}_{i,i+1}=(0,0,1)$. In this frame the vectors t_i and t_i+1 take the form t_i = (sinθ_icosϕ_i, sinθ_isinϕ_i, cosθ_i) and t_i+1 = (sinθ_i+1cosϕ_i+1, sinθ_i+1sinϕ_i+1, cosθ_i+1). Neglecting the cross terms in the scalar product t_i·t_i+1, we have $t_i·t_{i+1}\cong cos{\theta }_{i}cos{\theta }_{i+1}\approx 1-\frac{1}{2}{\theta }_i^2-\frac{1}{2}{\theta }_{i+1}^2$. The expression in the square brackets in Equation (2) reads $2-(t_i+t_{i+1})·{\hat{r}}_{i,i+1}=1-cos{\theta }_i-cos{\theta }_{i+1}$$\approx \frac{1}{2}{\theta }_i^2+\frac{1}{2}{\theta }_{i+1}^2$. Thus, in the limit of small deformations we have $2-(t_i+t_{i+1})·{\hat{r}}_{i,i+1}\approx {1-t}_i·t_{i+1}$.

The torsional energy can be written in the form (Maggs, 2000; Wada and Netz, 2012)

$\begin{array}{ll}{E}_t=\frac{C}{4a}{\displaystyle \sum _{i=1}^{N-1}}(1-n_i·n_{i+1}-b_i·b_{i+1}+t_i·t_{i+1}),& (3)\end{array}$

where C is the twisting modulus. The coefficients A and C have the units of energy times length, K has units of energy per area, while their respective magnitudes typically satisfy Ka³ ≫ A ≈ C.

The magnetic energy is owing to the interaction of the net magnetic moment $m=\sum _i{m}_i$ (see in Figure 1) with the external magnetic field H

$\begin{array}{ll}{E}_m=-m·H.& (4)\end{array}$

Magnetic moment of an individual bead, m_i, is affixed to the bead coordinate system formed by the triad {n_i, b_i, t_i}. Assuming these vectors to be the {x, y, z} axes of the local coordinate system, we characterize the local orientation of m_i by two spherical angles θ_mi and ϕ_mi, m_i = m_i{sinθ_micosϕ_mi, sinθ_misinϕ_mi, cosθ_mi}.

2.1.2. Force and Torque Balances

The elastic force $F_i^{el}$ acting on the ith bead is derived from Equations (1) and (2) as

$\begin{array}{ll}{F}_i^{el}=-\frac{\partial E_{st}}{\partial r_i}-\frac{\partial E_b}{\partial r_i}.& (5)\end{array}$

The force $F_i^{el}$ has two contributions $F_i^{el}=F_{i-1}^{el}+F_{i+1}^{el}$, where $F_{i-1}^{el}$ and $F_{i+1}^{el}$ are the elastic forces exerted on the ith bead by its neighbors, i.e., (i − 1)th and (i + 1)th beads, respectively. The corresponding torque L_i exerted on the ith bead is a sum of the magnetic $L_i^m$ and elastic $L_i^{el}$ contributions, $L_i=L_i^m+L_i^{el}$, where

$\begin{array}{ll}{L}_i^m=m_i\times H& (6)\end{array}$

and (Espanol, 1998; Wada and Netz, 2012)

$\begin{array}{l}{L}_i^{el}=-n_i\times \frac{\partial E^{el}}{\partial n_i}-b_i\times \frac{\partial E^{el}}{\partial b_i}-t_i\times \frac{\partial E^{el}}{\partial t_i}-\frac{r_{i-1,i}}{2}\times F_{i-1}^{el}\\ +\frac{r_{i,i+1}}{2}\times F_{i+1}^{el}.& (7)\end{array}$

The elastic forces and magnetic and elastic torques acting on each bead are balanced by the corresponding hydrodynamic (viscous) forces $F_i^h$ and torques $L_i^h$:

$\begin{array}{ll}{F}_i^{el}+F_i^h=0,L_i^{el}+L_i^m+L_i^h=0.& (8)\end{array}$

The viscous forces and torques applied to each bead are determined via the multipole expansion (ME) method (Filippov, 2000) that allows careful account of the mutual hydrodynamic interactions between different beads composing the filament. The method relies on multipole expansion of the Lamb's spherical harmonic solution of the Stokes equations in the unbounded fluid (Filippov, 2000). Given N spheres with radii a_i, the algorithm outputs the translational and angular velocities of the spheres for the input of prescribed forces and torques acting on each sphere. The no-slip condition at the surface of all spheres is enforced rigorously via the use of direct transformation between solid spherical harmonics centered at origins of different spheres. The ME method solves a system of ${O}(N{{L}}^2)$ linear equations for the expansion coefficients and velocities, whereas the accuracy is controlled by the truncation level, ${L}$, equal to the number of spherical harmonics retained in the expansion. The validity and accuracy of the ME algorithm have been previously tested for (i) the exact solution (in bi-spherical coordinates) for the flow past two close spheres and against (ii) a boundary element method numerical solution for the translation and rotation of straight chains of spheres (made of N = 2–30 spheres, Filippov, 2000). The method was previously applied for modeling self-locomotion of undulating filament (Berman et al., 2013), and externally driven propulsion of a microhelix (Walker et al., 2015), arc-shaped filament (Morozov et al., 2017; Sachs et al., 2018), and random fractal-like aggregates (Mirzae et al., 2018). In our calculations the truncation level was set to ${L}=2$–3, as it yielded sufficiently accurate results, while keeping the small size of the linear system to be solved at each time step. Notice that the extension of the ME method to wall-bounded domain is possible (Ozarkar and Sangani, 2008). The applied ME method valid for unbounded fluid domain can also be used for modeling propulsion near boundary by adding a stationary sphere of a large radius (e.g., larger than the filament length) in the vicinity of the micro-swimmer.

2.1.3. Non-dimensionalization

In the numerical simulations we distinguish between the parameters for the rigid and flexible segments of the filament. In particular, the elasticity constants of the rigid segments are much higher than the corresponding parameters for the flexible part. Typically the magnetic Ni segment is rigid, and thus its net magnetic moment m can be distributed equally between beads comprising it.

For simplicity we assume uniform thickness of the filament, i.e., the equal-sized beads with radius a composing the rigid and flexible parts. For characteristic scales of length, time and elasticity we use, respectively, the bead radius a, the period of the rotating field T = 2π/ω and the twisting modulus C of the flexible segment.

The non-dimensional form of the forces and torques balances in Equation (8) read,

$\begin{array}{ll}{\widehat{F}}_i^h=-\frac{1}{p}(C_{st}{\widehat{F}}_i^{st}+C_b{\widehat{F}}_i^b),& (9)\end{array}$

$\begin{array}{ll}{\widehat{L}}_i^h=-\frac{1}{p}(C_{st}{\widehat{L}}_i^{st}+C_b{\widehat{L}}_i^b+{\widehat{L}}_i^t+C_m{\widehat{L}}_i^m),& (10)\end{array}$

where the corresponding ${\widehat{F}}_i$ and ${\widehat{L}}_i$ stand for the dimensionless hydrodynamic (superscript h), bending (b), stretching (st), and magnetic (m) forces and torques, exerted on ith bead. Here p = ω/ω_f is a ratio of the frequency of the magnetic field, ω, and the characteristic frequency of the filament, ${\omega }_f=\pi C/\eta a^4$, where η is the dynamic viscosity of the fluid. C_m = 2amH/C measures the relative magnitude of the elastic and magnetic moments, $C_{st}=2K{a}^3/C$ stands for the ratio of stretching and twisting coefficients and C_b = A/C for the ratio of bending and twisting coefficients. Typically we assume that C_b ≈ 1 (Landau and Lifshitz, 1970) and C_st ≫ 1. The last condition implies nearly incompressible filament. Notice that the twist modulus C of the flexible part of the filament was chosen for non-dimensionalization. The observable deformation of the flexible filament in the external magnetic field takes place when the magnitudes of characteristic magnetic and elastic torques are close, i.e., C_m ~ 1. For the rigid segment(s) the corresponding value of C_m is smaller by a factor C/C_rigid ≪ 1.

The forces and torques exerted on a bead in Equations (9) and (10) are proportional to the rigidity of its links and to p⁻¹. Therefore, in cases where the composite nanowire has rigid segment(s) (with p ≪ 1) composed of multiple beads, the equations of time-evolution become stiff, forcing a small time step Δt and prolonged computation time. To overcome this difficulty, we assume that the rigid (small-p) section does not deform and participate in a rigid-body motion. The unknown velocities of the individual beads comprising the rigid section are expressed in terms of the translation and rotation velocity of the segment as a whole. In addition, the net hydrodynamic force and torque acting on the rigid segment are calculated as

$\begin{array}{ll}{\widehat{F}}_{rigid}={\displaystyle \sum _i}{\widehat{F}}_i^h,& (11)\end{array}$

$\begin{array}{ll}{\widehat{L}}_{rigid}={\displaystyle \sum _i}({\widehat{L}}_i^h+R_i\times {\widehat{F}}_i^h),& (12)\end{array}$

where i iterates over the beads composing the rigid part and R_i is the vector connecting a central point of the segment and the center of the ith bead. This modification removes the need to compute negligibly small deformations due to stiff elastic links. Furthermore, for every rigid segment made of M beads, the number of unknown variables in the linear system is reduced by 6(M − 1). Using this method we found acceleration of up to ~20 times in simulation speed depending on the value of p in comparison with the scheme that treats both rigid and flexible segments similarly based on Equations (9) and (10).

2.1.4. Initial Setup and Time-Evolution

Consider the laboratory coordinate system with the axes {X, Y, Z}. Initially the filament is at rest, parallel to the Z-axis, so as at t = 0, all tangent unit vectors are t_i = (0, 0, 1). We assume the initial orientation of the binormal b_i and normal n_i vectors along the axes X and Y, respectively: b_i = (1, 0, 0) and n_i = (0, 1, 0). As was mentioned above, the orientation of the magnetic moment of ith magnetized bead, ${\hat{m}}_i=m_i/m$, is prescribed by the two spherical angles θ_mi and ϕ_mi and initial orientation of the net magnetic moment is given by a superposition $\sum _i{\hat{m}}_i$.

The time-varying (e.g., rotating or oscillating) external magnetic field H is switched on at the moment t = 0. Since the initial geometry is prescribed, one can determine the forces F_i from (5) and torques L_i from (6–7) exerted on all beads. Then equating the respective elastic and hydrodynamic forces and torques one finds the translational U_i and angular Ω_i velocities of beads at first time step by using the ME algorithm.

Using these values we obtain the new positions and orientations of beads at the moment t = Δt:

$\begin{array}{ll}{r}_i(\Delta t)=r_i(0)+U_i\Delta t,& (13)\end{array}$

$\begin{array}{ll}{\xi }_i(\Delta t)={\xi }_i(0)+\left[{\Omega }_i\times {\xi }_i(0)\right]\Delta t,& (14)\end{array}$

whereas ξ_i stands for either b_i, n_i or t_i. Since the magnetic moment of all magnetic beads is affixed to the triad {b_i, n_i, t_i}, its time-evolution is also governed by an equation similar to (14) whereas ξ_i is replaced by ${\hat{m}}_i$ everywhere. The updated vectors according to Equation (14) are normalized with |ξ_i| to keep their magnitude of one. At the next time step the procedure is repeated.

2.2. Validation

We tested the validity of the proposed numerical model by comparing its results to two benchmark problems for which theoretical predictions exist: (i) magnetized rigid cylinder actuated by a rotating magnetic field; (ii) the twirling-whirling instability of an elastic filament rotated by its end. The first problem tests the hydrodynamic part of the algorithm, as filament is assumed to be rigid, while the second problem concerns the interplay between elasticity and hydrodynamics.

2.2.1. Magnetized Rigid Cylinder in a Rotating Field

Consider the dynamics of a cylinder with magnetic moment m driven by a rotating uniform magnetic field H = H(cosωt, sinωt, 0): at low frequency ω the cylinder tumbles in the plane of the field rotation in-sync with the field; at a certain critical frequency, ${\omega }_c^{(I)}$ (or the corresponding dimensionless parameter $p_c^{(I)}$), the tumbling switches to in-sync wobbling, where the precession angle (i.e., the angle between the field rotation axis and the long axis of the cylinder) gradually diminishes as the frequency is increased; then, at the step-out frequency ${\omega }_c^{(II)}$ ($p_c^{(II)}$) the synchronous regime switches to the asynchronous one. The analytical solution of the problem was given in Morozov and Leshansky (2014). The dynamics of a cylinder-like structure of a linear chain composed of 5 beads (see Supplementary Material for details) was simulated, yielding the numerical value of the angle α at steady-state between the vectors m and H. The ratio of the longitudinal and transverse components of the cylinder's magnetization m was set to m_||/m_⊥ = 2.33, and the dimensionless parameters were set to C_st = 2, C_b = 1, and C_m = 0.0015, resulting in the critical values $p_c^{(I)}\approx 2.48·10^{-7}$ and $p_c^{(II)}\approx 8.70·10^{-7}$. Table 1 shows a comparison of the theoretical values of α in the tumbling and wobbling regimes given, respectively by Equations (16) and (18) in Morozov and Leshansky (2014), to the numerical values obtained in this work. The agreement between the theoretical predictions and the numerical results is excellent.

TABLE 1

Table 1. The angle α between the magnetic moment m of a linear chain made of 5 beads and the rotating magnetic field H in the in-synch tumbling (upper table) and wobbling (lower table) regimes.

2.2.2. The Twirling-Whirling Instability

Consider an elastic rod of length L in a viscous fluid that at one end is forced to rotate about its long axis with frequency ω₀ and the other end is free. Upon increasing the driving frequency the twirling regime of the rod switches to the whirling regime and the initially straight rod buckles. This elasto-hydrodynamic instability occurs at the critical frequency (Wolgemuth et al., 2000; Powers, 2010)

$\begin{array}{ll}{\omega }_c\approx 8.98\frac{A}{{\zeta }_r{L}^2},& (15)\end{array}$

where ζ_r is the rotational drag coefficient (torque per unit length). At low angular frequency ω₀ < ω_c, the steady-state twist distribution (i.e., the twist angle per unit length of the filament) along the straight filament of length L reads (Wolgemuth et al., 2000; Powers, 2010),

$\begin{array}{ll}\tau (s)=\frac{{\zeta }_r{\omega }_0}{C}(s-L),& (16)\end{array}$

where s is the arc-length parameter.

We simulated the twirling-whirling instability numerically using bead-spring model as illustrated in Figure 2. Substituting p into Equation (16), the expression $\partial \tau /\partial s={\zeta }_r\pi p/\eta a^4$ is obtained. In the calculation the value of ${\zeta }_r\approx 10\eta a^2$ was used (Morozov and Leshansky, 2014). Table 2A shows an excellent agreement of the numerically determined values of the angle (∂τ/∂s)L² for a 25-bead filament at several driving frequencies ω₀ < ω_c to the analytical prediction in Equation (16). Substituting the critical frequency ω_c into (16) and integrating over the filament length leads to the critical twist angle between its ends, Δφ_c = 4.49C_b. The instability onset can be detected either by a sudden take-off of the bending energy or the maximum of the torsional energy. In the simulations p was gradually incremented until the torsional buckling (or whirling) of the filament occurs. Table 2B shows the critical value Δφ_c of the twist angle between filament ends at which the twirling-to-whirling transition occurs for a filament made of 25 beads, as calculated by Equation (15) and by using the bead-spring model upon varying C_b for a fixed value of C_st = 20. The agreement between the numerical results and the theoretical prediction is quite accurate, given the discrete nature of the bead-spring model. Figure 2B illustrates the near-critical shape of the filament at ω ≈ ω_c, while Figure 2C depicts the steady-state shape of the whirling filament at ω₀ > ω_c (see Supplementary Video 1).

FIGURE 2

Figure 2. The twirling-whirling instability: an initially straight flexible 25-bead filament with C_b = 0.5 immersed in a viscous liquid and rotated at constant rate ω₀ at one end, while the other end is free (A); the filament buckles at the critical frequency ω₀ = ω_c (B) acquiring the steady-state shape at higher frequency, ω₀ > ω_c (C). The color arrows stand for the binormal and normal vectors, b_i and n_i.

TABLE 2A

Table 2A. The steady-state twist distribution along the axis of a 25-bead filament at low driving frequencies ω₀ < ω_c.

TABLE 2B

Table 2B. The critical twist angle between the filament ends, Δφ_c, at the twirling-to-whirling transition.

2.3. The Two-Segment Nanowire Propeller

Pak et al. (2011) reported an experimental design of a flexible nanomotor that displays high propulsion speed. This nanowire propeller, composed of a 1.8 μm-long rigid magnetic head (1.5 μm nickel segment and 0.3 μm gold segment) with a diameter of 200 nm and a 4 μm-long flexible silver tail with a diameter of ~100 nm, was actuated by the conically rotating magnetic field H = H₀(hcosωt, −hsinωt, 1) where h = H₁/H₀ is the ratio of the rotating and constant components of the magnetic field whereas tan⁻¹h is the cone angle.

Experiments were performed to determine the dependence of the nanomotor propulsion speed on the sperm number (Pak et al., 2011),

$\begin{array}{ll}\text{Sp}=L{\left(\frac{f_{\perp }\omega }{A}\right)}^{1/4},& (17)\end{array}$

for several values of the field cone angle upon varying h. Here L is the length of the nanowire's flexible part and f_⊥ is the normal viscous drag coefficient (force density), given approximately for a slender filament by $f_{\perp }=\frac{4\pi \eta }{log(2/ϵ)+1/2}$, where ϵ = 2a/L ≪ 1 is the filament aspect ratio. In Pak et al. (2011) the value of the elastic constant A ≈ 3.6 × 10⁻²⁴ N·m² was best fitted to match the prediction of the approximate geometric theory assuming that the freely suspended nanowire can be approximated by an elastic filament with its one (clamped) end undergoing rotation with the cone angle equal to tan⁻¹h.

Using the method described in this work we simulated the nanowire propulsion assuming an initially straight filament, composed of a 5-bead rigid head and a 11-bead flexible tail (i.e., a head-to-tail ratio being similar to the experiment in Pak et al., 2011), with the magnetic moment of the head set along its longitudinal axis. Notice that simulations of very slender filaments (with ϵ ≈ 0.025 as in Pak et al., 2011) would require large number of beads and are computationally expensive. As one may expect for viscous hydrodynamics, the larger dimension (i.e., the length L) of the slender nanowire is important, while its diameter (or the aspect ratio ϵ) has a minor effect on the dynamics (as can also be seen from the definition of Sp in Equation 17). Therefore, we expect that using of a chubby nanowire with ϵ ≈ 0.09 would yet closely approximate the dynamics of a slender filament in experiments of Pak et al. (2011).

The dimensionless parameters were set to C_m = C_b = 1, C_st = 10 with p varying between 6·10⁻⁷ and 1.3·10⁻⁴. Figure 3 shows the dimensionless velocity, normalized by the length of the flexible tail, L, as a function of the sperm number Sp in (17) that can be rewritten in terms of dimensionless parameters as

$\begin{array}{ll}\text{Sp}=2{ϵ}^{-1}{\left(\frac{\pi f_{\perp }p}{\eta C_b}\right)}^{1/4}.& (18)\end{array}$

FIGURE 3

Figure 3. Dependence of the dimensionless swimming speed U/ωL of the 2-segment magnetic nanowire on the sperm number, Sp (upper panel). The experimental data (symbols), taken from Pak et al. (2011), and the numerical results for the bead-spring model (solid lines), are compared for three different actuating magnetic fields with h⁻¹ = 1.43 (blue); h⁻¹ = 1.18 (red); and h⁻¹ = 0.70 (black). The lower panel shows the steady shapes of the simulated nanowire made of an 11-bead flexible tail (Ag, gray) and a 5-bead rigid magnetic (Ni, red) head magnetized longitudinally for h⁻¹ = 0.70 and Sp = 1.1 (A), Sp = 1.5 (B), and Sp = 3.3 (C).

Best fitting of the numerical predictions to the experimental results requires ≈2.3 times larger value of the bending modulus, A ≈ 8.3 × 10⁻²⁴ N·m², when compared to the estimate of Pak et al. (2011). It can be readily seen that the agreement between the numerical predictions and the experimental findings in Figure 3 is excellent without any extra adjustable parameters¹. In particular, the maximum speed is attained at Sp ≈ 1.55 (see Supplementary Videos 2–4). The numerical method allows for finding the optimal value of h⁻¹ maximizing the propeller's velocity. Computing the velocity of the propeller at the optimal sperm number (Sp ≈ 1.55) reveals that the optimal value of h⁻¹ is ≈ 0.5, meaning the magnitude of the constant magnetic field is half the magnitude of the rotating field (see Figure 4).

FIGURE 4

Figure 4. The dependence of dimensionless velocity U/ωL on h⁻¹ at the (near) optimal value of Sp = 1.55. The data is based on the simulation of a nanowire propeller in Figure 3.

2.4. The Three-Segment Nanowire Propeller

Gao et al. (2010) demonstrated the propulsion of a nanowire motor with magnetic nickel head and golden tail (both rigid segments) connected by a flexible link made of porous silver. In the experiments they found that upon applying an external rotating magnetic field, two propulsion gaits are observed: slow backward motion (i.e., toward the passive tail), and fast forward motion (toward the magnetic head). We simulated this nanowire as a filament made of three segments each composed of 5, 10, and 5 beads, corresponding to the magnetic head, flexible link and passive tail, respectively. All simulations of the initially linear filaments resulted in the forward propulsion gait regardless of the magnetic moment orientation. However, it was found that some imperfection in the initial shape of the filament (e.g., slight intrinsic curvature of the Ni-head) the bidirectional propulsion gait can be realized in a qualitative agreement with the experiments by Gao et al. (2010), where the microfabricated nanowires at rest were not straight.

Morozov et al. (2017) predicted and it was later demonstrated experimentally (Sachs et al., 2018) that rigid planar objects (such as arc) can change their propulsion direction when the sign of the off-plane magnetization component is reversed. Similarly, it was observed in the simulation of the three-segment nanowire that reversal the off-plane component of the head's magnetic moment, results in propulsion reversal, i.e., switches between fast head-forward and slow tail-forward motion (see videos #5 and #6 in the Supporting Materials). In both cases the magnetic head follows the rotation of the actuating field. Figures 5A,B depict the nanowire steady-state shapes for both propulsion gaits. Setting the dimensionless parameters to C_m = C_b = 1, C_st = 10, and p = 10⁻⁴, the forward and backward dimensionless speeds U/Lω (here L is the length of the flexible part) were determined numerically to be 0.0011 and 0.0084, respectively, implying that the head-forward motion is about 8 times faster. These findings are in a qualitative agreement with the experimental observations of Gao et al. (2010).

FIGURE 5

Figure 5. The simulated three-segment nanowire, composed of Au rigid tail (yellow), a flexible middle Ag link (gray) and a slightly curved rigid magnetic Ni head (red). Two realizations of the nanowire were simulated, both have the same magnetization component in the swimmer's plane, whereas their off-plane components of magnetization have the same magnitude, but opposite sign. (A) The steady-state 3D shape of the nanowire when propelling head-forward (fast); (B) The steady-state shape of the nanowire when propelling tail-forward (slow).

3. Conclusions

We developed a numerical scheme based on a discrete bead-spring model for simulating soft nanowire-based propellers. The algorithm was favorably tested by simulating a well-known twirling-whirling instability of an elastic filament rotated in a viscous liquid. The comparison of the numerical prediction for the propulsion speed vs. the reported experimental results for the Ni-Ag nanowire actuated by a conically rotating magnetic field in Pak et al. (2011) as a function of actuation frequency (Sp) upon varying the ratio of the in-plane rotating and constant field components (h⁻¹), showed an excellent agreement (see Figure 3). At moderate actuation frequencies the Ni-Ag nanowire adopts the arc shape, while it develops 3D geometric chirality gradually as Sp increases above ≈ 2.5 The simulation results predict an optimal cone angle of the actuating field, ≈ 63°, that maximizes the value of the non-dimensional propulsion velocity U/ωL at Sp ≈ 1.5. Surprisingly, at the optimum the nanowire has a shape of a planar arc (see Figure 3B), indicating that efficient propulsion of magnetically driven soft nanowires does not require 3D geometric chirality as was previously suggested. This finding is in accord with Mirzae et al. (2018) where an efficient steering of rigid magnetic arc-shaped nanomotors was demonstrated.

Using our model we could substantiate important and previously unexplained details, such as bidirectional (fast head-forward and slow tail-forward) propulsion of the Ni-Ag-Au nanowire of Gao et al. (2010) powered by in-plane rotating magnetic field due to potential shape imperfection which is probably inevitable in nanowire nanofabrication. The proposed numerical model can be used as an efficient tool for design and optimization of previously proposed and also novel soft nanomachines. It can be also extended to modeling of polarizable (superparamagnetic) soft nanowires, biological and biohybrid microswimmers and other externally or self-propelled machines.

Data Availability Statement

The original contributions generated for the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding author/s.

Author Contributions

YM implemented the bead-spring algorithm and performed all the computations. BR and KM assisted in developing the key ideas, research design, and implementation. AL directed the research and wrote the paper. All authors contributed to the article and approved the submitted version.

Funding

This work was supported in part by the Israel Science Foundation (ISF) via the grant No. 1744/17 (AML) and by a joint grant from the Center for Absorption in Science of the Ministry of Immigrant Absorption and the Committee for Planning and Budgeting of the Council for Higher Education under the framework of the KAMEA Program (KIM).

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/frobt.2020.595777/full#supplementary-material

Supplementary Video 1. Twirling-whirling instability of an elastic filament for C_b = 0.5.

Supplementary Video 2. The two-segment (Ni-Ag) nanowire propeller at Sp = 1.1.

Supplementary Video 3. The two-segment (Ni-Ag) nanowire propeller at Sp = 1.5.

Supplementary Video 4. The two-segment (Ni-Ag) nanowire propeller at Sp = 3.3.

Supplementary Video 5. Head-forward (fast) motion of the three-segment (Ni-Ag-Au) nanowire propeller.

Supplementary Video 6. Tail-forward (slow) motion of the three-segment (Ni-Ag-Au) nanowire propeller.

Footnote

1. ^In Pak et al. (2011) the results shown by black diamonds in Figure 3 correspond to h⁻¹ = 0.95, while our numerical simulations suggest that the correct value is actually close to h⁻¹ = 0.7.

References

Alouges, F., DeSimone, A., Giraldi, L., Or, Y., and Wiezel, O. (2019). Energy-optimal strokes for multi-link microswimmers: Purcell's loops and taylor's waves reconciled. New J. Phys. 21:043050. doi: 10.1088/1367-2630/ab1142

CrossRef Full Text | Google Scholar

Berman, R. S., Kenneth, O., Sznitman, J., and Leshansky, A. M. (2013). Undulatory locomotion of finite filaments: lessons from caenorhabditis elegans. New J. Phys. 15:145. doi: 10.1088/1367-2630/15/7/075022

CrossRef Full Text | Google Scholar

Cheang, U. K., Meshkati, F., Kim, D., Kim, M. J., and Fu, H. C. (2014). Minimal geometric requirements for micropropulsion via magnetic rotation. Phys. Rev. E 90:033007. doi: 10.1103/PhysRevE.90.033007

PubMed Abstract | CrossRef Full Text | Google Scholar

Dreyfus, R., Baudry, J., Roper, M. L., Fermigier, N., Stone, H. A., and Bibette, J. (2005). Microscopic artificial swimmers. Nature 437, 862–865. doi: 10.1038/nature04090

PubMed Abstract | CrossRef Full Text | Google Scholar

Español, P. (1998). Fluid particle model. Phys. Rev. E 57, 2930–2948. doi: 10.1103/PhysRevE.57.2930

CrossRef Full Text | Google Scholar

Filippov, A. V. (2000). Drag and torque on cluster of N arbitrary spheres at low Reynolds number. J. Colloid Interface Sci. 229, 184–195. doi: 10.1006/jcis.2000.6981

PubMed Abstract | CrossRef Full Text | Google Scholar

Gao, W., Feng, X., Pei, A., Kane, C. R., Tam, R., Hennessy, C., et al. (2014). Bioinspired helical microswimmers based on vascular plants. Nano Lett. 14, 305–310. doi: 10.1021/nl404044d

PubMed Abstract | CrossRef Full Text | Google Scholar

Gao, W., Sattayasamitsathit, S., Manesh, K. M., Weihs, D., and Wang, J. (2010). Magnetically powered flexible metal nanowire motors. J. Am. Chem. Soc. 132, 14403–14405. doi: 10.1021/ja1072349

PubMed Abstract | CrossRef Full Text | Google Scholar

Gauger, E., and Stark, H. (2006). Numerical study of a microscopic artificial swimmer. Phys. Rev. E 74:021907. doi: 10.1103/PhysRevE.74.021907

PubMed Abstract | CrossRef Full Text | Google Scholar

Ghosh, A., and Fischer, P. (2009). Controlled propulsion of artificial magnetic nanostructured propellers. Nano Lett. 9, 2243–2245. doi: 10.1021/nl900186w

PubMed Abstract | CrossRef Full Text | Google Scholar

Gutman, E., and Or, Y. (2015). Simple model of a planar undulating magnetic microswimmer. Phys. Rev. E 90:013012. doi: 10.1103/PhysRevE.90.013012

PubMed Abstract | CrossRef Full Text | Google Scholar

Jang, B., Gutman, E., Stucki, N., Seitz, B. F., Wendel-Garcia, P. D., Newton, T., et al. (2015). Undulatory locomotion of magnetic multilink nanoswimmers. Nano Lett. 15, 4829–4833. doi: 10.1021/acs.nanolett.5b01981

PubMed Abstract | CrossRef Full Text | Google Scholar

Landau, L. D., and Lifshitz, E. M. (1970). Theory of Elasticity. Oxford: Pergamon Press.

PubMed Abstract

Li, T., Li, J., Morozov, K. I., Wu, Z., Xu, T., Rozen, I., et al. (2017). Highly efficient freestyle magnetic nanoswimmer. Nano Lett. 17, 5092–5098. doi: 10.1021/acs.nanolett.7b02383

PubMed Abstract | CrossRef Full Text | Google Scholar

Maggs, A. (2000). Twist and writhe dynamics of stiff polymers. Phys. Rev. Lett. 85, 5472–5475. doi: 10.1103/PhysRevLett.85.5472

PubMed Abstract | CrossRef Full Text | Google Scholar

Mirzae, Y., Dubrovki, O., Kenneth, O., Morozov, K. I., and Leshansky, A. M. (2018). Geometric constraints and optimization in externally driven propulsion. Sci. Robot. 3:eaas8713. doi: 10.1126/scirobotics.aas8713

CrossRef Full Text | Google Scholar

Morozov, K. I., and Leshansky, A. M. (2014). The chiral magnetic nanomotors. Nanoscale 6, 1580–1588. doi: 10.1039/C3NR04853E

CrossRef Full Text | Google Scholar

Morozov, K. I., Mirzae, Y., Kenneth, O., and Leshansky, A. M. (2017). Dynamics of arbitrary shaped propellers driven by a rotating magnetic field. Phys. Rev. Fluids 2:044202. doi: 10.1103/PhysRevFluids.2.044202

CrossRef Full Text | Google Scholar

Ozarkar, S. S., and Sangani, A. S. (2008). A method for determining Stokes flow around particles near a wall or in a thin film bounded by a wall and a gas-liquid interface. Phys. Fluids 20:063301. doi: 10.1063/1.2930667

CrossRef Full Text | Google Scholar

Pak, O. S., Gao, W., Wang, J., and Lauga, E. (2011). High-speed propulsion of flexible nanowire motors: theory and experiments. Soft Matter 7, 8169–8181. doi: 10.1039/c1sm05503h

CrossRef Full Text | Google Scholar

Peters, C., Ergeneman, O., Nelson, B. J., and Hierold, C. (2013). “Superparamagnetic swimming microrobots with adjusted magnetic anisotropy,” in 2013 IEEE 26th International Conference on Micro Electro Mechanical Systems (MEMS) (IEEE), 564–567. doi: 10.1109/MEMSYS.2013.6474304

CrossRef Full Text | Google Scholar

Powers, T. R. (2010). Dynamics of filaments and membranes in a viscous fluid. Rev. Mod. Phys 82, 1607–1631. doi: 10.1103/RevModPhys.82.1607

CrossRef Full Text | Google Scholar

Sachs, J., Morozov, K. I., Kenneth, O., Qiu, T., Segreto, N., Fischer, P., et al. (2018). Role of symmetry in driven propulsion at low Reynolds number. Phys. Rev. E 98:063105. doi: 10.1103/PhysRevE.98.063105

CrossRef Full Text | Google Scholar

Schamel, D., Mark, A. G., Gibbs, J. G., Miksch, C., Morozov, K. I., Leshansky, A. M., et al. (2014). Nanopropellers and their actuation in complex viscoelastic media. ACS Nano 8, 8794–8801. doi: 10.1021/nn502360t

PubMed Abstract | CrossRef Full Text | Google Scholar

Smith, E. J., Makarov, D., Sanchez, S., Fomin, V. M., and Schmidt, O. G. (2011). Magnetic microhelix coil structures. Phys. Rev. Lett. 107:097204. doi: 10.1103/PhysRevLett.107.097204

CrossRef Full Text | Google Scholar

Tottori, S., and Nelson, B. J. (2018). Controlled propulsion of two-dimensional microswimmers in a precessing magnetic field. Small 14:1800722. doi: 10.1002/smll.201800722

PubMed Abstract | CrossRef Full Text | Google Scholar

Tottori, S., Zhang, L., Qiu, F., Krawczyk, K. K., Franco-Obregon, A., and Nelson, B. J. (2012). Magnetic helical micromachines: fabrication, controlled swimming, and cargo transport. Adv. Mater. 22, 811–816. doi: 10.1002/adma.201103818

CrossRef Full Text | Google Scholar

Vach, P. J., Brun, N., Bennet, M., Bertinetti, L., Widdrat, M., Baumgartner, J., et al. (2013). Selecting for function: solution synthesis of magnetic nanopropellers. Nano Lett. 13, 5373–5378. doi: 10.1021/nl402897x

PubMed Abstract | CrossRef Full Text | Google Scholar

Vach, P. J., Fratzl, P., Klumpp, S., and Faivre, D. (2015). Fast magnetic micropropellers with random shapes. Nano Lett. 15, 7064–7070. doi: 10.1021/acs.nanolett.5b03131

CrossRef Full Text | Google Scholar

Wada, H., and Netz, R. R. (2012). “Actively twisted polymers and filaments in biology,” in Non-equilibrium Soft Matter Physics, eds S. Komura and T. Ohta (Berlin: World Scientific Publishing Company), 275–363. doi: 10.1142/9789814360630_0006

CrossRef Full Text | Google Scholar

Walker, D., Kübler, M., Morozov, K. I., Fischer, P., and Leshansky, A. M. (2015). Optimal length of low Reynolds number nanopropellers. Nano Lett. 15, 4412–4416. doi: 10.1021/acs.nanolett.5b01925

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, J. (2013). Nanomachnines: Fundamentals and Applications. Weinheim: Wiley-VCH.

Google Scholar

Wolgemuth, C. W., Powers, T. R., and Goldstein, R. E. (2000). Twirling and whirling: viscous dynamics of rotating elastic filaments. Phys. Rev. Lett. 84, 1623–1626. doi: 10.1103/PhysRevLett.84.1623

PubMed Abstract | CrossRef Full Text | Google Scholar

Yu, Y., Shang, L., Gao, W., Zhao, Z., Wang, H., and Zhao, Y. (2017). Microfluidic lithography of bioinspired helical micromotors. Angew. Chem. Int. Ed. 56, 12127–12131. doi: 10.1002/anie.201705667

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhang, L., Abbott, J. J., Dong, L., Kratochvil, B. E., Bell, D., and Nelson, B. J. (2009). Artificial bacterial flagella: fabrication and magnetic control. Appl. Phys. Lett. 94:064107. doi: 10.1128/JB.00587-09

CrossRef Full Text | Google Scholar