Capillary assemblies in a rotating magnetic field

TL;DR

This paper investigates how rotating magnetic fields influence capillary-driven assemblies of metallic spheres, revealing various rotational behaviors and mechanisms that depend on magnetic properties, symmetry, and field parameters.

Contribution

It provides a detailed analysis of the effects of in-plane magnetic field rotation on capillary assemblies, highlighting different modes and underlying mechanisms involved.

Findings

Multiple rotational modes observed with different mechanisms.

Magnetic properties induce individual particle rotation.

Assembly rotation driven by non-reciprocal deformations.

Abstract

Small objects floating on a fluid have a tendency to aggregate due to capillary forces. This effect has been used, with the help of a magnetic induction field, to assemble submillimeter metallic spheres into a variety of structures, whose shape and size can be tuned. Under time-varying fields, these assemblies can propel themselves due to a breaking of time reversal symmetry in their adopted shapes. In this article, we study the influence of an in-plane rotation of the magnetic field on these structures. Various rotational modes have been observed with different underlying mechanisms. The magnetic properties of the particles cause them to rotate individually. Dipole-dipole interactions in the assembly can cause the whole structure to align with the field. Finally, non-reciprocal deformations can power the rotation of the assembly. Symmetry plays an important role in the dynamics, as well as the frequency and amplitude of the applied field. Understanding the interplay of these effects is essential, both to explain previous observations and to develop new functions for these assemblies.

Equations

$$\vec{m}=\frac{\chi V\vec{B}}{{\mu }_0}+{\vec{m}}_0,$$

$$\begin{array}{rl}{\vec{B}}_z& =B_z{\vec{e}}_z,\\ {\vec{B}}_h& =B_h\cos \left(\omega t\right){\vec{e}}_x+B_h\sin \left(\omega t\right){\vec{e}}_y,\end{array}$$

$$\frac{z}{D}\propto K_0\left(\frac{d}{l_c}\right),$$

$$U_c=-2\pi \gamma q^2{K}_0\left(\frac{d}{l_c}\right),$$

$$U_m=-\frac{{\mu }_0}{4\pi d^3}\left(3(\vec{m}\cdot \vec{e})(\vec{m}\cdot \vec{e})-\vec{m}\cdot \vec{m}\right).$$

$$U_m=\frac{{\mu }_0{m}_z^2}{4\pi d^3}+\frac{{\mu }_0{m}_h^2}{4\pi d^3}\left(1-3{\cos }^2\varphi \right),$$

$${\tau }_m=|\vec{m}\times \vec{B}|,$$

$${\tau }_{\eta }=8\pi \eta R^3\Omega .$$

$${\Omega }_{\max }\propto {\tau }_{m,\max }=m_{\mathrm{perm}}{B}_{h,0}\propto B_{h,0}.$$

$${\Omega }_{\max }\propto {\tau }_{m,\max }=m_{\mathrm{ind}}{B}_{h,0}\propto B_{h,0}^2.$$

$$U_{m,x}^{\mathrm{pair}}=\frac{{\mu }_0{m}^2}{4\pi d^3}\left(1-3{\cos }^2\varphi \right)=U_0\left(1-3{\cos }^2\varphi \right),$$

$$U_{m,x}^{\mathrm{iso}}=U_0\left(3\cos \alpha \cos 2\varphi -1+\frac{1-3{\cos }^2\varphi }{8{\sin }^3\frac{\alpha }{2}}\right).$$

$$\sin \theta =\frac{8\pi \eta R^3}{m{B}_h}\dot{\theta }=\tau \dot{\theta },$$

$$\theta =2{\cot }^{-1}\left(e^{-t/\tau +C_0}\right),$$

$$f_{\mathrm{cut}}=\frac{1}{\tau }=\frac{m{B}_h}{8\pi \eta R^3}.$$

$$\Omega \propto \sum _{i=1}^{N}2\pi f{A}_{d,i}{A}_{\theta ,i}\sin \left({\varphi }_i\right),$$

$$\Omega \propto \sum _{i=1}^{N}2\pi f{S}_i=2\pi fS,$$

$${\vec{F}}_{h,i}^t=\sum _{j\ne i}\left(\frac{9\pi \eta R^2}{2d}\left({\vec{v}}_j+\left(\vec{e}.{\vec{v}}_j\right)\vec{e}\right)+\frac{3\pi \eta R^4}{2d^3}\left({\vec{v}}_j-3\left(\vec{e}.{\vec{v}}_j\right)\vec{e}\right)\right),$$

$${\vec{F}}_{h,i}^r=\sum _{j\ne i}\frac{6\pi \eta R^4}{d^3}\left(\vec{d}\times {\vec{\Omega }}_j\right).$$

$$m_i{\vec{a}}_i=-\vec{\nabla }{U}_c-\vec{\nabla }{U}_m-6\pi \eta R{\vec{v}}_i+{\vec{F}}_{h,i}^t+{\vec{F}}_{h,i}^r.$$

$$I{\vec{\alpha }}_i=-{\vec{m}}_{0,i}\times \vec{B}-8\pi \eta R^3{\vec{\Omega }}_i,$$

$$\Omega =2{\Omega }_i{R}^3/d^3.$$

$$I{\alpha }_i=-m_{0,i}{B}_h\sin \left({\theta }_i-\omega t\right)-8\pi \eta R^3{\Omega }_i.$$

$${\ddot{\theta }}_i+\frac{8\pi \eta R^3}{I\omega }{\dot{\theta }}_i+\frac{m_{0,i}{B}_h}{I{\omega }^2}\sin \left({\theta }_i-u\right)=0.$$

$$\mathrm{Re}\sim \frac{{\rho }_f{R}^2\omega }{\eta }.$$

$$\xi =\frac{I\alpha }{8\pi \eta R^3\eta }\sim \frac{{\rho }_s{R}^2\omega }{\eta }.$$

Full text

\makeFNbottom

[TABLE]

††footnotetext: a* GRASP Lab, CESAM Research Unit, University of Liège, B-4000 Liège, Belgium; E-mail: ga.grosjean@uliege.be*††footnotetext: b* Helmholtz Institute Erlangen-Nürnberg for Renewable Energy (IEK-11), Forschungszentrum Jülich, Fürther Straße 248, 90429 Nürnberg, Germany.††footnotetext: c PULS Group, Institute for Theoretical Physics and Cluster of Excellence: Engineering of Advanced Materials, Friedrich Alexander University Erlangen-Nürnberg, Cauerstraße 3, 91058 Erlangen, Germany.††footnotetext: d Group for Computational Life Sciences, Division of Physical Chemistry, Institute Ruđer Boškovic, Zagreb, Croatia.††footnotetext: e Department of Applied Physics, Eindhoven University of Technology, P.O. box 513, NL-5600MB Eindhoven, The Netherlands.*

1 Introduction

As a way of simplifying the handling of microscopic objects, self-assembly has been intensely studied in recent years 1, 2, 3, 4, 5, 6. In particular, supramolecular self-assembly has allowed to undertake an impressive number of complex tasks including low-Reynolds-number mixing 7, drug delivery 8, assay of drugs 9, cargo transport 10, microfabrication 11, controlled self-propulsion 12, photonics 13 or sensoring 14. By requiring no direct intervention, self-assembly can offer clever solutions at length scales, ranging from the nanometer to the centimeter, that are difficult to access using traditional fabrication techniques, motors and actuators. One particular system uses magnetic spheres suspended at an interface to form reversible clusters called magnetocapillary self-assemblies 15, 16. These assemblies differ from other magnetic floating crystals 17, 18, 19, 20 in that they can form outside of a confining potential. In this case, the particles experience a pairwise lateral capillary attractive force 21, 22 as well as a magnetic dipole-dipole repulsive force. The resulting interaction potential resembles the Derjaguin-Landau-Verwey-Overbeek (DLVO) potential used in colloidal science and adsorption 23, with a primary minimum that corresponds to the generally irreversible contact between the particles, and a secondary minimum that corresponds to an equilibrium distance where the particles can be reversibly trapped 15.

The resulting structures can be made to perform several tasks, such as remote-controlled swimming 24, 12, transport and fluid mixing 25. They can also serve to study fundamental principles of low-Reynolds-number locomotion 26, 27, 28. Compared to other interface-bound systems 29, they do not use surface waves 30, 31 or surface tension gradients 32, 33, 34, 35 for locomotion. Instead, the role of the interface is to provide the attractive force and confinement necessary to form the assembly. The motion itself is powered by the non-reciprocal deformation sequence adopted by the assemblies when subjected to time-dependent magnetic fields.

However, several essential questions remain unanswered. First, the mechanism behind the rotational motion of the individual particles is not well understood, despite being necessary to explain the non-reciprocal deformation required for propulsion 12, 36. Secondly, the rotational dynamics of an assembly has never been experimentally studied, despite being used to change the swimming direction 12 or locally stir fluid without turbulence 25. For instance, Fig. 1 shows the stirring of fluids by an assembly under a rotating field. Dyes allow to visualize the flow, showing both a rotation of the individual particles and a slower rotation of the whole structure. In general, a trio of rotating particles can cause complex dynamics in the Stokes regime, despite the fact that the Stokes equations are linear and time-reversible 37. This can be an advantage for micromixing, particularly if there are multiple length and time scales involved. Finally, studying the influence of a rotating field on a magnetocapillary assembly can teach us more about the many swimming regimes observed 38 and the nontrivial influence of the oscillating field’s amplitude, the field’s frequency and the geometry of the assembly 12, 25.

2 Experimental setup and protocol

A glass container filled with water is placed at the center of a tri-axis Helmholtz coil system, as illustrated in Fig. 2. A transverse Hall Effect probe was used to obtain the ratio between the amplitude of the magnetic field and the current applied in each coil, as well as measure the spatial uniformity of the field. An integrated triaxis magnetoresistive sensor was also used to test the orthogonality between the coils. As the coils are order of magnitudes larger than the assemblies, about $50\text{\,}\mathrm{cm}$ in diameter, the nonuniformity of the magnetic field is negligible. Currents in the $x$ and $y$ coils are injected using a computer-controlled multichannel arbitrary function generator and a pair of amplifier, while the $z$ coils are fed by a regulated DC power supply.

Metallic spheres are deposited on the water surface with a pair of titanium tweezers. To minimize contamination, the tank is covered by a glass lid and the typical experimental runs are kept short, typically below an hour. Before each experimental run, the tank and lid are cleaned, and the water and particles replaced. The bottom side of the lid is coated with a transparent conductive oxide and connected to earth to prevent the accumulation of electric charges. A video camera fitted with a macro lens records the motion of the particles from the top. An LED array illuminates the bath from below. The particles are tracked using an algorithm based on circle Hough transforms.

The spheres used have a diameter $D$ of $500\text{\,}\mathrm{\SIUnitSymbolMicro m}$, are highly monodisperse, and are made of either steel or alloy steel, most commonly UNS G52986 and UNS S42000. Although they are 7.8 times denser than water, they float thanks to surface tension. Each particle is described by its polar coordinates $r_{i}$ and $\theta_{i}$ in the center-of-mass frame, as shown in Fig. 2. We will not consider displacements along the vertical direction, as they are typically very small. Indeed, for the range of particle-to-particle distances used herein, we can estimate using the Laplace equation that vertical displacements are typically around $1\text{\,}\mathrm{\char 37\relax}$ of $D$.

Under a magnetic induction field $\vec{B}$, the particles display a highly linear magnetization 39. This means that their magnetic moment is aligned with, and proportional to the field such that $\vec{m}\approx\chi V\vec{B}/\mu_{0}$ with $V$ the volume of the spheres and $\chi\approx 3$ their magnetic susceptibility. We consider here that the particles are perfectly isotropic and display no residual magnetism, or remanence. Indeed, the remanence measured in 39 is about two orders of magnitude smaller than $\chi VB/\mu_{0}$ for typical values of $B$. However, as can be seen in Fig. 1, the individual particles can rotate around their center. This can only be explained by taking into account some anisotropy in the properties of the beads, as the induced component is always parallel to the external field. Its origin will be discussed later on. We can create the following ansatz for the magnetization of the particles

[TABLE]

where $\vec{m}_{0}$ is a vector of fixed magnitude and orientation with respect to the bead that accounts for the anisotropy in magnetic properties. In general, though, the contribution of the remanence $\vec{m}_{0}$ is negligible.

The magnetic induction field used throughout this paper takes the form of a constant vertical field $\vec{B}_{z}$ and a horizontal field $\vec{B}_{h}$ of constant magnitude rotating parallel to the interface at an angular frequency $\omega=2\pi f$. The total magnetic field vector is therefore precessing around the $z$ axis at a frequency $f$. We have

[TABLE]

and the amplitude of the vertical field is kept constant at a value $B_{z}=$4.9\text{,}\mathrm{mT}$$. Changing $B_{z}$ would likely change the relative importance of the various effects discussed here, and therefore the frontiers of the regimes observed. However, to limit the number of control parameters in the experiments, we chose to mainly vary the excitation frequency $f$ and amplitude $B_{h}$.

3 Static assembly

In the absence of a magnetic field, the particles aggregate thanks to a lateral capillary force, a phenomenon known as the Cheerios effect 21, 40. Floating particles are surrounded by a meniscus, which in the case of a spherical particle of diameter $D$ is given by

[TABLE]

where $l_{c}=\sqrt{\gamma/\rho g}$ is the capillary length that relates surface tension $\gamma$ and unit weight $\rho g$, and $K_{0}$ is the zeroth order modified Bessel function of the second kind. This assumes small deformations compared to the capillary length.

In general, a solution of the Laplace equation cannot be found in the case of two particles. However, if the particles are separated by a large distance $d$, one can assume that the profile is given by the sum of both individual menisci. If a particle experiences a net vertical force, such as the resultant force of gravity and buoyancy, then an interaction potential arises from the product of such vertical force and the vertical displacement caused by the second particle. In other words, the slope of the interface can cause the particle to rise or fall along the meniscus, leading to a horizontal displacement. The potential energy associated with this interaction is given by

[TABLE]

where $q=L\sin\Psi$ is called the capillary charge, with $L$ the radius of the contact line and $\Psi$ the meniscus slope angle. This is assuming that the particles have the same capillary charge, which leads to an attractive force. The capillary charge is a measurement of a typical deformation length. We typically have $q\approx$10\text{,}\mathrm{\SIUnitSymbolMicro m}$$ for a $500\text{\,}\mathrm{\SIUnitSymbolMicro m}$ particle both from theory 22, 39 and from direct observations of the contact line. The validity of the linear and superposition approximations has been experimentally tested in 39. These approximations reproduce the experimental findings quantitatively granted that the particles are separated by a distance $d\gtrsim 2D$, which is always the case in the present paper.

The role of $B_{z}$ is to oppose this capillary attraction with a dipole-dipole magnetic repulsion. The pair potential for two particles of magnetic moment $\vec{m}$ is

[TABLE]

If we neglect the contribution of the permanent magnetization $\vec{m}_{0}$, this can be rewritten as

[TABLE]

where $m_{z}$ and $m_{h}$ represent the vertical and horizontal components of $\vec{m}$ and $\phi$ is the angle between the pair of particles and $\vec{B}_{h}$. This means that the contribution of the vertical field $\vec{B}_{z}$ is an isotropic repulsion.

The competition of the lateral capillary force and this dipole-dipole repulsion can lead to an equilibrium distance between particles. In the absence of a horizontal field, three particles assemble into an equilateral triangle 15. The addition of a horizontal field $\vec{B}_{h}$ can break the symmetry of the assembly, which deforms into an isosceles under a quasi-static increase of $B_{h}$ 12.

4 Dynamics of a triangular assembly

The angular frequency of the field $\omega=2\pi f$ is chosen to be positive, corresponding to a counterclockwise rotation. If we denote the average rotation speed of the assembly $\Omega$, we have $0\leq\Omega\leq\omega$. Figure 3 shows the overall motion of the assembly over 5 periods of the rotating field, as well as the trajectory of one of the three particles in more detail. These four trajectories are also visible in video format, as part of the online supplementary information accompanying this article. First, Figs. 3, 3 and 3 show the trajectories for three increasing values of the horizontal amplitude $B_{h}$ with the same applied frequency $f=$1\text{,}\mathrm{Hz}$$. It can be seen that an increase in $B_{h}$ leads to an increase in the rotation speed. The trajectory of each particle resembles an epicycloid, globally following a circle but punctuated by periodical cusps corresponding to each passing of the field. The instantaneous rotation speed can briefly dip below zero during these cusps. For larger values of $B_{h}$, the frequency of the cusps is doubled, suggesting a departure from the low-amplitude, linear regime that is discussed later on.

Figure 3 shows the trajectory of one particle for even higher values of $B_{h}$ and higher $f$. The rotation speed $\Omega$ is dramatically increased in this regime, performing a full turn in only about 6 periods of oscillation. Both the rotation speed and the amplitude of deformation dramatically increase. The trajectory is also quantitatively different, as there is no retrograde motion anymore. This regime is characterized by much larger deformations, somewhat resembling a juggling pattern. Typically, the particles oscillate sinusoidally, each separated by a phase of $2\pi/3$, as seen in the temporal evolution of $r$ and $\theta$ in Fig. 4. The radial and angular components are roughly in quadrature, so that the particles follow open cycles in the ($r$, $\theta$) plane. Note that for legibility, $\theta-\Omega t$ is used instead of $\theta$, which shows the evolution of $\theta$ in the referential of the rotating assembly. The cycles have been filtered using a Fourier transform to remove non-periodic parts of the signal.

Cycles in the configuration space are synonymous with non-reciprocal motion, a necessary condition for locomotion at low Reynolds number 41, 42. The Reynolds number in the experiment spans several orders of magnitude depending on $f$ and $B_{h}$. For instance, we have $\mathrm{Re}\sim 10^{-2}$ in Fig. 3 and up to $\mathrm{Re}\sim 1$ in the juggling regimes of Fig. 3. The fact that the individual particles in Fig. 4 each describe a cycle, with a phase difference between neighboring particles, might be reminiscent of the swimming mechanisms adopted by some living organisms. Ciliates are covered in slender organelles called cilia. They each individually describe a cyclic trajectory, with a phase difference between neighboring cilia such that a wave is propagated along the body of the bacteria 42. Recently, theoretical studies of ciliated motion have relied on particles describing cyclic trajectories 43, 44. Using these principles, magnetocapillary self-assemblies could serve as a simpler, self-assembled alternative to microfabricated systems 45, 46. With simple cilia composed of one particle each 47, 48 could be beneficial to improve our fundamental understanding of ciliated locomotion. Future work could expand this principle to larger, many-particle systems. However, fully characterizing the dynamics and mechanisms at play in the simpler, three-particle system is an essential first step.

Figures 5 and 6 show the influence of the two forcing parameters, respectively frequency $f$ and amplitude $B_{h}$, on the average rotation speed $\Omega$ of the system. The value of $\Omega$ is measured from a linear regression over typically about 100 revolutions of the magnetic field. Measurement errors are typically smaller than the symbols thanks to the high spatial and temporal resolution of the video acquisition. When error bars are included, they correspond to the deviation from the linear regression. This takes into account the oscillations of $\Omega$ over each passing of the field as well as slight variations of $\Omega$ that can occur over longer periods of time.

In Fig. 5, the frequency $f$ is varied for three values of the forcing amplitude. One can see that the rotation speed reaches a maximum that is dependent on the forcing amplitude $B_{h}$, as evidenced in the linear plot of $\Omega$. The bottom of Fig. 5 shows the same experimental data in a logarithmic plot and with rotation speed $\Omega$ nondimensionalized by the forcing frequency. This means that $\Omega/2\pi f=1$ corresponds to the assembly simply following the external field. This locking with the external field can only be observed for extremely low forcing frequencies, of the order of ${10}^{-2}\text{\,}\mathrm{Hz}$. For larger values of $f$, $\Omega/2\pi f$ follows a broken power law with a transition that corresponds to a maximum of $\Omega$. Dotted and dashed lines help to underline the power laws on each side of the maximum. In terms of $\Omega$, we have an increase $\Omega\sim f^{1/3}$, followed by the maximum, then a decrease $\Omega\sim f^{-3}$. The latter is discussed in Section 6. The position of the maximum depends on the forcing amplitude $B_{h}$ in a non-trivial way.

Figure 6 shows the influence of the forcing amplitude $B_{h}$ on rotation speed $\Omega$, for three values of the forcing frequency $f$. Three regimes have been distinguished. At low forcing amplitude (I), $\Omega$ increases with $B_{h}$. For the lowest values of the forcing frequency $f$, this increase is quadratic and independent on $f$, as shown by the logarithmic plot in Fig. 6. For intermediate amplitudes (II), $\Omega$ reaches a plateau. This intermediate regime is complex and highly nonlinear, with local extrema and a stronger dependence in $f$. At combined high amplitude and high frequency (III), $\Omega$ increases rapidly. This corresponds to the ‘juggling’ regime from Fig. 4.

This regime is discussed in more detail in Fig. 7, where the dimensionless rotational speed is measured as a function of $B_{h}$. At lower amplitudes, we have roughly $\Omega/2\pi f\approx 0.2$. Right before the particles come into contact, the speed drastically increases. This is accompanied by much larger deformations like those shown in Fig. 4. The theoretical value for the collapse of the assembly is also shown. One can see that, for excitation frequencies larger than $3\text{\,}\mathrm{Hz}$, the ‘juggling’ mode can persist above the theoretical collapse. Indeed, a static field $\vec{B}_{h,\;0}=B_{h,\;0}\;\vec{e}_{x}$ of more than $3.2\text{\,}\mathrm{mT}$ would cause the particles to come into contact 36. This indicates that the particles are not at equilibrium, and that the inertia of the particles therefore plays a role. Indeed, the Reynolds number in this case approaches unity. The bottom of Fig. 7 compares the rotation speed with the area $S$ of the cycles identified in Fig. 4. The relation between $S$ and $\Omega$ is discussed in Section 5.3.

5 Discussion

Three mechanisms can potentially play a role in the rotation of the assembly. First, if the particles are not perfectly isotropic, the rotation could be driven by the individual rotations of the particles. Secondly, the whole structure could reorient in an external field due to the pairwise magnetic dipolar interactions between particles. Thirdly, non-reciprocal deformations could power the rotation in a similar fashion to microswimmers.

5.1 Individual rotations

To unravel the origin of the ferromagnetic contribution to the total magnetization of the particles discussed in Eq. 1, we performed series of micromagnetic simulations which are based on the solution of the Landau-Lifshitz-Gilbert equation 49, 50 for the magnetization motion for relevant internal and external magnetic energy contributions in a bulk. We use the mumax3 simulation package optimized for graphics processing units (GPU) 51.

Though modern computational power allows for reaching dozens of millions of interacting magnetic cells, the physical size of simulated systems is limited to several (2-3) micrometers. At this scale we witness a rapid drop of the remanent magnetization (below one percent) of a single perfectly spherical bead mainly because of the formation of magnetic domains. A further increase in size up to desired $500\text{\,}\mathrm{\SIUnitSymbolMicro m}$ should result in even lower net magnetization. In addition, thermal fluctuations of the magnetic order parameter are known to enhance the effect of the magnetization reduction. Further effects like polycrystalline structure of the beads at such sizes should also lead to a vanishing remanent magnetization, as the average over all directions of polycrystalline axes in space should be zero. We conclude that the performed micromagnetic simulations suggest zero remanent magnetization for particles with diameters $500\text{\,}\mathrm{\SIUnitSymbolMicro m}$ and for parameters related to bulk steel alloy.

Non-sphericity of the beads might also explain the existence of the permanent internal magnetic moment. However, for particles used herein (precision ball bearings of AMBA grade 50), the deviation from a perfect sphere is at most $1.2\text{\,}\mathrm{\SIUnitSymbolMicro m}$ or $0.25\text{\,}\mathrm{\char 37\relax}$. Using known expressions of the shape anisotropy energy of an ellipsoid and assuming a monodomain particle 52, 53, we estimate that such a deviation in shape would only induce a magnetization of about $0.1\text{\,}\mathrm{\char 37\relax}$ of the saturation one. A more plausible explanation for the magnetic anisotropy is that the internal structure of the particles might be affected by the manufacturing process. Indeed, the spheres are formed by cutting the steel wire into small cylinders, which are then pressed into a spherical die and finally rounded. The process of pressing the cylinders into a die (heading) might create an anisotropy in the spheres, which could be at the origin the permanent component of the magnetization.

This permanent component is essential to account for the individual rotation of the particles. Evidence of such rotation can be found by tracing a spot on a particle with a marker pen. By correlation, we can then calculate the angle $\theta$ of the bead at any time. If we consider that we are in the Stokes regime, then the magnetic torque

[TABLE]

where $\vec{m}$ is the magnetic moment of the particle, has to be equal at all times to the viscous torque, which for a fully immersed sphere rotating at a constant speed $\Omega$ reads

[TABLE]

We can therefore expect a direct proportionality between $\tau_{m}$ and $\Omega$. In Fig. 8, a single sphere is exposed to a horizontal constant field $B_{h,\;0}$ that changes sign abruptly, leading the particle to rotate a half turn. The maximal rotation speed $\Omega_{\mathrm{max}}$ is measured for several values of $B_{h,\;0}$. Two possible evolutions of $\Omega_{\mathrm{max}}$ with $B_{h,\;0}$ can be expected. If the magnetic moment $\vec{m}$ responsible for the rotation is permanent and of magnitude $m_{\mathrm{perm}}$, we have

[TABLE]

On the other hand, if the moment responsible is an induced moment of magnitude $m_{\mathrm{ind}}=\chi VB_{h,\;0}/\mu_{0}$, then we have

[TABLE]

As can be seen in Fig. 8, the evolution of $\Omega_{\mathrm{max}}$ with regards to $B_{h,\;0}$ is linear, suggesting that the rotation of the individual particles is indeed an effect of a small remanent magnetization.

The top left inset of Fig. 8 shows the angle $\theta$ as a function of time during the inversion. One can see that a stronger field leads to a faster rotation. For higher values of $B_{h,\;0}$, there is a slight overshoot due to inertia. This is not surprising considering that, in order to facilitate the measurement, a larger particle of $1\text{\,}\mathrm{mm}$ in diameter was used. It was observed that the axis of rotation is always perpendicular to the interface. This might be surprising considering that the direction of $\vec{m}_{0}$ depends on the random orientation of the particle when dropped. However, even the strongest horizontal fields used herein, about $50\text{\,}\mathrm{mT}$, fail to produce any measurable out-of-plane rotation. These types of rotations might be prevented due to the additional energy cost required to move the contact line 28. The bottom right inset of Fig. 8 shows a side view of a particle. The contact angle is approximately $80\text{\,}\mathrm{\SIUnitSymbolDegree}$. The contact line is generally axisymmetric, though its position on the particle can vary depending on how the particle was placed at the interface.

These individual rotations likely power the rotational motion at the lower excitation amplitudes, where the effect of the permanent component of the magnetization dominates. This corresponds to regime (I) in Fig. 6. As the amplitude increases, the effect of the induced component of the magnetization becomes important through dipole-dipole interactions.

5.2 Dipole-dipole alignment

Let us consider only induced magnetic dipoles, which are always aligned with the external field. From Eq. 6, the potential energy associated with the horizontal field for a pair of particles can be written as

[TABLE]

where $U_{0}$ corresponds to the maximum energy for a given $d$, reached when the pair is oriented perpendicular to the field. The minimum of energy is reached when the pair is aligned with the field and is equal to $-2\,U_{0}$. This means that under a horizontal field $\vec{B}_{h,\;0}$ a pair of particles experiences a torque and tends to align with the field.

Let us now consider the case of three particles. For simplicity and because this situation is prevalent in the experiment 12, 25, we suppose that the particles lay on the vertices of an isosceles triangle. Let $\alpha$ denote the vertex angle of the isosceles and $\phi$ the angle between its base and the external field. Summing the pairwise energy contributions leads to

[TABLE]

For two values of $\alpha$ between 0 and $\pi$, this expression is independent on $\phi$, meaning that all the orientations of the assembly with regard to the field are energetically equivalent. These values are (i) $\alpha=\pi/3$, corresponding to an equilateral triangle, and (ii) $\alpha\approx 1.2605$, which is the value of $\alpha$ such that $\cos\alpha\sin^{3}(\alpha/2)=1/16$. In fact, by symmetry, the independence of the total magnetic interaction energy on the orientation of a rotating field $\vec{B}_{h}$ is verified for any regular polygon. This means that, for low values of $B_{h}$, the dynamics is driven by the individual rotations of the particles for which the torque is proportional to $B_{h}$. Conversely, we can expect the dipole-dipole interaction, quadratic in $B_{h}$, to play a role at high amplitudes, for which the deformation of the triangle is significant.

Considering that both individual rotations and deformations of the structure can be observed in regime (II), the complex dynamics observed in this intermediate regime might be due to a combination of the two effects. Both the induced and the permanent components of the magnetization expressed in Eq. 1 compete, though it is difficult to isolate the influence of each effect experimentally. Other effects might add to the complexity of the intermediate regime. For instance, self-assemblies of rotating particles have been shown to form due to hydrodynamic repulsions 19, though the typical magnitude of this force $\mathcal{F}\sim\rho\omega^{2}a^{7}/d^{3}$ is orders of magnitude smaller in our case. At even higher amplitudes, we can expect the contribution of $\vec{m}_{0}$ to become negligible. As we will show next, the dipole-dipole alignment mechanism also becomes ineffective at higher frequencies, such as in regime (III). A third mechanism is therefore needed to explain this high-deformation regime.

5.3 Rotational locomotion

Under a vertical magnetic field, magnetocapillary self-assemblies often adopt the shape of a regular polygon 24, with notable exceptions at $N=2$ and $N=5$. These regular assemblies only have a preferred orientation when the horizontal field is strong enough to significantly alter their shape. However, both this effect and the effect of the individual rotations of the particles should become inefficient at high frequencies, as the dipoles struggle to follow the external field. For instance, equating the magnetic torque from Eq. 7 to the viscous torque from Eq. 8 and neglecting inertia, we find

[TABLE]

which has a solution of the type

[TABLE]

where $C_{0}$ is an integration constant. This means that there is a typical reaction time $\tau$ needed for a dipole to align with an external field. In fact, looking at Eq. 14, we can see that it takes about $5\tau$ to cover $95\text{\,}\mathrm{\char 37\relax}$ of the half-turn. We can therefore define a typical cut-off frequency

[TABLE]

which is equal to $10\text{\,}\mathrm{Hz}$ for a field of about $0.6\text{\,}\mathrm{mT}$. This would mean that there is an observable decrease in the oscillation amplitude as early as $2\text{\,}\mathrm{Hz}$. Dipole-dipole interactions are usually weaker than dipole-field interactions in magnetocapillary assemblies, as the field generated from a particle is given by $\mu_{0}m/4\pi d^{3}$, which for an induced dipole is smaller by a factor $R^{3}/d^{3}\approx 1/64$ for particles at a typical distance $2D$. This means that we can expect an even more drastic cut-off in the case of the dipole-dipole alignment. This is consistent with the extremely low excitation frequency needed to observe a rotation of the whole assembly locked with the external field, as shown in Fig. 5.

This means that both previously discussed mechanisms cannot explain the sudden increase in rotational speed at combined high frequency and high amplitude described in Fig. 4. As deformations are high in this regime, one might wonder if the motion is powered by non-reciprocal deformations. In other words, the cycles shown in Fig. 4 could power the rotational motion, similarly to the translational motion of the Najafi-Golestanian swimmer 54. Let us assume that each particle, represented by its polar coordinates $d_{i}/D$ and $\theta_{i}$, oscillates harmonically both radially and angularly at a frequency $f$. We denote the corresponding amplitudes $A_{d,i}$ and $A_{\theta,i}$ and the phase between both signals $\phi_{i}$. By analogy with 54, we construct the ansatz of an expression for the rotational speed

[TABLE]

where we consider each particle to act as an independent motor whose contributions can be summed. As the trajectories in the $(d,\theta)$ plane are not necessarily ellipses, we can assume a more general expression

[TABLE]

where $S_{i}$ is the dimensionless surface enclosed by the trajectory of particle $i$. The expressions from Eq. 16 and Eq. 17 are equivalent in the case of an ellipse. Figure 7 compares the rotational speed and the dimensionless surface of the cycles for the data of Fig. 7. Most of the points are concentrated around $\Omega/2\pi f\approx 0.04$ and $S\approx 0$, corresponding to the plateau in Fig. 7, where deformations are small and the rotation is explained by the two previously discussed mechanisms. However, for fast rotations $\Omega/2\pi f>0.04$ we see a linear relationship between the combined dimensionless area $S$ and the dimensionless rotational speed $\Omega/2\pi f$, confirming our hypothesis.

6 Model

In order to address the question of the importance of each dynamics described in Section 5, we develop a numerical model of the experiments presented in this article. The algorithm integrates Newton’s equations for each bead. As the beads both translate and rotate on the surface, an inventory of all the forces and torques is required.

To describe the translation, one has to take into account the Cheerios effect introduced in Eq. 4, the magnetic dipole-dipole interaction introduced in Eq. 6 as well as the hydrodynamic interactions. The translation of a bead generates a stokeslet, i.e. the flow field originating from a point force 55. This flow drags neighboring beads with the following force

[TABLE]

where the labels $i$ and $j$ refer to two beads on the surface and $\vec{v}$ denotes their speed. To express the forces as a function of the translation and rotation speeds, the substitutions $v_{i}\rightarrow f_{i}/6\pi\eta R$ and $f_{j}\rightarrow 6\pi\eta Rv_{j}$, which are typically valid up to $\mathcal{O}(d^{-4})$, have been applied 56. Since the beads rotate in the plane of the interface, the flow field generated by a point torque, a rotlet, must also be taken into account 55. Using the same substitution 56, we find

[TABLE]

For simplicity, we considered fully immersed spheres. Previous investigations showed that the effect of the partial immersion is essentially to reduce the Stokes drag by a prefactor which was determined to be about $0.86$ for the particles used herein 39. A theoretical prediction for this factor can also be obtained from 57, yielding a prefactor of $0.55$ for a contact angle of $80\text{\,}\mathrm{\SIUnitSymbolDegree}$. This model assumes a perfectly flat interface and does not take the meniscus into account, which might explain the difference.

The forces associated with the Cheerios effect and the dipole-dipole interaction are obtained by taking the derivative of the corresponding potentials, i.e. Eqs. 4 and 6. Therefore the equations of motion are

[TABLE]

In this model, only pair interactions have been considered therefore justifying the expansion of the hydrodynamic coupling to the third order in $1/d$.

As stated previously the rotation of the beads comes from the magnetic torque arising from the residual magnetization. No torque can arise from the Cheerios effect under the assumption that the beads create a symmetric wetting line when dropped on the liquid interface. To describe the magnetic torque, we need to take into account the anisotropy in magnetic properties. Using the ansatz of Eq. 1, the rotation of the beads is described via the equation

[TABLE]

where $I$ is the moment of inertia and $\vec{\alpha}$ the angular acceleration. This equation is projected along the $z$ direction as the rotation of each bead lies in the plane defined by the surface. The equations of motion presented above are integrated numerically. Note that each parameter appearing in those equations can be measured experimentally, leaving no free fitting parameter. In order to illustrate the results provided by the numerics, Fig. 9 gives the dynamics of the swimmer for the same parameters as in Fig. 3. The dynamics obtained in simulations is essentially similar to the what is measured experimentally. For small amplitude and frequency, the deformation of the triangle is small and the rotation slow. When the amplitude is increased, the rotation becomes faster while the deformation of the structure becomes more important. Finally, the nonlinear juggling regime for high amplitude and frequency is also observed.

Some differences with the experiment can also be noted. First, a more pronounced retrograde motion is observed, as seen in Figs. 9 and 9 when compared with Figs. 3 and 3. Recovering these finer structures accurately might require a more precise description of the various terms in Eq. 20, including higher-order hydrodynamic terms and taking into account the partial immersion of the spheres. Secondly, the rotation speed is typically much lower than in the experiment, a notable exception being the high-deformation regime. Taking into account the inertia of the flow might increase the numerical prediction for $\Omega$. This explanation is insufficient, though, considering that the model is the most accurate when the Reynolds number is the highest. The solution might reside in the description of the magnetic properties of the particles, as it becomes more complex at low amplitude. The anisotropy of the particles’ properties comes into play, and a perfectly soft-ferromagnetic description becomes insufficient to describe what is observed. While the ansatz of Eq. 1 allows to recover the general behavior of the particles, the model would certainly benefit from a finer, three-dimensional study of their magnetic properties. A better understanding of the magnetic properties of the particles could also improve reproducibility. Indeed, the orientation of the anisotropy, modelled by $\vec{m}_{0}$, is not controlled. This can lead to variations in the response of a particle to the field, depending on how it was placed at the interface. Combined with the variations in contact line radius, this effect is likely the main cause of variations between experimental runs.

The numerical model allows to turn off some effects to determine their influence. As expected, removing $\vec{m}_{0}$ leads to no detectable rotation in the low-amplitude regime of Figs. 9 and 9. In the intermediate regime of Fig. 9, the overall speed is reduced but the motion qualitatively similar, showing that both the effect of $\vec{m}_{0}$ and the dipole-dipole interaction matter. Finally, the high-amplitude, high-frequency regime is almost not affected by the removal of $\vec{m}_{0}$, which confirms that the anisotropy of the particles is not at play there.

This model also gives the possibility to describe the power laws observed in Figs. 5 and 6. Let us assume that the triangle remains equilateral once exposed to the rotating magnetic field with one bead at each vertex. This hypothesis corresponds to the observation in the linear regime. As the beads rotate due to the external magnetic field, each of them is the source of a rotlet at the surface of the fluid. Since the Stokes equation is linear, the global flow field is a rotlet centered at the geometrical center of the assembly. This feature can be observed in Fig. 1 by following the dyes around each bead. The triangle therefore rotates with this induced flow. Interestingly, again thanks to the linearity of the Stokes equation, the speed of rotation of the triangle is proportional to the speed of rotation of each bead. Considering that each particle rotates at the same rate $\Omega_{i}$ and using the definition of a rotlet 55, one can find that

[TABLE]

As a consequence, quantifying the rotation of one bead allows to understand the global dynamics and to extract scaling which can be compared to experiments. Equation 21 can be written as

[TABLE]

Introducing the dimensionless time $u=\omega t$ and the corresponding time derivative with dotted symbols, one has

[TABLE]

This equation gives two dimensionless parameters $c_{1}=8\pi\eta R^{3}/I\omega$ and $c_{2}=m_{0}B_{h}/I\omega^{2}$, each depending on the frequency of oscillation of the magnetic field. The dashed line in Fig. 5 gives the average rotation speed $\Omega_{i}$ of a single bead as a function of the frequency $f$, using typical experimental parameters. Two regimes are recovered. Below $1\text{\,}\mathrm{Hz}$, the particle is locked with the external field so that $\Omega_{i}/2\pi f=1$. At higher frequency, we have roughly $\Omega_{i}\propto f^{-3}$, as the particle struggles to follow the field. Compared to the three-particle assembly, the rotation is generally faster and the uncoupling with the rotating field happens at a much higher frequency. The intermediate regime is absent, confirming that it is likely influenced by the dipole-dipole interactions in the assembly. In Fig. 6, $\Omega_{i}$ is shown as a function of $B_{h}$, which is proportional to $c_{2}$ as $c_{1}$ is kept constant. The curve has been adjusted on the experimental points, showing a very similar quadratic dependence in $B_{h}$. This means that in the low-amplitude linear regime, the assembly as a whole essentially behaves like a single effective dipole.

In general, it is not trivial to link the rotation of each particle $\Omega_{i}$ to the rotation of the assembly $\Omega$, as including inertia is necessary to reproduce the experimental observations in the higher-frequency regimes. Indeed, if we neglect the left-hand side of Eq. 23 and combine it with Eq. 22, we find $\Omega\propto f$ at low frequency and $\Omega\propto f^{-1}$ at high frequency. By taking the inertia of the particles into account, we find a third regime where $\Omega\propto f^{-3}$, which corresponds to what is observed experimentally at high frequency in Fig. 5. The single-particle study from Fig. 8 also shows a visible overshoot due to inertia, although this study was performed on a larger particle of $1\text{\,}\mathrm{mm}$ in diameter.

Nonetheless, inertia in Eq. 23 must be distinguished from the inertia in the flow. If we compare the relative importance of inertia over viscous damping in the flow, we find the Reynolds number defined as $Re=\rho_{f}UL/\eta$. Using $L\sim R$ and $U\sim R\omega$, one has

[TABLE]

Conversely, in Eq. 23, we can compare the inertial term $I\alpha$ with the damping one $8\pi\eta R^{3}\omega$. Using $\alpha\sim\omega^{2}$, one finds the dimensionless number

[TABLE]

This means that depending on the density of the fluid $\rho_{f}$ and the density of the particle $\rho_{s}$, it is possible to remain in a viscosity-dominated flow while having a non-negligible influence of inertia on the bead’s rotation. In the experiment, we have $\rho_{s}/\rho_{f}\approx 7.8$. Considering that $\mathrm{Re}$ can approach $1$ from below at the highest amplitudes and frequencies used herein, this density ratio means that we often have $\mathrm{Re}<1$ but $\xi>1$. Compared to the Stokes regime, this intermediate region where both inertia and viscosity matter is seldom discussed in the literature. Further studies on this aspect using magnetocapillary assemblies could therefore be beneficial to our general understanding of small-scale flows.

Particles of a few tens of microns in diameter could in theory still be bound by the magnetocapillary interaction 39, assuming that the magnetic properties of the particles scale properly. Indeed, both the magnetic and the capillary force scale like $\sim D^{6}$. Below this limit, thermal agitation start playing a significant role and affecting the translational and rotational degrees of freedom of the particles. The typical particle size below which thermal agitation could free a particle from a magnetocapillary bond was estimated in 39 to be around $3\text{\,}\mathrm{\SIUnitSymbolMicro m}$. In addition, the magnetic order will be greatly influenced by the thermal noise and the question of the presence of the permanent magnetic moment is hard to answer. How the rotational regimes would be affected by a downscaling thus remains an open question.

7 Conclusion

To conclude, rotating fields cause magnetocapillary assemblies to rotate thanks to three distinct mechanisms. First, individual particles rotate, similarly to a compass needle in a magnetic induction field. This effect can be attributed to a permanent component in the magnetization of the spheres, probably linked to their fabrication process. Because the individual particles rotate at a different speed than the whole assembly, this motion introduces multiple length and time scales which is the ideal condition for efficient micromixing. Secondly, dipole-dipole interactions can contribute to a reorientation. This effect is expected to dominate at vanishing excitation frequency, where the rotation of the assembly is locked with the external field. Furthermore, dipoles on a regular polygon have no reason to rotate as their interaction energy is independent on the orientation with regards to the field. This means that this effect requires a deformation of the assembly. Under high frequency and high amplitude, a very fast rotation is observed that cannot be explained by the previous two mechanisms. In this regime characterized by high deformations, the rotation speed is proportional to the area of the cycles in the configuration space, similarly to how nonreciprocal deformations power the translational motion of microswimmers. This fast mode is particularly suitable for propulsion, and could serve as the basis for efficient ciliate-like locomotion in more complex geometries.

Conflicts of interest

There are no conflicts of interest to declare.

Acknowledgements

This work was financially supported by FNRS PDR grant T.0129.18, as well as the DFG Priority Programme SPP 1726 “Microswimmers - From Single Particle Motion to Collective Behaviour” . GG thanks FRIA for financial support. Computational resources have been provided by the Consortium des Équipements de Calcul Intensif (CÉCI), funded by the Fonds de la Recherche Scientifique de Belgique (F.R.S.-FNRS) under Grant No. 2.5020.11.

Notes and references

• Whitesides and Grzybowski 2002

G. M. Whitesides and B. Grzybowski, Science, 2002, 295, 2418–2421.

• Whitesides and Boncheva 2002

G. M. Whitesides and M. Boncheva, Proc. Natl. Acad. Sci. USA, 2002, 99, 4769–4774.

• Zhang et al. 2009

J. Zhang, Z. Sun and B. Yang, Curr. Opin. Colloid Interface Sci., 2009, 14, 103 – 114.

• Thorkelsson et al. 2015

K. Thorkelsson, P. Bai and T. Xu, Nano Today, 2015, 10, 48 – 66.

• Boles et al. 2016

M. A. Boles, M. Engel and D. V. Talapin, Chem. Rev., 2016, 116, 11220–11289.

• Dou et al. 2017

Y. Dou, B. Wang, M. Jin, Y. Yu, G. Zhou and L. Shui, J. Micromechanics Microengineering, 2017, 27, 113002.

• Rida and Gijs 2004

A. Rida and M. A. M. Gijs, Anal. Chem., 2004, 76, 6239–6246.

• Majedi et al. 2013

F. S. Majedi, M. M. Hasani-Sadrabadi, S. Hojjati Emami, M. A. Shokrgozar, J. J. VanDersarl, E. Dashtimoghadam, A. Bertsch and P. Renaud, Lab Chip, 2013, 13, 204–207.

• Wu et al. 2008

L. Y. Wu, D. Di Carlo and L. P. Lee, Biomed. Microdevices, 2008, 10, 197–202.

• Martinez-Pedrero and Tierno 2015

F. Martinez-Pedrero and P. Tierno, Phys. Rev. Applied, 2015, 3, 051003.

• Chung et al. 2008

S. E. Chung, W. Park, S. Shin, S. A. Lee and S. Kwon, Nature mater., 2008, 7, 581.

• Grosjean et al. 2015

G. Grosjean, G. Lagubeau, A. Darras, M. Hubert, G. Lumay and N. Vandewalle, Sci. Rep., 2015, 5, 16035.

• Kim et al. 2011

S.-H. Kim, S. Y. Lee, S.-M. Yang and G.-R. Yi, NPG Asia Mater., 2011, 3, 25.

• Edel et al. 2013

J. B. Edel, A. A. Kornyshev and M. Urbakh, ACS Nano, 2013, 7, 9526–9532.

• Vandewalle et al. 2012

N. Vandewalle, L. Clermont, D. Terwagne, S. Dorbolo, E. Mersch and G. Lumay, Phys. Rev. E, 2012, 85, 041402.

• Vandewalle et al. 2013

N. Vandewalle, N. Obara and G. Lumay, Eur. Phys. J. E, 2013, 36, 1–6.

• Golosovsky et al. 1999

M. Golosovsky, Y. Saado and D. Davidov, Appl. Phys. Lett., 1999, 75, 4168–4170.

• Wen et al. 2000

W. Wen, L. Zhang and P. Sheng, Phys. Rev. Lett., 2000, 85, 5464–5467.

• Grzybowski et al. 2000

B. A. Grzybowski, H. A. Stone and G. M. Whitesides, Nature, 2000, 405, 1033–1036.

• Golosovsky et al. 2002

M. Golosovsky, Y. Saado and D. Davidov, Phys. Rev. E, 2002, 65, 061405.

• Kralchevsky and Nagayama 1994

P. A. Kralchevsky and K. Nagayama, Langmuir, 1994, 10, 23–36.

• Vella and Mahadevan 2005

D. Vella and L. Mahadevan, Am. J. Phys., 2005, 73, 817–825.

• Ninham 1999

B. Ninham, Adv. Colloid Interface Sci., 1999, 83, 1 – 17.

• Lumay et al. 2013

G. Lumay, N. Obara, F. Weyer and N. Vandewalle, Soft Matter, 2013, 9, 2420–2425.

• Grosjean et al. 2018

G. Grosjean, M. Hubert and N. Vandewalle, Adv. Colloid Interface Sci., 2018, 255, 84 – 93.

• Grosjean et al. 2016

G. Grosjean, M. Hubert, G. Lagubeau and N. Vandewalle, Phys. Rev. E, 2016, 94, 021101.

• Pande et al. 2017

J. Pande, L. Merchant, T. Krüger, J. Harting and A.-S. Smith, New J. Phys., 2017, 19, 053024.

• Sukhov et al. 2019

A. Sukhov, S. Ziegler, Q. Xie, O. Trosman, J. Pande, G. Grosjean, M. Hubert, N. Vandewalle, A.-S. Smith and J. Harting, J. Chem. Phys., 2019, 151, 124707.

• Grosjean et al. 2018

G. Grosjean, M. Hubert, Y. Collard, S. Pillitteri and N. Vandewalle, Eur. Phys. J. E, 2018, 41, 137.

• Snezhko et al. 2006

A. Snezhko, I. S. Aranson and W.-K. Kwok, Phys. Rev. E, 2006, 73, 041306.

• Snezhko and Aranson 2011

A. Snezhko and I. S. Aranson, Nat. Mater., 2011, 10, 698–703.

• Nakata et al. 2013

S. Nakata, M. Hata, Y. S. Ikura, E. Heisler, A. Awazu, H. Kitahata and H. Nishimori, J. Phys. Chem. C, 2013, 117, 24490–24495.

• Bassik et al. 2008

N. Bassik, B. T. Abebe and D. H. Gracias, Langmuir, 2008, 24, 12158–12163.

• Bormashenko et al. 2015

E. Bormashenko, Y. Bormashenko, R. Grynyov, H. Aharoni, G. Whyman and B. P. Binks, J. Phys. Chem. C, 2015, 119, 9910–9915.

• Okawa et al. 2009

D. Okawa, S. J. Pastine, A. Zettl and J. M. Fréchet, J. Am. Chem. Soc., 2009, 131, 5396–5398.

• Chinomona et al. 2015

R. Chinomona, J. Lajeunesse, W. H. Mitchell, Y. Yao and S. E. Spagnolie, Soft Matter, 2015, 11, 1828–1838.

• Lushi and Vlahovska 2015

E. Lushi and P. M. Vlahovska, J. Nonlinear Sci., 2015, 25, 1111–1123.

• Hubert et al. 2013

M. Hubert, G. Grosjean, Y.-E. Corbisier, G. Lumay, F. Weyer, N. Obara and N. Vandewalle, arXiv preprint arXiv:1310.3094.

• Lagubeau et al. 2016

G. Lagubeau, G. Grosjean, A. Darras, G. Lumay, M. Hubert and N. Vandewalle, Phys. Rev. E, 2016, 93, 053117.

• Vella 2015

D. Vella, Annu. Rev. Fluid Mech., 2015, 47, 115–135.

• Purcell 1977

E. M. Purcell, Am. J. Phys., 1977, 45, 3–11.

• Lauga and Powers 2009

E. Lauga and T. R. Powers, Rep. Prog. Phys., 2009, 72, 096601.

• Maestro et al. 2018

A. Maestro, N. Bruot, J. Kotar, N. Uchida, R. Golestanian and P. Cicuta, Commun. Phys., 2018, 1, 28.

• Meng et al. 2019

F. Meng, D. Matsunaga, J. Yeomans and R. Golestanian, Soft Matter, 2019, 15, 3864–3871.

• Belardi et al. 2011

J. Belardi, N. Schorr, O. Prucker and J. Rühe, Adv. Funct. Mater., 2011, 21, 3314–3320.

• den Toonder and Onck 2013

J. M. den Toonder and P. R. Onck, Trends Biotechnol., 2013, 31, 85–91.

• Vilfan et al. 2010

M. Vilfan, A. Potočnik, B. Kavčič, N. Osterman, I. Poberaj, A. Vilfan and D. Babič, Proc. Natl. Acad. Sci. USA, 2010, 107, 1844–1847.

• Kokot et al. 2011

G. Kokot, M. Vilfan, N. Osterman, A. Vilfan, B. Kavčič, I. Poberaj and D. Babič, Biomicrofluidics, 2011, 5, 034103.

• Landau and Lifshitz 1935

L. D. Landau and E. Lifshitz, Phys. Z. Sowjet., 1935, 8, 153.

• Gilbert 2004

T. L. Gilbert, IEEE Trans. Magn., 2004, 40, 3443–3449.

• Vansteenkiste et al. 2014

A. Vansteenkiste, J. Leliaert, M. Dvornik, M. Helsen, F. Garcia-Sanchez and B. Van Waeyenberge, AIP Adv., 2014, 4, 107133.

• Beleggia et al. 2006

M. Beleggia, M. De Graef and Y. Millev, J. Phys. D, 2006, 39, 891.

• Coey 2010

J. M. Coey, Magnetism and magnetic materials, Cambridge university press, 2010.

• Golestanian and Ajdari 2008

R. Golestanian and A. Ajdari, Phys. Rev. E, 2008, 77, 036308.

• Blake and Chwang 1974

J. R. Blake and A. T. Chwang, J. Eng. Math., 1974, 8, 23–29.

• Mazur and van Saarloos 1982

P. Mazur and W. van Saarloos, Physica A, 1982, 115, 21.

• Dörr et al. 2016

A. Dörr, S. Hardt, H. Masoud and H. A. Stone, J. Fluid Mech., 2016, 790, 607–618.