### Self-assembly

We focus on self-assembled rafts constituted by *N* soft ferromagnetic beads of different diameters placed at a water–air interface and immersed in a vertical magnetic field **B**_{z} perpendicular to the interface. Figure 1 presents a few structures from regular to more complex assemblies. The total number of particles *N* ranges between 7 and 19, and the diameters *D*_{1}, *D*_{2}, and *D*_{3} used herein are 800, 500, and 400 μm. To simplify the notations, we will designate by *D* the diameter of the peripheral beads of an assembly in the following sections. We introduce the notation {*N*_{1}, *N*_{2}, *N*_{3}}, where the *N*_{i} are the number of large, intermediate, and small particles, respectively. For instance, Fig. 1a, b shows a {1,6,0} assembly. Figure 1c–f shows several other assemblies we will discuss in this paper. While the equilibrium distance between particles is independent of their size, larger particles are naturally located toward the center of self-assemblies due to gravity^{20}. The combination of long-range capillary pull and short-range magnetic repulsion typically leads to compact hexagonal structures, such as in Fig. 1b, d–f. However, the weight of the assembly can curve the water surface downward, which explains why fivefold symmetries are sometimes observed^{20}. With a larger particle in the center, sevenfold, and even eightfold symmetries can become stable owing to its stronger capillary pull, as shown in Fig. 1c.

### Mimicking ciliated locomotion

The main objective of this work is to develop a universal strategy to move self-assembled rafts composed of many particles, drawing inspiration from the metachronal rhythms used by ciliates and some arthropods. A common example of ciliated organism is the algae colony *Volvox*, which has a roughly spherical shape. Each individual cilium on *Volvox* follows a periodic beat with an effective and a recovery stroke, as shown in Fig. 2a. The blue cycle shows the position of the tip of the cilium over time. We exploit the fact that peripheral beads in magnetocapillary assemblies can describe a similar flow pattern (Fig. 2b). Ciliary beats are generally quite complex, asymmetric and can either be planar or nonplanar. However, the motion of a rigid sphere near a surface can produce a similar flow pattern, at least in the far field, and is therefore often used to model a cilium^{39,45,46,47,48}. In the case of an elliptical trajectory, one can show that the flow far from the sphere is proportional to the area of the ellipse projected on the plane perpendicular to the surface^{45}.

Remarkably, *Volvox* is able to change its swimming direction by changing the symmetry of the strokes on its body^{49}. When it swims in a straight line, metachronal waves propagate between two poles on its body, in the same direction as the effective strokes of the individual cilia^{47}. To change direction, for instance when exposed to light, *Volvox* can revert the strokes on part of its body^{49}. Figure 2c, d illustrates symmetries on *Volvox* that would lead to pure rotation and pure translation, respectively. The direction of the effective strokes is represented by the gray arrows, while the resulting motions are represented by the black arrow. This is a greatly simplified view, as the configurations adopted by *Volvox* lie on a continuum. Nonetheless, changing the symmetry of how the strokes are distributed on the body produces different behaviors. By using different symmetries for the variations of the magnetic field, we report in this paper how to change the direction of both the metachronal waves and the strokes, causing different behaviors. This is illustrated in Fig. 2e, f.

Focusing first on the {1,6,0} assembly from Fig. 1b, the central bead is surrounded by a shell of six neighboring beads forming a hexagon. These peripheral particles will play the role of cilia, as sketched in Fig. 2b. To generate metachronal locomotion, these simple cilia must become motile. The energy required to set the cilia in motion, which in eukaryotes is provided by the hydrolysis of adenosine triphosphate, will here come from the dipole–dipole interaction between the beads. Indeed, the addition of a horizontal field, for example, along the *x*-axis, can change the interaction energy. This causes a distortion in the assembly, as it rearranges to find a new minimum of energy. The interaction potential between neighboring beads is now given by

$${u}_{ij}=-{K}_{0}({x}_{ij})+frac{{rm{Mc}}(1+{beta }^{2})}{{x}_{ij}^{3}}left[1-3{cos }^{2}{theta }_{ij}right],$$

(2)

where *β* = *B*_{x}/*B*_{z}, and *θ*_{ij} is the angle formed by **r**_{ij} and the external field. When the horizontal component *B*_{x} vanishes, *β* goes to zero and *θ*_{ij} reaches *π*/2, such that Eq. (2) simplifies into Eq. (1). Additional information about this interaction can be found in Supplementary Note 1 and Supplementary Fig. 1.

In Fig. 3, quasi-static experiments have been performed by increasing the horizontal field *B*_{x} step by step, in order to observe the distortion of the structure. By tracking the distances *r*_{i} between the six beads and the central one (see Fig. 1a), one can build the dimensionless shape distortion vector

$${boldsymbol{sigma }}=frac{1}{6langle rrangle }mathop{sum }limits_{i=1}^{6}{{boldsymbol{r}}}_{i},$$

(3)

which is zero for a regular hexagon. Figure 3a shows *σ* = ∣** σ**∣ as a function of

*B*

_{x}/

*B*

_{z}. As expected, the distortion increases with the field ratio

*β*meaning that the shape deviates from a regular polygon. At a critical ratio

*β*

_{c}≈ 0.17, the system collapses, i.e., some beads come into contact. Figure 3b presents four snapshots of the assembly for

*β*= 0.0285,

*β*= 0.0857,

*β*= 0.1428, and the collapse at

*β*>

*β*

_{c}. An asymmetric deformation appears along the

*x*-axis, causing

*σ*to become nonzero as the particles in the direction of

**B**

_{x}(

**r**

_{i}and

**B**

_{x}parallel) move closer to the central one. This causes

**to orient in the direction opposite to**

*σ***B**

_{x}. A slight azimutal reorientation of the assembly is also observed. One could wonder why the interaction is not mirror-symmetric with respect to the

*y*-axis, as it is a function of the magnetic field squared. The origin of this asymmetry lies in the vertical positions of the particles. The center of the larger central particle sits deeper in the liquid, so that the angle

*θ*in Eq. (2) depends on the sense of the vector

**B**

_{x}. The continuous curve in Fig. 3a is a simplified model developed from Eq. (2), taking into account the fact that the central bead is larger than the others. More details on the model and the discussion on asymmetry can be found in Supplementary Note 2 and Supplementary Fig. 2.

### Rotational motion

Because of this asymmetry, the distortion of the assembly follows the direction of the magnetic field. By tracking it, we can follow the metachronal waves around the periphery of the assembly. As with ciliates, both the individual strokes and the phase shift between neighboring particles can produce a net flow. As will be shown later on, both effects would lead to a motion in the direction of the metachronal wave. We can therefore expect different behaviors as a function of the symmetry of the waves, as sketched in Fig. 2. In order to propagate the distortion around the assembly, we will first consider a counterclockwise rotating field along the horizontal plane, while keeping a constant field in the vertical direction. This can be written

$${B}_{x} = , {B}_{h}cos (2pi ft),\ {B}_{y} = , {B}_{h}sin (2pi ft),\ {B}_{z} = , {B}_{v},$$

(4)

with a typical rotation frequency *f* between 10^{−1} Hz and 2 Hz, leading to periodic distortions with period *T* = 1/*f*. An example trajectory is shown in Fig. 4a, as well as Supplementary Movie 1. Frequencies up to 10 Hz have been tested, as discussed in Supplementary Note 3 and Supplementary Fig. 3. One expects that dipole–dipole interactions along the horizontal plane will induce a distortion propagating on the periphery of the assembly^{29}. This distortion is shown in Fig. 4b, where a color scale from yellow to red is used to characterize the distance *r*_{i} between each peripheral bead *i* and the central one. The particle closest to the central one is colored yellow, and the farthest red. By following the yellow bead, one can see that the distortion is rotating along the periphery of the assembly, as denoted by the arrow. In a frame of reference rotating with the assembly, the trajectory of each peripheral bead describes a cycle as sketched in Fig. 2b.

Figure 4c shows the dimensionless distances (*r*_{i}(*t*) − 〈*r*〉)/*D*, where 〈*r*〉 is the distance in the absence of distortion. Each of the six peripheral beads is shown over two periods of the precessing magnetic field. To improve readability, the measurements for each bead were filtered at the frequency corresponding to the maximum of their Fourier transform (see “Methods” section). Each particle oscillates at the frequency of the magnetic field. The minimum of each *r*_{i}(*t*) over one period is denoted by vertical dot dashed lines, emphasizing that the oscillations of the beads are successively shifted by *T*/6. The small variations in amplitude are induced by the small differences in magnetic properties and wetting of the beads. This propagating distortion is the signature of a metachronal wave rotating around the body, in the direction opposite to the effective strokes as sketched in Fig. 2e. This is similar to the symmetry responsible for the rotational motion of *Volvox*, except in that case the wave is in the direction of the working stroke (Fig. 2c).

The resulting motion of the structure is a rotation at an angular speed <2*π**f*. This is illustrated in Fig. 4a, where the trajectory of the beads is shown over four periods *T*. This rotational motion comes from the cooperative motion of all beads. Each bead in the perimeter of the structure is roughly displaced by *D*/2 over a period *T*. Note that depending on the frequency and amplitude, different regimes can be observed under a rotating field, as explorer for three particles in ref. ^{29}. For instance, very low frequencies can lead to a locking with the external field, and for higher frequencies, inertia can start to play a role. A sweep in frequency can be found in Supplementary Fig. 3. Here, we limit our study to the regime where locomotion is due to nonreciprocal deformations.

### Translational motion

While using a precessing field to generate a rotation of the body seems natural, the symmetry needed to generate a translation is not trivial. In the typical ciliate depicted in Fig. 2d, metachronal waves appear and disappear at two poles on the body. A 2D equivalent would be two waves starting from the back and traveling on each side along the periphery of the assembly to the front of the system. The horizontal field should cause two waves of deformation propagating in the same direction, and vanishing at the poles. Starting from the field from Eq. (4), we can double the excitation frequency along one direction, or

$${B}_{x} = {B}_{h}cos (2pi ft),\ {B}_{y} = {B}_{h}sin (4pi ft),\ {B}_{z} = {B}_{v}.$$

(5)

This is the equation of a Lissajous figure with a ratio of 1:2, as sketched in Fig. 1a. Such a field describes a figure-8 pattern in the horizontal plane, or lemniscate. The resulting trajectory is shown in Fig. 5a, over 20 periods *T* = 1/*f* of the horizontal magnetic field. A similar trajectory is shown in Supplementary Movie 2 for the larger assembly {2,8,0}. For translational motion, the typical precession frequency *f* is taken between 10^{−1} Hz and 2 Hz. In Fig. 5b, two successive waves of deformation on each side of the assembly can be seen. The yellow circle, which denotes the peripheral particle closest to the central one, moves from back to front. The system is translating at a speed of ~10^{−1}*D*/*T*. This is comparable to the speeds that were obtained in the much simpler three-particle systems, ranging from 10^{−2}*D*/*T* in the collinear case to five 10^{−1}*D*/*T* in the triangular one^{26}.

Figure 5c displays the dimensionless distances (*r*_{i}(*t*) − 〈*r*〉)/*D* as a function of time, for two periods *T* = 1/*f*. The frequency spectrum of the distances *r*_{i}(*t*) shows peaks at *f* and 2*f*. To better visualize the movement of each bead, a filter at the frequency *f* has been applied (see “Methods” section). Two groups of three particles can be distinguished that correspond to the left and right sides of the assembly, their trajectories shown respectively in green and in blue. Each side is successively affected by the passage of the magnetic field. The deformation wave moves from back to front, with a shift of ~*T*/8 between the oscillations of the three particles. We use the same color code in Fig. 5a, c, where the beads from the back to the front of the assembly are colored light to dark, allowing to see the propagation of the wave.

The waves shown in Fig. 5 cause a translation toward the positive *y*-axis. To cause the assembly to swim toward the negative *y*-axis, the motion of the field must proceed in the opposite direction, for instance by inverting the sign of *B*_{y} in Eq. (5). It is possible to swim in any direction by changing the orientation of the lemniscate in the plane of the interface. For instance, inverting the *x* and *y* components in Eq. (5) causes a translation along the *x*-axis. By construction, any trajectory in the plane is therefore possible, following the same remote control strategy as in previously studied, low-particle-count cases^{24}.

Since the particles are mostly immersed in water, i.e., 80% of their volume is immersed, we will consider their Reynolds number ({rm{Re}}=rho vD/eta), with *ρ* and *η* the density and viscosity of water. For all motions discussed above, ({rm{Re}}) is between 10^{−2} and 10^{−1}. To eliminate any possible influence of inertia, additional experiments were performed on a water–glycerol mixture, increasing the viscosity by a factor 10 (10 ± 0.5 mPa). This allows to reach Reynolds numbers down to 10^{−3}. The behavior of the assembly is generally unchanged, as seen in Supplementary Movie 3. The only noticeable difference is in the swimming speed, which is ~0.8 times the speed on water. Such a small effect of a tenfold increase in viscous drag might seem surprising; however, this is not uncommon in microswimmers, as the swimming stroke can also be affected by the change in viscosity^{50}. Comparatively, few-particle magnetocapillary swimmers were considerably more affected by a change in viscosity^{15}.

### Locomotion mechanism

As with ciliates, both individual and collective effects between the peripheral particles could lead to locomotion. First, consider a single peripheral particle *i* and its position relative to the central sphere **r**_{i}, as defined in Fig. 1a, that we will express in terms of polar coordinates *r*_{i} and *α*_{i}. From Eq. (2), we see that under a horizontal external field **B**_{x}, the minimum of energy corresponds to the situation, where **r**_{i} and **B**_{x} are parallel. Because of the presence of other particles in the assembly, particle *i* can only move around its equilibrium position. Therefore, *α*_{i} will oscillate, following the direction of the horizontal field, when **r**_{i} and **B**_{x} are close to being aligned. Concerning distance *r*_{i}, as discussed in Fig. 3, it becomes smaller when **B**_{x} and **r**_{i} are parallel. The combination of these two effects is a trajectory that resembles the effective and recovery strokes of cilia, as shown in Fig. 6a. The exact shape of the cycle in (*r*, *α*) can vary, as the interactions are nonlinear and the neighboring particles also have an influence. However the basic characteristics of an effective stroke far from the central sphere (**B**_{x} and **r**_{i} antiparallel) and a recovery stroke closer to it (**B**_{x} and **r**_{i} parallel) are consistently seen. Because the recovery stroke always follows the direction of the wave, so does the induced motion. In that case, the direction of the wave is said to be antiplectic. By contrast, *Volvox* typically generates symplectic waves^{47}.

To visualize how a phase shift between neighboring particles can also lead to locomotion, Fig. 6b compares the angle *α*_{i} of a particle with the angles of its neighbors. One could similarly measure the distance between neighboring particles; however, the angle is slightly easier to interpret as it is decoupled from the variation of *r*_{i}. When three spheres on a line oscillate, they produce a net flow proportional to the sine of the phase difference Δ*ϕ* between the oscillations of each neighboring pair, as well as the product of their amplitudes^{51,52}. In the plane defined by their interdistances, this corresponds to the area of the enclosed cycle. The direction of the induced motion of the assembly is determined by the sign of (sin Delta phi) (ref. ^{51}), which in our case is imposed by the direction of the metachronal wave. As this phase difference, typically given by 2*π*/*N*_{t}, where *N*_{t} is the number of peripheral particles, is always smaller than *π*, the induced motion follows the direction of the wave.

In the case of an assembly, such as {1,6,0}, the peripheral particles are, of course, not aligned. The angle between two neighboring pairs is, on average, 60 degrees. The effect that this angle has on a collinear three-particle swimmer is to change the orientation of the velocity, causing a rotation whose radius depends on the deformation amplitude and the angle, with a local maximum of angular displacement close to 60 degrees^{53}. The effect of this angular displacement likely depends on the symmetry considered. In the translation case, the contributions from each side of the assembly of this angular displacement will cancel out over one period. However, the strokes left and right of the assembly are in phase opposition. This might explain the slight periodic rocking motion of the central particle that can be seen in Fig. 5a, as the rotation radius likely does not match the radius of the assembly. Concerning the rotation regime, both the tangential speed and this angular displacement would cause a rotation in the same direction^{53}. Whether this situation is favorable compared to a purely tangential contribution remains an open question.

One might wonder whether the individual stroke and the phase shift between neighboring particles contribute equally to locomotion, or if one effect dominates over the other. As discussed above, assuming a simplified geometry, the net fluid flow is in both situations proportional to the area enclosed by the cycles shown in Fig. 6. When expressed in units of angles and length ratios, the areas of these two types of cycles are comparable. However, the prefactor that links the nonreciprocal cycles with the fluid flow is typically given by the geometry of the system^{45,54}. To quantitatively compare the net quantity of fluid displaced in each case would require a comprehensive theoretical study taking fully into account the geometry of the system. Furthermore, in both situations, the motion always follows the direction of the metachronal wave. On the one hand, the recovery stroke corresponds to the moment where *r*_{i} is the smallest, which happens when **B**_{x} and **r**_{i} are parallel. On the other, the phase difference between neighboring particles, which determines the swimming direction in a three-particle swimmer, is also imposed by the course of the wave and would lead to a motion in the same direction. Therefore, while the motion observed likely is, as in ciliates, a result of both the individual strokes and the phase difference between neighbors, identifying their relative contribution would require a quantitative theoretical study, as the effects cannot be separated experimentally. However, looking at assemblies of different sizes might give additional insights into the mechanism, as we will see in the next section.

### Universality of the approach

While the assembly {1,6,0} has been used as an example throughout this paper, the wave-based approach has the advantage of remaining general. For instance, Supplementary Movie 2 shows the translational motion of the asymmetric raft {2,8,0} from Fig. 1f. In Fig. 7, we measure the speed *v* of seven different structures, which depends on the total number of particles *N*. Structures from *N* = 6 to *N* = 19 are shown, corresponding respectively to {1,5,0} and {1,6,12}. These structures all have a core of one or several larger particles, and one or two layers of smaller particles at their periphery. A much wider range of configurations has been successfully tested, including assemblies where *N* < 6, monodisperse assemblies, and asymmetric ones. However, to study the relation between *N* and *v*, we will only consider the most common and stable configurations. Particles in a magnetocapillary assembly tend to pack in a hexagonal or pentagonal symmetry^{19}. While a purely hexagonal lattice would be expected on a flat surface, the particles cause a local curvature of the interface. Other assemblies include ones with incomplete layers, and rarer configurations such as the octagonal {1,8,0} in Fig. 1c. These structures are often metastable, causing the particles to reorganize under variations of the applied field. As abrupt configuration changes and crowding effects can affect the average speed, only the more stable assemblies where the central particles are surrounded by five or six neighbors were considered in Fig. 7.

The average speed *v* is given in diameter (bar{D}) per period of oscillation *T*, where (bar{D}=sum {N}_{i}{D}_{i}/N) corresponds to the average diameter of the beads. The experimental parameters were identical for the various assemblies, with two exceptions. First, the vertical field *B*_{z} was adjusted to keep the distance between particles constant, as large assemblies tend to be more compact. Secondly, to minimize the influence of resonances that can appear in magnetocapillary interactions^{25}, every data point in Fig. 7 was averaged over several values of *f* between 0.5 and 2 Hz. When *v* is expressed in diameters per period, it is roughly independent of frequency *f* for the frequencies considered herein. Insets show typical assemblies and their trajectories. The red and purple dots distinguish assemblies where *N*_{3} = 0 and *N*_{3} ≠ 0, respectively. For instance, the two dots on the top right corner of Fig. 7 represent {1,0,5} for the purple one, and {1,5,0} for the red one. The reason for changing the size of the peripheral particles is to change the magnitude of some forces in the system, most notably the hydrodynamic and viscous forces. However, we can see that the points overlap for the same *N*, if we rescale the speed by the average diameter of the assembly. The speed *v* grows linearly with 1/*N*, meaning that the largest rafts tend to be less efficient under fixed conditions.

Though describing analytically, the dynamics of large assemblies would be complex due to the many degrees of freedom, it is possible to justify the 1/*N* scaling from Fig. 7. The thrust provided by the cooperative motion of peripheral particles can be written *F*_{t} ~ *N*_{t}*e*, where *N*_{t} < *N* is the number of these particles acting as motors (or cilia) and *e* is the efficiency of each motor. To simplify the discussion, we will consider that *N*_{t} is given by the number of beads in the outermost layer of the assembly. This is equivalent to considering that neighboring particles on successive layers act as a single cilium, as their motion is coupled. Moreover, we will make the assumption that, of the two mechanisms described in the previous section, the effect of the phase shift between neighboring particles is dominant. In this case, we expect the efficiency *e* of each motor to be proportional to the area of the nonreciprocal cycles as presented in Fig. 6b, and therefore to the sine of the phase shift between adjacent beads. For small phase differences, we find *e* ~ 1/*N*_{t}. In the particular case of precessing fields (Eq. (4)), the phase shift is simply given by 2*π*/*N*_{t}. As a result, the total thrust *F*_{t} is independent both of the total number of particles *N*, and of the number of motors *N*_{t}. At low Reynolds number, the total thrust *F*_{t} should be counterbalanced by the viscous forces *F*_{η} acting on all particles of the structure, i.e., *F*_{t} = *F*_{η} with *F*_{η} ~ *N**v*. As a result, the average translation speed of the assembly is given by

$$v(N)={v}_{0}/N .$$

(6)

This scaling is in agreement with the data from Fig. 7, and a fit provides the value *v*_{0} = 4.6 ± 0.2 in average diameter (bar{D}) per period *T*. Had we made the assumption that the stroke of a single peripheral particle was the dominant mechanism, efficiency *e* would have been independent on *N*_{t}. This would have led to the scaling (v=v_0^{prime} {N}_{t}/N). Indeed, even if crowding could affect the dynamics of a particle, the number of first neighbors is relatively constant for the motor particles in the structures considered here, i.e., with a complete outermost layer of particles on a hexagonal or pentagonal lattice. The results from Fig. 7 therefore suggest that the effect of the phase shift is dominant for the assemblies considered here.

The slope *v*_{0} is likely dependent on the experimental parameters, including particle-to-particle distance and amplitude *B*_{h}. Future studies could aim at maximizing *v*_{0} by changing the properties of the horizontal field. In Eq. (5), one could vary the relative amplitude as well as the phase difference between *B*_{x} and *B*_{y}, change the frequency ratio, or change the waveform altogether. For instance, a quadrifolium curve might generate two successive waves on each side of the assembly. One could also create a piecewise function such that the horizontal field never goes through zero. Indeed, when the field is interrupted, the assembly travels by a dimensionless coasting distance given by (d/D sim {rm{Re}} {rho }_{s}/rho approx 1{0}^{-1}), where *ρ*_{s} is the density of the particles^{12}. This means that any interruption in the waves is essentially lost, as coasting is very small at low Reynolds number.