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.

There are three contributions into the elastic energy: stretching, bending and torsion, Eel=Est+Eb+Et. The stretching energy reads

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 Aa(1-ti·r^i,i+1), where A is the bending modulus and r^i,i+1 is the unit vector along the line of centers of i and i + 1 beads, r^i,i+1=ri-ri+1|ri-ri+1|. Summing up all bead-spring contacts and symmetrizing the bending energy relative to both filament ends, one obtains

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 Eb=(A/2a)∑i(1-ti·ti+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 r^i,i+1, i.e., 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 ti·ti+1≅cosθicosθi+1≈1-12θi2-12θi+12. The expression in the square brackets in Equation (2) reads 2-(ti+ti+1)·r^i,i+1=1-cosθi-cosθi+1≈12θi2+12θi+12. Thus, in the limit of small deformations we have 2-(ti+ti+1)·r^i,i+1≈1-ti·ti+1.

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

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=∑imi (see in Figure 1) with the external magnetic field H

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}.

The elastic force Fiel acting on the ith bead is derived from Equations (1) and (2) as

The force Fiel has two contributions Fiel=Fi-1el+Fi+1el, where Fi-1el and Fi+1el 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 Lim and elastic Liel contributions, Li=Lim+Liel, where

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

The elastic forces and magnetic and elastic torques acting on each bead are balanced by the corresponding hydrodynamic (viscous) forces Fih and torques Lih:

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(NL2) 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.

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,

where the corresponding F^i and 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, ωf=πC/ηa4, where η is the dynamic viscosity of the fluid. C_m = 2amH/C measures the relative magnitude of the elastic and magnetic moments, Cst=2Ka3/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

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).

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, m^i=mi/m, is prescribed by the two spherical angles θ_mi and ϕ_mi and initial orientation of the net magnetic moment is given by a superposition ∑im^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:

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 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.

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, ωc(I) (or the corresponding dimensionless parameter pc(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 ωc(II) (pc(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 pc(I)≈2.48·10-7 and pc(II)≈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.

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)

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),

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 ∂τ/∂s=ζrπp/ηa4 is obtained. In the calculation the value of ζr≈10ηa2 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).