## Abstract

Micro/nanomachines capable of propulsion through fluidic environments provide diverse opportunities in important biomedical applications. In this paper, we present a theoretical study on micromotors steered through liquid by an external rotating magnetic field. A purely geometric tight upper bound on the propulsion speed normalized with field frequency, known as propulsion efficiency, δ, for an arbitrarily shaped object is derived. Using this bound, we estimate the maximum propulsion efficiency of previously reported random magnetic aggregates. We introduce a complementary definition of the propulsion efficiency, δ*, that ranks propellers according to their maximal speed in body lengths per unit time and that appears to be preferable over the standard definition in a search for fastest machines. Using a bead-based hydrodynamic model combined with genetic algorithms, we determine that δ*-optimal propeller deviates strongly from the bioinspired slim helix and has a surprising chubby skew-symmetric shape. It is also shown that optimized propellers with preprogrammed shape are substantially more efficient than random magnetic aggregates. We anticipate that the results of the present study will provide guidance toward prospective experimental design of more efficient magnetic micro/nanomachines.

## INTRODUCTION

Development of artificial micro- and nanomachines that can either self-propel or be steered through fluidic environments in a controllable fashion is at the cutting edge of nanotechnology. Different approaches—ranging from catalytic nanowires to thermally, light-, and ultrasound-driven nanomachines—are actively being investigated [see (*1*) for review]. One of the promising techniques is engine-less and fuel-free navigation powered by magnetic actuation. Inspired by bacteria that spin their helical flagella to swim at low Reynolds number, researchers have recently proposed artificial magnetic micro- and nanohelices powered by uniform rotating magnetic fields, and various methods—such as “top-down” approach (*2*), glancing angle deposition (GLAD) (*3*), delamination of magnetic stripes (*4*), and direct laser writing (*5*)—have been developed for microfabrication of micrometer- and submicrometer-size (*6*) helical structures. The obvious advantage of using uniform rotating magnetic field for steering microscopic objects over the standard gradient methods is low field strength, typically of the order of a few millitesla.

In cases of both biological and artificial helices, the propulsion owing to rotation-translation coupling arises due to the symmetry-broken chiral shape. Some recent findings indicate that geometrically achiral planar structures (*7*) and randomly shaped aggregates (*8*, *9*) can be efficiently steered through liquid by a rotating magnetic field similarly to helices. One can argue that symmetries (e.g., chirality) of the driven object cannot be determined solely by its geometry and that symmetry properties of the magnetic (or electric) dipolar moment affixed to it must be considered (*10*, *11*). In particular, the dipolar moment affixed to a planar achiral shape can render it chiral, resulting in a unidirectional motion closely resembling propulsion of, for example, right-handed helix. However, for certain orientations of the dipolar moment, the structure remains achiral, whereas it can yet propel quite efficiently owing to a spontaneous symmetry breaking (*11*). These findings may have important practical implications, because planar micro- and nanostructures are easier to fabricate (using, e.g., standard lithography techniques) than three-dimensional helical motors, whereas randomly shaped aggregates in (*8*, *9*) form spontaneously and do not require sophisticated microfabrication techniques. Therefore, an interesting question is how do the planar and randomly shaped motors compare with propellers with a preprogrammed shape (e.g., bioinspired microhelices). In a broader context, how should one define the propulsion efficiency, and what is the optimal shape of an externally driven propeller?

Notice that the earlier optimization studies (*12*, *13*) relied on an a priori assumption of precession-free rotation of the magnetic object, which implies a prescribed magnetization transverse to the rotation axis. Thus, the search for the best structure involved only optimizing the chiral shape. The present approach does not make such an assumption and optimizes both the shape and the magnetization. The resultant optimal structures are quite different from those in (*12*, *13*); in particular, we demonstrate that achiral planar shapes are nearly optimal propellers.

## RESULTS

### Magnetically driven propulsion

Let us consider dynamics of an arbitrary shaped and magnetized rigid microscopic object immersed in an incompressible viscous liquid and actuated by uniform rotating (in *XY* plane) magnetic field , where *H* and ω are, correspondingly, the field amplitude and its angular frequency. We further assume that the permanent magnetic moment **m** is affixed to the object.

The motion of the magnetized body is force-free and driven solely by the magnetic torque **L**_{m} = **m** × **H**. In the Stokes (zero Reynolds number) approximation, the condition of the balance of forces and torques on the object by an incompressible Newtonian fluid reads(1)Here, **U** and **Ω** are the translational and angular velocities of the magnetic object, and and **F** are the coupling and rotation viscous mobility tensors, respectively. The triad of unit eigenvectors, {**e**_{1}, **e**_{2}, **e**_{3}} of **F** corresponding to the respective eigenvalues satisfying _{1} ≤ _{2} ≤ _{3}, defines the body-frame principal rotation axes. The lab-frame unit vectors are related to the body-frame axes {**e**_{1}, **e**_{2}, **e**_{3}} by a rotation matrix **R** parameterized by the Euler angles φ, θ, and ψ describing the instantaneous orientation of the object in the lab frame.

We are interested in the regime, whereas the propeller turns in sync with the actuating field (provided ω is not too high), rotating about the *Z* axis with angular velocity **Ω** = ω. This condition turns the second equation in Eq. 1 into a system of equations for the three angles ψ, θ, and . It has been demonstrated in (*10*) that the number of stable in-sync solutions corresponding to constant values of ψ, θ, and is at most two.

Knowing the dynamic orientation of the object and thus the magnetic torque, **L**_{m}, the translational velocity **U** can be readily found from the first equation in Eq. 1 as . By symmetry, the average velocity in an in-sync regime is along the *Z* axis. Taking a scalar product of both sides of this equation with **Ω** = ω, we readily obtain the projection of the propulsion in-sync velocity on the field rotation *Z* axis in a compact covariant form as(2)where **Ch** is a dimensionless chirality matrix given by the symmetric part of , with *ℓ* being the propeller size and being the normalized (unit) angular velocity. It is most convenient to write the right-hand side (r.h.s.) of Eq. 2 in the body frame in which **Ch** is fixed and determined solely by the geometry of the object and is expressed via the Euler angles aswhere we used the compact notation *c*_{ψ} ≡ cos ψ, *s*_{θ} ≡ sin θ, etc. Notice that both the diagonal (owing to the object’s chirality) and off-diagonal (which do not necessitate chirality) terms of **Ch** can contribute to net propulsion in Eq. 2 [see (*10*, *11*) for the role of symmetry in driven propulsion]. One may demonstrate that **Ch** (in contrast to ) is independent of the choice of coordinate origin, and it transforms as a (symmetric) pseudo-tensor under rotation of coordinate frame.

### Propulsion efficiency

The absolute value of the propulsion velocity scaled with ν*ℓ*, where ν = ω/2π is the cyclic frequency, is known as propulsion efficiency, δ ≡ |*U*_{Z}|/ν*ℓ* (*8*, *9*, *14*), and quantifies the displacement compared with the propeller’s size per one revolution of the field (*15*). The expression **Ω** ⋅ **Ch** ⋅ **Ω** /|**Ω**|^{2} in the r.h.s. of Eq. 2 is a Rayleigh quotient (*16*), and it is bounded by λ_{max}, the largest (by an absolute value) eigenvalue of **Ch**. Thus, for the propulsion efficiency, we have(3)where the equality holds when is aligned with the eigenvector corresponding to λ_{max}. Therefore, the maximum propulsion efficiency, δ_{max}, of an arbitrary magnetic propeller is a purely geometric property that can be determined without solving the accompanying problem of driven rotation. The bound on δ in (*3*) is tight, meaning that for an arbitrary in-sync actuation frequency ω, one can always find ω-dependent orientation of **m**, resulting in a rotation with .

Propulsion efficiency δ has a clear physical interpretation, and Rayleigh quotient provides an elegant geometric means of estimating its maximum value—magnetic propellers can be ranked by the maximal eigenvalue of **Ch**. However, this definition does not provide the complete information about the magnetic propeller. In particular, λ_{max} sets an upper bound on the slope of velocity-frequency curve without taking into account (i) potential nonlinearity and (ii) the applicable range of in-sync actuation frequency. As we will show below, maximizing δ at different frequencies produces rather different velocity-frequency relationships. Moreover, the propeller with a higher value of δ may have lower step-out frequency, ω_{s-o} (*17*), than the one with a lower value of δ. As a result, the latter propeller not only may be actuated in sync at higher frequency but also may exhibit higher speed at the step-out than the former seemingly more efficient propeller with a higher value of δ.

### Modified propulsion efficiency

Here, we put forward an alternative definition of the propulsion efficiency that ranks magnetic propellers by their speed at the step-out. Using Eq. 2, the dimensionless propulsion velocity at the step-out reads(4)where is the step-out frequency scaled with magneto-viscous frequency ω_{*} = *mH*_{*}, with _{*} being a representative rotational mobility and stands for the normalized angular velocity at the step-out. The modified propulsion efficiency is thus defined as δ* ≡ |*U*_{Z, s-o}|/ν_{*}*ℓ*, where ν_{*} = ω_{*}/2π. The propellers exhibiting a higher value of δ* have higher propulsion velocity in body lengths per unit time at the step-out. Nevertheless, weakly chiral (compact) propellers may have a high value of δ* at ω_{s-o} well beyond experimentally accessible frequency range. Thus, in general, the definitions of δ and δ* are complementary and should supplement each other to characterize the efficiency of the externally driven propeller. There is no simple geometrical constraint, similar to Rayleigh quotient, on δ*. Moreover, magnetization of random aggregates in (*8*, *9*) (and thus ) cannot be controlled, and thus, the use of δ* for their ranking is not practical. However, both definitions can be applied to optimize and rank microfabricated [e.g., using the GLAD method (*6*, *18*)] propellers with preprogrammed shape and magnetization.

### Random-shaped aggregates

We use particle-cluster (PC) and cluster-cluster (CC) aggregation algorithms (*19*–*21*) to generate multiple random fractal-like clusters made of *N* spherical beads with radii *a*. The geometry of the clusters satisfies the scaling relationship , where *D*_{f} and *k*_{f} are the aggregate fractal dimension and the fractal prefactor, respectively, and *R*_{g} is the aggregate’s gyration radius. In general, the exponent *D*_{f} controls the global morphology of the aggregate: *D*_{f} ~ 1 corresponds to filamentous aggregates, whereas the values *D*_{f} ~ 3 correspond to sphere-like object and the prefactor *k*_{f} controls local morphology. The representative geometry of the clusters made of 100 beads with different values of *D*_{f} is depicted in Fig. 1.

The viscous mobility tensors and and the chirality matrix **Ch** of the random aggregates are determined via multipole expansion method (*22*). The analysis of the driven rotation of an arbitrary shaped magnetic object is isomorphic to that of a triaxial ellipsoid (*10*) and requires knowledge of just two hydrodynamic quantities: the longitudinal *p* = _{3}/_{⊥} ≥ 1 and the transverse ε = (_{2} − _{1})/(_{2} + _{1}) ≥ 0 rotational anisotropy parameters. Here, _{1} ≤ _{2} ≤ _{3} are the eigenvalues of , and _{⊥} is the harmonic average of the two minor mobilities, _{1} and _{2}.

The computed values of *p* and ε are shown in Fig. 2 (A and B, respectively) in a wide range of *D*_{f} between 1.5 and 2.5 for two distinct values of *k*_{f} = 1.2 and 2.3. As expected, *p* is maximal for linear aggregates with the lowest value of *D*_{f} = 1.5. Transverse rotational anisotropy is quite small: on the average, meaning that, rotation-wise, all aggregates are cylinder-like to the first approximation with _{1} ≈ _{2} (*23*). The individual cluster exhibiting the highest value of ε ≈ 0.24 is shown in Fig. 2D (cluster #1). Notice that linear aggregates (*D*_{f} = 1.5 to 1.75) with larger prefactor *k*_{f} = 2.3 acquire lower values of *p* and higher values of ε when compared with aggregates with *k*_{f} = 1.2, because the former are locally thicker and thus less branched. Given the values *D*_{f} and *k*_{f}, the results are largely insensitive to the choice of aggregation algorithm (i.e., PC versus CC) and to the number of beads *N* composing the cluster (comparison not shown).

The major result is the maximal propulsion efficiency, δ_{max} plotted versus *D*_{f}, as shown in Fig. 2C. It can be seen that linear aggregates are better propellers than compact clusters, whereas on the average the best propellers show . The local structure of the linear aggregates has a moderate effect on the average value of δ_{max} (up to 20% difference), whereas compact aggregates with *k*_{f} = 2.3 perform significantly worse (by ~50%) than the corresponding clusters with *k*_{f} = 1.2.

The overall best individual (linear with *D*_{f} = 1.5) propellers showing δ_{max} ≈ 0.214 and 0.177 are arc-shaped (see clusters #2 and #3, respectively, in Fig. 2D). These findings are in agreement with experimental observations (*9*) where the “champion” random aggregate had an arc shape attaining δ ≈ 0.19. As we shall demonstrate below, maximizing δ for propellers possessing two mutually perpendicular planes of reflection symmetry results in arc shapes with similar values of δ_{max}. However, since not only the shape but also the magnetization of individual aggregates is random, the expected average propulsion efficiency δ will be much lower than not only the optimal value of δ_{max} ≈ 0.2 but also the mean value of . For example, in the experiments (*9*), the separation between the average value of propulsion efficiency and that of the fastest sample was about 10σ, with σ being the SD, whereas in our calculation this difference is smaller: ~3.6σ to 3.8σ. This discrepancy is due to the fact that random structures with uncontrolled (random) magnetization in (*9*) should perform much worse than , the ensemble average propulsion efficiency of randomly shaped but optimally magnetized clusters. The effect of magnetization on propulsion of a typical random aggregate is provided in the Supplementary Materials (see figs. S4 and S5).

### Optimization of magnetic micromachines

We further use genetic algorithms (GAs) combined with multipole expansion method [see (*10*) for details] to search for δ-optimal shapes of driven propellers under various symmetries. The genetic search was limited to propellers with filamentous morphology (i.e., using chains of beads as building blocks) because our results indicate that allowing for a more general morphology while being computationally expensive and time-consuming eventually converges to filamentous structures. In particular, we optimized δ for achiral shapes with (i) a single plane of reflection symmetry or (ii) two mutually perpendicular planes of symmetry, (iii) chiral skew-symmetric (*24*) shapes, and (iv) no symmetry (unconstrained optimization). The symmetry class governs the structure of the chirality matrix, for example, (i) zero-diagonal **Ch** matrix with two distinct off-diagonal entries, in (ii) a single off-diagonal entry, (iii) **Ch** with elements on the main diagonal plus a single off-diagonal term, etc. (*25*). The results of the optimization are gathered in Table 1 upon varying the aspect ratio of the filament the propeller is made of. Notice that, in all cases, increasing the slenderness of the filament improves the propulsion efficiency, and λ_{max} increases moderately with the slenderness.

For the twofold reflection symmetry, the corresponding eigenvalues of **Ch** are just , where (equal to either Ch_{13} or Ch_{23}, depending on the aspect ratio) is the only off-diagonal term (*10*). Maximizing yields nearly circular arcs. The optimal value of δ_{max} ≈ 0.19 found for a slender (1:12) arc using GA is close to that found for best fractal-like aggregates #2 and #3 in Fig. 2C. Notice that maximizing δ subject to a lower symmetry of a single plane of reflection yet converges to arc shapes with twofold reflection symmetry.

Maximizing δ subject to skew-symmetry results in chiral shapes resembling arc with twisted ends rather than a helix. The skew-symmetric propellers have higher values of δ_{max} than the arc-shaped motors (see Table 1). Moreover, unconstrained optimization converges to skew-symmetric shapes with practically the same values of δ_{max} as when skew symmetry was imposed explicitly (see Fig. 3 and table S1 for details). These findings indicate that skew-symmetric shapes are overall best δ-optimal propellers. For comparison, the reported values of δ for magnetic micro/nanohelices range from 0.03 (*2*) to 0.13 (*3*) [see the Supplementary Materials of (*9*) for a more detailed comparison of magnetic micro/nanohelices].

We further apply GA to search for δ*-optimal structures under the same symmetry constraints as before. It should be emphasized that the choice of in definition of magneto-viscous frequency ω_{*} in Eq. 4 is somewhat arbitrary. Here, we use the rotational mobility = (η*a*^{3})^{−1} based on the filament radius *a* and where η is a dynamic viscosity so that ω_{*} = *mH*/η*a*^{3}. Notice that the choice of affects the scaling of δ* with *a* and *ℓ* and thus the results of the optimization. The present choice based on the filament radius aims at comparison of propellers with a fixed thickness. However, one could have alternatively chosen to compare propellers of the same size *ℓ* assuming = (η*ℓ*^{3})^{−1}, or propellers with the same value of the rotational mobility _{⊥} and thus the same characteristic magneto-viscous frequency by taking ω_{*} = ω_{0} = *mH*_{⊥} (*10*). Whereas optimizing δ is purely geometric, optimizing δ* may also require the search for optimal orientation of the magnetic dipole moment **m**. The exact closed-form expressions for the optimal efficiency can be obtained under certain symmetries, for example, (ii) or (iii) (*10*), whereas for an arbitrary shape, the analytical expression for δ* is not known. Nevertheless, the transverse rotational anisotropy is typically small, ε ≪ 1, and in Eq. 4 can be closely approximated [with accuracy] by the cylindrical solution by setting ε = 0. The analytical solution for a driven magnetic cylinder is readily available [e.g., (*26*)], allowing derivation of an approximate expression for as a function of magnetization orientation and **Ch** (*27*).

The results of the optimization of δ* are shown in Table 2 upon varying slenderness . The representative shapes of the δ*-optimal propellers with twofold reflection symmetry and skew symmetry are shown in Fig. 4 (A and B) (see also figs. S2 and S3 and tables S2 and S3 for tabulated shapes and magnetization). The δ*-optimal shapes are similar to δ-optimal shapes, although the former are slightly more compact. This resemblance is reflected in the similarity of corresponding values of δ: for example, the structures in Fig. 4 attain near-optimal values of δ = 0.126 (A) and δ = 0.141 (B) (compare with δ_{max} = 0.129 and δ_{max} = 0.153, respectively, for δ-optimal propellers with the same slenderness in Table 1). The unconstrained search opted for shapes that are nearly skew-symmetric with similar values of δ* (see the last column of Table 2), although the discrepancy here is larger than that in Table 1 owing to cylindrical approximation used in unconstrained optimization of δ*. Some δ*-optimal symmetric shapes are depicted in Fig. 5 for several values of the aspect ratio. An important observation is that, although optimal propulsion efficiency, δ, increases with slenderness, maximizing δ* opts for plump structures (*28*). For structures with a twofold reflection symmetry, for instance, it can be shown that (*10*), where is the only nonzero dimensionless coupling coefficient. Because shows rather weak dependence on the slenderness (*a*/*ℓ*), it is always possible to maximize δ* by selecting a plump structure over a slim one. These findings explain why the three-bead cluster in (*7*) could be quite an efficient propeller. Notice that an alternative scaling of δ* with, for example, ω_{*} = ω_{0} = *mH*_{⊥} or ω_{*} = *mH*/η*ℓ*^{3} results in slightly different optimal shapes due to different penalization on the propeller’s size *ℓ*. Nevertheless, we found that the optimal speed *U*_{Z, s-o}/*ℓ* (in body lengths per unit time) diminishes with the increase in aspect ratio, confirming that chubby structures are faster propellers than the slim ones, regardless of the choice of representative magneto-viscous frequency ω_{*}.

To demonstrate the implication of the proposed concepts, in Fig. 6, we plot the scaled velocity *U*_{Z}/ω_{0}*ℓ* versus the frequency ω/ω_{0} for a particular δ*-optimal skew-symmetric propeller with an aspect ratio of 1:4 (see also the Supplementary Materials and table S4 for tabulated shape and magnetization). It can be seen that magnetization maximizing δ* yields the overall fastest propulsion that takes place at the step-out (thick red curve). Maximizing δ at different actuation frequencies 0 < ω < ω_{s-o} yields rather different binal velocity-frequency dependences (color solid and dashed curves, respectively). The maximum speed is attained when δ is optimized at the step-out when is aligned with the eigenvector corresponding to the maximal eigenvalue λ_{max} of **Ch**, whereas paired propulsive modes degenerate into a single solution (thick dashed blue curve). The thin (black) line with a constant slope equal to λ_{max} provides a tight upper bound to all δ-optimal velocity curves.

Purely transverse (to the easy rotation axis) magnetization typically applied to microhelices results in much slower speeds (dotted line), emphasizing the fact that optimal plumb propellers rely on large off-diagonal (pseudo-chiral) component of **Ch** that does not contribute to propulsion under purely transverse magnetization.

## CONCLUSIONS

Bioinspired helical micro/nanomachines that can be driven through fluids by a rotating magnetic field have been extensively studied during the past decade. This technology offers a remote, fuel-free, and engineless propulsion platform in a variety of fluidic environments, and it has been tested for various biomedical applications. For many such applications, the matter of speed, or efficiency, of the propellers is crucial. A number of recent studies raised the intriguing question of whether the bioinspired magnetic corkscrew is the optimal propeller, suggesting that propellers with achiral planar shapes [e.g., (*7*)] or randomly shaped magnetic aggregates (*19*, *20*) can compete with helical structures speed-wise. The present paper aims at answering this question.

Standard definition of the propulsion efficiency δ quantifies the displacement of the driven propeller per one revolution of the field in the regime, whereas the propeller rotates in sync with the field. Based on the compact covariant representation of the steady-state solution, we derived a tight bound on the propulsion efficiency in Eq. 3 given by the largest eigenvalue, λ_{max}, of the chirality matrix **Ch** that depends solely on the geometry of the driven object. Thus, one can now estimate the maximal value of δ for a prescribed geometry of the propeller without solving the corresponding dynamical problem. Using this bound, we numerically estimated the optimal performance of previously reported random magnetic propellers.

The conventional propulsion efficiency δ alone, however, does not fully characterize the propeller’s performance because it does not account for potential nonlinearity of velocity-frequency relationship and for the range of in-sync actuation frequency. Motivated by this observation, we put forward an alternative, or complementary, definition of the propulsion efficiency δ* based on maximal speed at the highest (or step-out) frequency of synchronous actuation. Because the ultimate aim is the search for the fastest propellers, one may argue that optimizing δ* should guide the experimental design, rather than δ. However, the limitation can be experimentally inaccessible high step-out frequency.

Using a bead-based hydrodynamic model combined with GAs, we numerically determined δ- and δ*-optimal shapes and magnetization. Apparently, optimal (by both definitions) propellers have skew-symmetric shape, resembling an arc with twisted ends that considerably deviate from a bioinspired helical geometry. Optimal planar structures are arc-shaped and fall short (up to ~30%) of the optimal skew-symmetric structures. Optimized slim structures exhibit higher values of δ than the corresponding optimal plumb structures, meaning that the latter exhibit a higher value of displacement per revolution. However, optimization of δ* shows that the chubby structures exhibit higher values of the propulsion velocity, in body lengths per unit of time, than the slim structures. This explains the efficient propulsion of the three-bead cluster reported in (*7*). Our results also reveal that, on average, random fractal-like aggregates are lousy propellers in comparison with motors with preprogrammed shape and magnetization. At the same time, the champion aggregates are arc-shaped with δ ≈ 0.2 similar to the δ-optimal planar structures in agreement with available experimental findings.

## MATERIALS AND METHODS

### Generation of random fractal-like aggregates

Random fractal-like aggregate of particles can be described by the following scaling law: , where *N* is the number of particles in the aggregate, *k*_{f} is the fractal prefactor, *D*_{f} is the fractal dimension, *a* is the radius of each particle, and *R*_{g} is the radius of gyration or the aggregate given by , with *R*_{i} being the distance between the aggregate center and center of the *i*th particle (*22*). *R*_{g} is the root mean square distance between the centers of *N* particles and the aggregate center. We used two different algorithms for generating aggregates that satisfy the scaling law for prescribed values of fractal dimension and fractal prefactor: PC and CC aggregation methods. The PC method creates a cluster by iteratively adding one spherical particle to a growing aggregate, whereas the CC method starts with an initial set of clusters, each consisting of very few particles, and iteratively combines the two smallest clusters. The PC method is simpler and allows for the generation of clusters that have a high fractal dimension (e.g., *D*_{f} = 2.5), which is not possible to accomplish using the more complex and time-consuming CC method (*29*). However, it was found that the CC method yields fractal-like structures with morphology similar to that of real (soot) aggregates (*29*). In this work, we used an implementation of the aggregation algorithms (*21*) based on the principles in (*19*).

### Multipole expansion and particle-based hydrodynamic model

We apply particle-based method (*22*) for computing mobility tensors and and the resulting chirality matrix **Ch** of random propellers and preprogrammed optimal propellers. This technique is based on multipole expansion of the Lamb’s spherical harmonic solution of the Stokes equations. The object is constructed from *N* (nearly) touching rigid spheres (see, e.g., Fig. 1), of radius *a*. 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 method yields a system of linear equations for the expansion coefficients, where the accuracy of calculations is controlled by the number of spherical harmonics (i.e., the truncation level), *L*, retained in the series. The validity and accuracy of the multipole expansion algorithm were previously tested in (*22*) against (i) the exact solution (in bispherical 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 to 30 spheres). For most of our calculations, the truncation level *L* = 3 yielded very accurate results. This particle-based approach was previously applied for modeling propulsion of, for example, a rotating helical filament (*18*, *30*), undulating filament (*31*), etc.

### Genetic search algorithm

To find optimal structures that maximize the target functions (see the next section), we applied a GA (*32*), which is a nondeterministic optimization method based on the principle of natural selection or “survival of the fittest.” This method can often find optimal or nearly optimal solutions for problems that deterministic approaches fail to solve in a reasonable time because of the exponential size of the search space. We considered structures made of spherical beads, for which the value of the target functions can be determined using the multipole expansion method. The GA involves the following steps: (i) generate a random initial population of filaments, each is composed of a fixed number of beads and meets the required symmetry constraints (if any); (ii) evaluate each propeller by the value of the corresponding target function; (iii) apply the mutation and crossover operators to the half of the population, which are the fittest individuals to create the next generation of propellers; and (iv) repeat steps (ii) and (iii) until convergence to an optimal solution. The structure of a random filament composed of spherical beads is determined by the sequence of vectors connecting the centers of every two adjacent beads. The mutation operator in step (iii) is implemented by randomly varying each of the vectors in this sequence with a fixed probability *p* (i.e., a vector remains unchanged with probability 1 − *p*). The crossover operator in step (iii) takes as input two filaments and combines them into a new filament by merging their two sequences of vectors. In the cases where the geometry and magnetization are coupled, each bead is also prescribed a magnetization vector that is randomly varied under the mutation operator. Constructing a structure that belongs to the desired symmetry class is done by applying the relevant symmetry operations. For example, taking the union of a nonplanar filament and its rotation by π about an axis passing through the first (last) bead results in a skew-symmetric shape (*24*).

### Target functions for δ* optimization

#### Twofold reflection symmetry

The permanent magnetic moment in the body-frame axes {**e**_{1}, **e**_{2}, **e**_{3}} is given by **m** = *m* (*s*_{Φ}*c*_{α}**e**_{1} + *s*_{Φ}*s*_{α}**e**_{2} + *c*_{Φ}**e**_{3}), where Φ and α are the spherical polar and azimuthal magnetization angles, respectively, and *s*_{Φ} ≡ sin Φ, *c*_{α} ≡ cos α, etc.

For a magnetic object made of beads of radius *a* with two mutually orthogonal symmetry planes, it can be demonstrated using the solution in (*10*) that the maximal value of the velocity is attained for polar magnetization angle , where *f* = _{3}/_{i} is defined as the ratio of the largest eigenvalue _{3} of the rotational mobility to that along the axis (*i* = 1 or *i* = 2) orthogonal to the body plane. The velocity then takes the form(5)Here, (no summation over the repeating indices) is the only nontrivial entry of **Ch**, where *i* = 1 or 2, depending on slenderness of the object. For , the optimal azimuthal magnetization angle α = π when *i* = 1 or α = −π/2 for *i* = 2. The optimal velocity in Eq. 5 can also be rewritten in terms of the (only) nonzero element of the coupling matrix as .

#### Skew symmetry

For the objects with skew symmetry [i.e., having a twofold skew-symmetry axis and no symmetry plane either normal to or containing this axis (*25*)], the closed-form expression for the maximal velocity at the step-out as a function of the magnetization polar angle Φ can be found (*10*):(6)where (no summation over the repeating indices) stand for the diagonal entries of **Ch** matrix, is the only nontrivial off-diagonal element of Ch, and λ = *f* tan Φ with *f* = _{3}/_{i}. Here, *i* = 1 or 2 depending on whether, respectively, **e**_{2} or **e**_{1} is aligned with the skew-symmetry axis. In deriving Eq. 6, we assumed that the shape is right-handed, that is, Ch_{3} > 0, and that orientation of the principal mobility axes is chosen in such a way that . Thus, the optimal azimuthal magnetization angle α = π and the optimal value of Φ is determined from maximizing the (positive) propulsion velocity in Eq. 6.

#### Cylindrical approximation for unconstrained optimization

Rotational mobility of arbitrary rigid object can be rather well approximated by that of a cylinder, as typically the transverse rotational anisotropy parameter ε ≪ 1 (*23*). Because driven rotation of a cylinder for ε = 0 has a closed-form solution [see, e.g., (*10*, *26*)], it can be used to derive an approximate analytical target function for optimizing δ* without any underlying symmetry assumptions. The normalized (unit) angular velocity can be written via Euler angles as . Within the cylindrical approximation, the step-out frequency is given by , where ω_{0} = *mH*_{⊥} stands for the characteristic frequency of the driven object and γ = *p* tan Φ, with *p* = *F*_{3}/*F*_{⊥} ≥ 1 being a longitudinal anisotropy parameter (*10*). Here, *F*_{⊥} stands for the harmonic average of the minor mobilities, . Substituting the limiting Euler angles sinθ|_{s-o} = (1 + γ^{2})^{−1/2} and ψ|_{s-o} = −π/2 − α at the step-out in the above expression for and then substituting it together with ω_{s-o} into Eq. 4 gives, after some algebra, the approximate target function for unconstrained optimization:(7)where, as before, (no summation over the repeating indices) stand for the diagonal entries of **Ch** matrix and denote the off-diagonal entries of **Ch**. The optimal magnetization angles Φ and α are found from maximizing the propulsion velocity in Eq. 7.

## SUPPLEMENTARY MATERIALS

robotics.sciencemag.org/cgi/content/full/3/17/eaas8713/DC1

Supplementary Text

Fig. S1. δ-Optimal propeller (Fig. 3).

Fig. S2. δ*-Optimal propeller (Fig. 4A).

Fig. S3. δ*-Optimal propeller (Fig. 4B).

Fig. S4. Typical random aggregate.

Fig. S5. Effect of magnetization on propulsion of a random aggregate.

Table S1. Tabulated shape of the δ-optimal propeller (Fig. 3).

Table S2. Tabulated shape of the δ*-optimal propeller (Fig. 4A).

Table S3. Tabulated shape of the δ*-optimal propeller (Fig. 4B).

Table S4. Tabulated shape of the δ*-optimal propeller (Fig. 6).

This is an article distributed under the terms of the Science Journals Default License.

## REFERENCES AND NOTES

**Acknowledgments:**

**Funding:**This work was supported, in part, by the German-Israeli Foundation for Scientific Research and Development (grant no. I-1255-303.10/2014), the Israel Science Foundation (grant no. 1744/17), and 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 (to K.I.M.).

**Author contributions:**Y.M. implemented the GA and performed all computations. O.D. assisted with implementation of the hydrodynamic algorithm. O.K. and K.I.M. assisted in developing the key ideas, research design, and implementation. A.M.L. directed the research and wrote the paper.

**Competing interests:**The authors declare that they have no competing interests.

**Data and materials availability:**All data needed to evaluate the conclusions in the paper are present in the paper or the Supplementary Materials.

- Copyright © 2018 The Authors, some rights reserved; exclusive licensee American Association for the Advancement of Science. No claim to original U.S. Government Works