Stability of Soft Magnetic Helical Microrobots

: Nano/microrobotic swimmers have many possible biomedical applications such as drug delivery and micro-manipulation. This paper examines one of the most promising classes of these: rigid magnetic microrobots that are propelled through bulk ﬂuid by rotation induced by a rotating magnetic ﬁeld. Propulsion corresponds to steadily rotating and translating solutions of the dynamics of such microrobots that co-rotate with the magnetic ﬁeld. To be observed in experiments and be amenable to steering control, such solutions must also be stable to perturbations. In this paper, we analytically derive a criterion for the stability of such steadily rotating solutions for a microrobot made of soft magnetic materials, which have a magnetization that depends on the applied ﬁeld. This result generalizes previous stability criteria we obtained for microrobots with a permanent magnetization.

A common group of such microrobots are rigid with nanoscale features that lead to chiral helical geometries [22,23], and are propelled by exerting a magnetic torque on the body through a rotating magnetic field. Similar to a corkscrew, such a torque rotating a helical tail can result in a translational motion along the rotation axis of the motor. Typically, they are made of a helical tail that is either magnetic itself or attached to a magnetic head. Such helical microrobots take inspiration from biological microorganisms such as bacteria which swim by rotating a helical flagellar filament [24] and are particularly promising for possible biomedical applications [25]. Over the past decade, several studies have investigated the dynamics and control of such helical shaped magnetic swimmers experimentally and theoretically [25][26][27][28][29][30][31][32][33][34][35].
In experiments, it is most common for these magnetic microrobots to be actuated by a uniform magnetic field h that rotates in a plane. The rotation of the field is described by its angular velocity ω, which is perpendicular to the plane of rotation ( Figure 1). Steady propulsion is achieved when the microrobot co-rotates with the field, meaning that the microrobot's angular velocity Ω is equal to the field angular velocity, Such co-rotating solutions exist for small enough rotation frequencies such that the magnetic field is able to exert enough torque on the microrobot to rotate it at the rotation frequency. In the limit of very small rotation frequencies, the magnetization of the microrobot follows the direction of the field, while, if the rotation frequency becomes too large, the microrobot can no longer rotate fast enough to keep up with the field, leading to unsteady rotation dynamics [32,36]. When steady co-rotating solutions exist, the average propulsion direction is in the direction of the rotation axis, although the instantaneous translational velocity of the microrobot U may be in some other direction [32]. For a steady solution in the presence of a rotating applied magnetic field h, the soft magnetic head is magnetized with a magnetization vector m, resulting in a torque that rotates the body together with the field, with constant angular velocity Ω in the body frame. The microswimmer is propelled with constant translational velocity U in the body frame.
To achieve controllable propulsion, however, it is not enough for a steady co-rotating solution to exist; it also must be stable to perturbations from the control system and environment, for example from thermal Brownian motion. An unstable co-rotating solution would never be observed experimentally as any perturbation would result in the microrobot departing from the steady rotational dynamics. The region of stability, i.e., the size of perturbations to which a steady solution is stable, also sets limitations on steering, since steering requires that the rotation axis of the field be changed slowly enough that the microrobot maintains motion co-rotating about the new rotation axis. Therefore, previous works [28,32] have examined the stability of such microswimmers with a permanent (constant in time) magnetization both experimentally and through models. However, many such rotating magnetic swimmers are made of soft ferromagnetic materials [23,37], which have magnetizations that depend on the applied field. It is known that magnetization direction of permanent magnets affects their propulsion characteristics [35,38], thus we expect that in general the difference in magnetization behavior of different magnetic materials will result in different responses to applied magnetic fields. Soft magnets are especially desirable in the design of microrobots due to their availability and compatibility with manufacturing techniques. In addition, unlike hard magnets, soft magnets do not require initial magnetization. Therefore, there is a need to evaluate the stability of such soft magnetic microswimmers. In this paper, we analytically evaluate the stability of a steadily rotating magnetic microswimmer based on its geometry and magnetic properties, eliminating the need to numerically propagate the model in time to see its long-term stability behavior.

Model of Soft Magnetism
To model the response of a swimmer with an ellipsoidal soft magnetic head, we use the formulation of soft magnetism for an ellipsoid developed and experimentally validated by Abbott et al. [39]. Briefly, for a soft magnet, the magnetization m is approximated as a linear function of the applied magnetic field h at small field strengths. As the applied field increases, the magnetization also increases, but only up to the saturation magnetization magnitude m s . For larger fields, the magnetization has magnitude m s but its direction gradually approaches that of the applied field.
In equations, in the linear regime, m = χh, (2) where χ is a diagonal matrix with diagonal elements (1/n a , 1/n r , 1/n r ), where n r and n a are demagnetization factors of the ellipsoid that are functions of its geometry; for a prolate ellipsoid with ratio of major to minor axis r (r ≥ 1), n a = 1 − 1 and n r = (1 − n a )/2 [39]. If the magnetic field amplitude h is large enough that the magnitude of the magnetization given by this equation is larger than m s , then the linear regime no longer applies and instead the magnetization is saturated. In this saturated regime, the magnetization has magnitude m s , and its direction is specified by two angles (Figure 2), the polar angle θ m from the major axis of the ellipsoid and the azimuthal angle φ m about the major axis, m = m s (cos θ m , sin θ m cos φ m , sin θ m sin φ m ) in Cartesian components. The applied field is likewise described by its magnitude h and two angles θ h and φ h specifying its direction relative to the major axis of the ellipsoid, Due to symmetry, the magnetization always lies in the plane formed by the applied field and the major axis of the ellipsoid, and therefore The polar angle of the magnetization θ m depends on the magnitude and direction of the applied field through The magnetization in the saturated regime given by Equation (6) has the property that for very large field strengths (h → ∞), θ h = θ m , i.e., the magnetization points along the field direction.
Once the magnetization of the ellipsoid is known, the torque τ exerted on the microrobot with magnetization m by the field is given by [32] where µ 0 is the magnetic permeability of free space and V is the volume of magnetic material in the microrobot with average magnetization m [32]. In a uniform magnetic field, the force on a magnetic dipole is zero.

Steady Solutions
These microrobots swim in a low Reynolds number regime, where the effect of inertial forces on the flow is negligible compared to the viscous forces. In such low Reynolds number regime, the translational velocity U and angular velocity Ω of the swimmer ( Figure 1) have a linear relation with the external force F and torque τ acting on the swimmer through a mobility matrix, written in terms of 3 × 3 submatrices K, C, and M. The linear and angular velocities are related to force and torque by K and M, respectively, while the linear velocity is related to the torque and vice versa by C and its transpose. The explicit forms of the submatrices used for our example swimmer in Section 3.4 are calculated using the method of regularized Stokeslets [40,41], as detailed in the next subsection.
To find steady co-rotating solutions, we follow the methods of our previous work [32]. We work in the body basis of the swimmer, so that, if a co-rotating swimmer satisfies Equation (15), then the fact that the magnetic field rotates in a plane perpendicular to its rotation axis implies the constraint For a field directionĥ = (cos θ h , sin θ h cos φ h , sin θ h sin φ h ) in the body basis, since φ h = φ m , the torque τ has no x component, i.e., it lies in the y − z plane. Therefore,τ = (0, − sin φ h , cos φ h ). Plugging these forms into the constraint in Equation (9) and using the mobility matrix in Equation (8) yields an equation relating θ h and φ h . Each pair of (θ h , φ h ) that satisfies this constraint gives a field direction that specifies a steady solution. Once the field direction is known, assuming an experimentally prescribed field magnitude h, we can calculate the magnetization using the model of the previous section, the torque from Equation (7), and then the swimmers angular (Ω) and translational (U) velocities from Equation (8).
Thus, for each steady solution, we know its corresponding magnetic field in the body basis specified by either Cartesian components or spherical parameters (θ h , φ h , h), as well as its magnetization in the body basis specified by either Cartesian components or spherical parameters (θ m , φ m , m). The main task of this paper is to determine whether such a steady solution is dynamically stable or unstable.

Calculation of Mobility Matrices
The method of regularized Stokeslets [40,41] represents the flow around an object as a superposition of flows due to spread-out point forces distributed on the surface area of the object. Our implementation of this method was described previously in [42,43]. In particular, our calculation of mobility matrices exactly follows the description in [43].
Briefly, rather than a singular point force, we use a force distribution which spreads the force out into a "blob". The precise form of the blob we use is given by Equation (9) of [41], φ = 15 4 /[8π(r 2 + 2 ) 7/2 ]. The flow velocity field v at position x due to a set of N of these force blobs is where the f α are unknown forces of each blob at position x α and the tensor S is the regularized Stokeslet, the precise form of which is Equation (10) of [41]. The forces are determined by requiring that the flow at the surface satisfies the no-slip condition at the blob locations, i.e., where the v α are the velocities of surface at the blob locations. If we prescribe rigid-body motion of the object with translational velocity U and angular velocity Ω, then the v α are known, yielding a linear set of equations that can be solved for the blob forces f α . Then, the total force and torque on the rigid object are Thus, for any translational U and angular velocity Ω, we calculate the total force and torque on a rigid object; the relation between the two sets of quantities is linear, which specifies the mobility matrix in Equation (8).
To evaluate the mobility matrix for the example microrobot of Section 3.4, we discretize the surface of the geometry with regularized Stokeslets as follows. The surface of the ellipsoidal head is discretized by placing uniformly separated Stokeslets around the perimeter of cross-sections of the ellipsoid perpendicular to the major axis. The regularization parameter of regularized Stokeslets on the ellipsoid surface is equal to the separation between successive cross-sections. The surface of the helical tail is discretized by placing uniformly separated Stokeslets around the perimeter of cross-sections perpendicular to the centerline of the helix filament. The regularization parameter of the regularized Stokeslets on the tail surface is equal to the separation between successive cross-sections. The Stokeslets on successive cross-sections are staggered to be more uniformly distributed, as described in [43]. We use the number of Stokeslets on each cross-section such that their separation is as close as possible to the separation between successive cross-sections. We choose the separation between successive cross-sections on the ellipsoid and tail to be equal.
Convergence is tested by evaluating the percentage change in ||K||, ||C||, and ||M||, where, e.g., ||K|| = ∑ ij K ij K ij , as the number of Stokeslets is increased by decreasing the separation between successive cross-sections. We use values for the mobility matrix obtained when the change between mobility matrices for the finest three discretizations is less than 1%.

General Stability Criterion
When the microswimmer is rotating steadily with the field, the field is constant in the body basis and the steady angular velocity in the body basis Ω s is where the mobility matrices, magnetization, and field on the right-hand side are also referenced in the body basis. Assume that at time t = t 0 the microswimmer is perturbed and hence rotated by an infinitesimal rotation vector σ(t 0 ). Over time, the perturbed body will rotate, but at a different angular velocity from the steady solution, thus the rotation vector connecting the two will evolve in time as σ(t). We reference σ to the frame of the steadily rotating body; thus, if a perturbed orientation also steadily co-rotates with the same angular velocity as the original solution, σ is constant in time.
Let R(t) define the rotation matrix corresponding to the infinitesimal rotation vector σ(t), such that any vector in the basis of the perturbed body x p (henceforth "perturbed basis") can be described in the basis of the steady body (henceforth "steady basis") as x s = Rx p . To first order in σ, R can be defined using indicial notation as where δ ij is the Kronecker delta and ikj is the Levi-Civita symbol. Since the cross products given by the second term on the right-hand side will occur frequently, we introduce the linear operator [σ] × such that [σ] × v = σ × v for any vector v. In indicial notation, [σ] × ij = ikj σ k ; thus, we can also write where 1 is the identity matrix. Now, we compute the time evolution of the perturbation. Due to the initial perturbation, the magnetic field in the perturbed basis will be h p = R T h and consequently the magnetization will also change to m p . Note that, although for a permanent magnet m p = m since the magnetization is constant, this is not true for a soft magnet. The angular velocity of the perturbed swimmer in the perturbed basis Ω p will then be, Since the vector σ describes an infinitesimal rotation due to perturbation, to lowest order in σ, its rate of change with time can be found as the difference between the perturbed angular velocity and steady state angular velocity, with both referenced to the steady basis: Next, we compute m p to leading order in σ. Since the magnetic field vector in the perturbed basis is h p = R T h, the change in the magnetic field due to the perturbation is ∆h = − [σ] × h. Then, the new magnetization will be, where the ijth Cartesian component of the matrix ∂m ∂h is ∂m ∂h ij = ∂m i ∂h j and depends on the magnetization model.
We use Equations (18) and (22) on the right-hand side of Equation (21) to obtaiṅ where in the second equation only terms to linear order in σ are retained to obtain the linearized dynamics of σ, and terms have been rearranged by using the properties of cross products. The right-hand side of Equation (24) can be written asσ = µ 0 VQσ for some matrix Q, and if all the real parts of the eigenvalues of Q are negative, it means that the perturbations of the swimmer for that solution decay and the solution is stable. In indicial notation, For the case of a permanent magnet, the magnetization is constant, thus ∂m ∂h = 0, and the result for Q agrees with the stability criterion previously derived in [32].

Stability in Linear Magnetization Response Regime
To evaluate the stability using Equation (24), we need to compute ∂m ∂h . In the linear phase, we simply obtain ∂m ∂h = χ, thus Q ij = −M kl lnq m n h q ikj + M ik kln m l nqj h q − M ik kln h l χ nq qwj h w .

Stability in Saturated Regime
In the saturated phase, our model relates the spherical coordinates of the magnetization µ = (θ m , φ m , m s ) to the spherical coordinates of the applied field η = (θ h , φ h , h), where we have introduced the vector parameters µ and η such that, with the index i = 1, 2, 3, the components µ i and η i refer to θ m , φ m , m, or θ h , φ h , h, respectively. Thus, to evaluate ∂m ∂h ij and ∂η ∂h Additionally, the derivatives of the magnetization spherical coordinates with respect to the field spherical coordinates are with obtained from Equation (6). Combining the above, the change in the magnetization due to changes in magnetic field in Cartesian components can be written using the chain rule as which can be readily evaluated for any steady solution specified by (θ h , φ h , h) and (θ m , φ m , m s ). Equation (26) and the combination of Equations (25) and (32) constitute the criteria for stability that are the main results of this paper.

Example of Stability Evaluation
In this section, we provide examples for the evaluation of stability of steady solutions for a soft magnetic microrobot (see Figure 3) with dimensions in the same order of magnitude of currently fabricated examples. In the next section, we validate these results against numerical simulations. The head is a prolate ellipsoid with minor axis diameter D 1 = 2 µm and major axis diameter D 2 = 4 µm. The tail is a left-handed helix with four turns and arc-length L t = 36 µm, with helix diameter D t 2 = 2.8 µm, thickness D t 1 = 0.8 µm, and pitch P t = 36/4 = 9 µm. Since the ratio of major to minor axis of the ellipsoid is 2, the demagnetization factors n a and n r are 0.1736 and 0.4132, respectively. The volume of the head (magnetic material) is V = πD 2 D 2 1 /6 = 8.38 µm 3 and magnetic permeability of free space is µ 0 = 1.2566 × 10 −6 N A −2 . The entire swimmer is a rigid object with no deformation, and only the head is a soft-magnetic material; the tail is not magnetic.
We used two magnetic field strengths of h = 796 A m −1 and h = 3580 A m −1 equivalent to B = µ 0 h = 1 mT and B = µ 0 h = 4.5 mT. We also assume that the saturation magnetization of the magnetic material is 10 4 A m −1 , in the same range as the common soft magnetic material YIG [44]. For this value of saturation magnetization and this geometry, most but not all of the steady solutions are in the linear regime of magnetic response for the smaller magnetic field and in the saturated regime of magnetic response for the larger magnetic field. Figure 3b,c shows directions ofĥ relative to the swimmer head for these two field magnitudes. The directionĥ can lie anywhere on the light gray sphere around the ellipsoidal head; the tail of the swimmer (not shown) points downwards. The red and green lines are all directions ofĥ, which are steady corotating solutions; i.e., directions ofĥ which satisfy Equation (9). The blue region shows possible stable solutions; i.e., directions ofĥ for which the real parts of eigenvalues of the matrix Q in Equation (25) are negative. These stable directions are found by determining the magnetization vector corresponding with each direction and computing the eigenvalues of Q for that pair of h and m vectors. The (un)stable solutions are highlighted with (red) green in the figure. Note that, although parts of the red curve may seem to be lying inside the blue region, there is a very narrow space visible separating the two halves of each blue region and hence only the part highlighted by green is actually stable. The tail (not shown) is oriented downwards. The light gray sphere is a unit sphere; any directionĥ lies on the sphere. Steady solutions for which the microrobot and field co-rotate are associated with specific field directions plotted in red and green. The blue areas are the region of stability, for which the real parts of all the eigenvalues of Q in Equation (25) are negative when evaluated for that magnetic field direction and its associated magnetization. Only a portion of the steady solutions (green) are within the blue areas and hence also stable.
As discussed above, only the steady solutions that are also stable can be expected to be seen in experimental observations-these stable and steady solutions are the green curves inside the blue region. For any initial conditions, and an arbitrary magnetic field vector rotating in a plane, as long as the rotation rate of the field is slow enough that a stable steady solution exists, the swimmer reorients itself until it rotates steadily with the field in a way corresponding to a stable steady solution.

Validation of Analytical Stability Criterion against Numerical Time Evaluation of Trajectories
We validated our analytic stability criterion against stability as determined by numerical time evolution of the trajectory of the swimmer through time. Each steady solution is identified by its magnetic field direction in the frame of the swimmer, as in Figure 3. For each such magnetic field and steady solution, the solution is either stable, if after an infinitesimal perturbation the field evolves in time back to the original direction, or unstable, if after an infinitesimal perturbation the field evolves in time to some other direction.
To numerically propagate the trajectory of the swimmer, we assume that at time zero the body basis is the same as the lab basis and hence the direction of the field in the lab basis is the same as its direction in the body basis. At every timestep, we use the magnetic field to calculate the magnetization in the soft magnetic model, the torque (Equation (7) on the swimmer, and then the translational and angular velocity of the swimmer (Equation (8). We rotate the field in the lab frame by its prescribed angular velocity; if the field corresponds to a steady solution, its angular velocity is given by Equation (15) in the body basis. We rotate the swimmer at its calculated angular velocity. We use the quaternion method [45] to compute the transfer matrix between the lab and body bases at each timestep. In Figure 4, we show three examples of time-evolution, as recorded by trajectories of the magnetic field direction in time in the body frame. In Figure 4a, the initial field direction corresponds to a steady and stable solution, and under time evolution the field direction remains constant since the body and field corotate. In Figure 4b, the initial field is a steady but unstable soltuion, and under time it evolves to a different steady solution that is stable. In Figure 4c, the initial field does not correspond to a steady solution, and it evolves under time to a steady and stable solution.
To test stability, we start with an initial field direction which corresponds to a steady solution. To introduce an infinitesimal perturbation, we rotate the magnetic field in the body frame 3 s after the start of the simulation along a random direction with an angle of 1 • . To check the stability of that solution, we compare the starting position of the field in the body frame to the final position of the field in the body frame after the time evolution of the trajectory continues for some time (typically between 20 and 300 s). If the difference between the starting and ending position of the field |h 0 −h 1 | |h 0 | is less than 1%, we consider the solution to be stable.
We tested the predictions of our analytical stability criterion with this numerical time evolution method for the two magnetic field strengths used in Figure 3 and the results are shown in Figure 5. We picked some of the steady solutions based on a coarser grid and performed the numerical time evaluation for those solutions. We found that all the points that were predicted to be stable based on our analytic criteria were also shown to be stable through time evaluation (shown in green), and none of our analytically predicted stable points were shown to be unstable numerically. When we tested the steady solutions that were analytically predicted to be unstable, we found that all the solutions were also unstable (shown in blue). However, some of those points (shown in red) required increased simulation times (200-300 s) to be found numerically unstable; we believe this is due to them being close to marginally unstable, with slow growth in perturbations.

Effect of Head and Tail Geometries on Stability
We used our analytic criterion to investigate the effect of head and tail geometry on which types of solutions are stable. In general, steady solutions can be categorized as propulsive, if the rotation axis is along the helical axis, or tumbling, if the rotation axis is not along the helical axis [28,30,35]. Here, we describe the types of solutions obtained and whether they are stable or not for four different cases corresponding to combinations of two different head shapes (prolate and oblate ellipsoidal heads) with two different tail lengths. The prolate head is identical to that of the previous sections, with major to minor axis ratio of 2, and the oblate head has minor to major axis ratio of 2 but the same volume as the prolate head. The longer tail geometry is identical to that in the previous section, with four turns, while the shorter tail has the same pitch, radius, and diameter but only one turn.
The results are shown in Figure 6. In each panel, we categorize the steady solutions as either propulsive or tumbling according to whether the angle between the angular velocity and helical axis is less than or greater than 20 degrees. We further color-code the steady solutions as either stable or unstable. We find that, for shorter tails, tumbling solutions are seen over a greater rotational frequency range, consistent with previous investigations which show that swimmers with lower aspect ratio have more tumbling solutions [35]. Whether the ellipsoidal head is prolate or oblate has a large effect on stability; prolate ellipsoidal heads can have stable propulsive and tumbling solutions, while oblate ellipsoidal heads have very few stable solutions, only in the propulsive branch very close to zero rotation frequency.

Discussion
We have derived and validated an analytical criterion that can be used to test the stability of steady solutions for rigid microrobotic swimmers made of soft magnetic materials that are actuated by a rotating magnetic field to propel through bulk fluid. Although soft magnetic microrobots are commonly fabricated, such a criterion for the stability of their dynamics has not yet been presented. This analytic criterion is much more convenient than numerical methods used to test stability by directly time-evolving the dynamics of a steady solution after subjecting it to a small perturbation.
Our results show how the stability criterion we previously derived for microrobots with a permanent magnetization is altered by the response of the magnetization to the applied field. Additionally, superparamagnetic swimmers [29] have also been fabricated, and the result we obtained for the linear regime of soft magnets in our model is more broadly applicable to paramagnetic microrobots as well. We investigated the effects of tail length and head geometry on stability of swimmers, and our methods can be used to guide the rational design of microswimmer geometry. Funding: This work was funded by National Science Foundation awards CMMI-1650968, CBET-1252182 and CBET-1805847.

Conflicts of Interest:
The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.