Optimal motion of triangular magnetocapillary swimmers

TL;DR

This paper numerically investigates the optimal motion of triangular magnetocapillary swimmers at liquid interfaces, revealing how external magnetic field frequency and inter-particle distance influence swimmer speed, supported by simulations and bead-spring models.

Contribution

It introduces a hybrid simulation approach to analyze magnetocapillary swimmer dynamics and identifies optimal frequencies for maximum speed based on particle relaxation times.

Findings

Maximum swimmer speed depends sharply on magnetic field frequency.

Optimal frequency is centered around the particle's relaxation time.

Increasing inter-particle distance reduces maximum speed.

Abstract

A system of ferromagnetic particles trapped at a liquid-liquid interface and subjected to a set of magnetic fields (magnetocapillary swimmers) is studied numerically using a hybrid method combining the pseudopotential lattice Boltzmann method and the discrete element method. After investigating the equilibrium properties of a single, two and three particles at the interface, we demonstrate a controlled motion of the swimmer formed by three particles. It shows a sharp dependence of the average center-of-mass speed on the frequency of the time-dependent external magnetic field. Inspired by experiments on magnetocapillary microswimmers, we interpret the obtained maxima of the swimmer speed by the optimal frequency centered around the characteristic relaxation time of a spherical particle. It is also shown that the frequency corresponding to the maximum speed grows and the maximum average speed decreases with increasing inter-particle distances at moderate swimmer sizes. The findings of our lattice Boltzmann simulations are supported by bead-spring model calculations.

Figures

a, suggesting that this could be the relevant process dominating the dynamics. However, our test simulations with a significantly larger box do not reproduce this result, while the coasting time scale reappears as the relevant one in the system.

a, suggesting that this could be the relevant process dominating the dynamics. However, our test simulations with a significantly larger box do not reproduce this result, while the coasting time scale reappears as the relevant one in the system.

Equations

$$f_i^c(\mathbf{x}+{\mathbf{c}}_i\Delta t,t+\Delta t)=f_i^c(\mathbf{x},t)+{\Omega }_i^c(\mathbf{x},t),$$

$${\Omega }_i^c(\mathbf{x},t)=-\frac{f_i^c(\mathbf{x},t)-f_i^{\mathrm{eq}}({\rho }^c(\mathbf{x},t),{\mathbf{u}}^c(\mathbf{x},t))}{{\tau }^c/\Delta t}.$$

$${\mathbf{F}}_{\mathrm{C}}^c(\mathbf{x},t)=-{\psi }^c(\mathbf{x},t)\sum _{c^{\prime }}{g}_{c{c}^{\prime }}\sum _{{\mathbf{x}}^{\prime }}{\psi }^{c^{\prime }}({\mathbf{x}}^{\prime },t)({\mathbf{x}}^{\prime }-\mathbf{x}).$$

$$f_i^c(\mathbf{x}+{\mathbf{c}}_i,t+1)=f_{\bar{i}}^c(\mathbf{x}+{\mathbf{c}}_i,t)+{\Omega }_{\bar{i}}^c(\mathbf{x}+{\mathbf{c}}_i,t)+C\text{mbox},$$

$${\mathbf{F}}_p$$

$${\mathbf{D}}_p$$

$${\mathbf{F}}_{pn}=-\sum _c{\rho }^c(\mathbf{x},t){\mathbf{u}}^c(\mathbf{x},t)\text{mbox}.$$

$${\bar{\rho }}^c(\mathbf{x},t)=\frac{1}{N_{FN}}\sum _{i_{FN}}{\rho }^c(\mathbf{x}+{\mathbf{c}}_{i_{FN}},t)\text{mbox},$$

$$f_i^c(\mathbf{x},t)={\rho }_{\mathrm{new}}^c\cdot f_i^{eq}({\mathbf{u}}_{\mathrm{surface}}(\mathbf{x},t),{\rho }_{\mathrm{new}}(\mathbf{x},t))\text{mbox},$$

$$U_{ij}=-\frac{{\mu }_0}{4\pi r_{ij}^3}\left[3({\mathbf{\mu }}_i\cdot {\mathbf{e}}_{ij})({\mathbf{\mu }}_j\cdot {\mathbf{e}}_{ij})-({\mathbf{\mu }}_i\cdot {\mathbf{\mu }}_j)\right].$$

$${\mathbf{B}}_i=-\frac{\delta U_{ij}}{\delta {\mathbf{\mu }}_i}=\frac{{\mu }_0}{4\pi r_{ij}^3}\left[3{\mathbf{e}}_{ji}({\mathbf{\mu }}_j\cdot {\mathbf{e}}_{ji})-{\mathbf{\mu }}_j\right].$$

$${\mathbf{F}}_i=\frac{3{\mu }_0}{4\pi r_{ij}^4}$$

$${\mathbf{T}}_i=\frac{{\mu }_0}{4\pi r_{ij}^3}\cdot \left(3\left({\mathbf{\mu }}_j\cdot {\mathbf{e}}_{ji}\right)\left[{\mathbf{\mu }}_j\times {\mathbf{e}}_{ji}\right]-\left[{\mathbf{\mu }}_i\times {\mathbf{\mu }}_j\right]\right)+\left[{\mathbf{\mu }}_i\times \mathbf{B}\right].$$

$$\begin{array}{rl}\rho {}_{\mathrm{p}}g& \frac{4}{3}\pi R^3=2\pi \gamma R\sin \phi {}_{\mathrm{c}}\sin \psi +\pi R^3(\rho {}_{fl1}-\rho {}_{fl2})g\times \\ & \left(\frac{2}{3}+\cos \phi {}_{\mathrm{c}}-\frac{1}{3}{\cos }^3\phi {}_{\mathrm{c}}+\frac{h{}_{\mathrm{cl}}}{R}{\sin }^2\phi {}_{\mathrm{c}}\right)+\rho {}_{fl2}g\frac{4}{3}\pi R^3,\end{array}$$

$$Bo=\frac{(\rho {}_{\mathrm{p}}-\rho {}_{fl1})g{R}^2}{\gamma },$$

$${\cos }^3\phi {}_{\mathrm{c}}-3\cos \phi {}_{\mathrm{c}}-\frac{3}{2k_{0}B{o}^{\mathrm{LB}}}\sin (2\phi {}_{\mathrm{c}})+\frac{4}{k_0}-2=0,$$

$$U^{\mathrm{tot}}(r\ll l{}_{\mathrm{c}})=C^{\mathrm{cap}}\ln \left(\frac{r-a}{l{}_{\mathrm{c}}}\right)+\frac{C^{\mathrm{mag}}}{(r-a{)}^3},$$

$$\begin{array}{rl}\frac{L_{12}}{2R}(r\ll l{}_{\mathrm{c}})=& 1+\frac{1}{2R}\sqrt[3]{\frac{3C^{\mathrm{mag}}}{C^{\mathrm{cap}}}}=\\ & 1+\frac{1}{2}\sqrt[3]{\frac{\left(\frac{3{\mu }_0{\mu }_{\mathrm{S}}^2}{4\pi R^4}\right)}{(2\pi \gamma R)(B{o}^{\mathrm{LB}}{)}^2{\sigma }^2}}{\left(\frac{{\mu }_i}{\mu {}_{\mathrm{s}}}\right)}^{2/3},\end{array}$$

$$\frac{L_{12}}{2R}(r\ll l{}_{\mathrm{c}})=1+\frac{1}{2}{\left(B{o}^{\mathrm{LB}}\sigma \right)}^{-2/3}{\left(\frac{{\mu }_i}{\mu {}_{\mathrm{s}}}\right)}^{2/3}.$$

$$U^{\mathrm{tot}}(r\approx l{}_{\mathrm{c}})=-0.45\frac{l{}_{\mathrm{c}}{C}^{\mathrm{cap}}}{r-a}+\frac{C^{\mathrm{mag}}}{(r-a{)}^3}.$$

$$\begin{array}{rl}\frac{L_{12}}{2R}(r\approx l{}_{\mathrm{c}})=1+0.68\frac{\sqrt[4]{\rho {}_{\mathrm{b}}-\rho {}_{\mathrm{r}}}}{2\rho {}_{\mathrm{p}}-\rho {}_{\mathrm{r}}}& {\left(\rho {}_{\mathrm{p}}-\rho {}_{\mathrm{r}}\right)}^{3/4}\times \\ & {\left(B{o}^{\mathrm{LB}}\right)}^{-3/2}\left(\frac{{\mu }_i}{\mu {}_{\mathrm{s}}}\right).\end{array}$$

$$⟨\mathbf{v}⟩=\frac{1}{3}\sum _{i=1}^3\frac{{\mathbf{r}}_i(t+n2\pi /\omega )-{\mathbf{r}}_i(t)}{n2\pi /\omega }=\frac{1}{3}\sum _{i=1}^3{\mathbf{v}}_i,$$

$${\tau }_{\mathrm{HO}}=\sqrt{\frac{m_{\mathrm{R}}}{k}}.$$

$${\omega }_{\mathrm{St}}(L)\sim \frac{1}{L}\sqrt{C_1+C_{2}L+C_3{L}^2},$$

$${\tau }_{\mathrm{St}}=\frac{6\pi \eta R}{k}.$$

$${\tau }_{\mathrm{cs}}=\frac{m}{6\pi \eta R},$$

$${\tau }_{\mathrm{int}}=\frac{R\eta }{\gamma },$$

$${\varphi }_{ij}\equiv \varphi ({\mathbf{r}}_i-{\mathbf{r}}_j)=\frac{1}{2}k(|{\mathbf{r}}_i-{\mathbf{r}}_j|-L{)}^2.$$

$$\hat{\mathbf{T}}({\mathbf{r}}_j-{\mathbf{r}}_k)=\frac{1}{8\pi \eta |{\mathbf{r}}_j-{\mathbf{r}}_k|}\left(\hat{\mathbf{I}}+\frac{({\mathbf{r}}_j-{\mathbf{r}}_k)\otimes ({\mathbf{r}}_j-{\mathbf{r}}_k)}{({\mathbf{r}}_j-{\mathbf{r}}_k{)}^2}\right)$$

$$\begin{array}{rl}\frac{d{\mathbf{r}}_j}{dt}& ={\nu }_j\left({\mathbf{F}}_j(t)+\sum _{k\ne j}\mathbf{G}({\mathbf{r}}_j-{\mathbf{r}}_k)\right)\\ & +\sum _{k\ne j}\hat{\mathbf{T}}({\mathbf{r}}_j-{\mathbf{r}}_k)\left({\mathbf{F}}_k(t)+\sum _{l\ne k}\mathbf{G}({\mathbf{r}}_k-{\mathbf{r}}_l)\right),\end{array}$$

Full text

††preprint: AIP/123-QED

I Introduction

Microswimmers are a paradigmatic example of active matter, and have come into focus in the last two decades with the development of non-equilibrium physics Lauga and Powers (2009). Typically, microswimmers are associated with motion at low Reynolds numbers, when self-propulsion dominates over stochastic diffusion, at least in the case of artificial devices. Actually, because self-propulsion requires breaking of the time-reversal symmetry in the stroke, the design of microswimmers often involves a compromise between engineering a device with the smallest possible degrees of freedom (at least two) and obtaining a time-irreversible stroke Elgeti, Winkler, and Gompper (2015). These difficulties are circumvented in a number of devices Purcell (1977); Avron, Kenneth, and Oaknin (2005); Downton and Stark (2009), however, the devices based on interacting beads gained particular attention starting with the analysis of a linear arrangement of three spherical beads by Najafi and Golestanian Najafi and Golestanian (2004). In this design, the stroke involves a phase-shift between the contractions and elongations of the two arms connecting the central bead and the external ones, propelling the swimmer along the line Golestanian and Ajdari (2008). This model was used to identify some of the fundamental properties of microswimming including the difference between pushers and pullers Pande and Smith (2015); Daddi-Moussa-Ider et al. (2018), interaction between two swimmers Pooley, Alexander, and Yeomans (2007), and of the swimmer with a wall Zargar, Najafi, and Miri (2009); Daddi-Moussa-Ider et al. (2018).

An advantage of the bead-based design is its amenability to experimental investigations. For example, a linear swimmer was achieved employing three ferromagnetic beads placed on a surface of water, which due to the balancing of repulsive magnetic interactions between the induced dipole moments and the capillary attraction, adopt a metastable, linear configuration Grosjean et al. (2016). The oscillation of the arms is generated by oscillating magnetic fields, such that the overall force applied to the device remains zero on average Lagubeau et al. (2016). This design is fundamentally different to the original theoretical counterpart in that the stroke is not imposed but emerges from the balance of forces on the bead, an effect that was concomitantly incorporated into the theoretical modelling and simulations Felderhof (2006); Pande and Smith (2015); Pickl et al. (2017); Pande et al. (2017a, b).

A more stable configuration is a triangular arrangement of magnetic particles at the interface, a real experimental situation of which is demonstrated in Fig. 1Grosjean et al. (2015). The particles form an equilateral triangle at an air-water interface. The size of the triangular arrangement is controlled by an external vertical static magnetic field, while the swimmer is driven by a time-dependent magnetic field aligned in the plane of the swimmer.

Numerous efforts involving bead-based devices focused on triangular arrangements that allow for the rotation of the swimmer Rizvi, Farutin, and Misbah (2018), and a richer interaction with external flows Küchler, Löwen, and Menzel (2016). The first models were concerned with modelling the swimming stroke of Chlamydomonas reinhardtii Friedrich and Jülicher (2012); Polotzek and Friedrich (2013); Bennett and Golestanian (2013); Jibuti et al. (2017). Alternative propulsion mechanisms included modelling a stroke in which the angle between the extendable arms is fixed Ledesma-Aguilar, Löwen, and Yeomans (2012), and a stroke in which both arm lengths and the enclosed angle vary Earl et al. (2007). Further efforts involve a design imposing motion of three beads on a circle Dreyfus, Baudry, and Stone (2005) or rotating beads Lobaskin, Lobaskin, and Kulić (2008). In the context of force-based description, internal driving was recently considered based on internal forces acting along the lines connecting the beads arranged in a triangular geometry Rizvi, Farutin, and Misbah (2018).

Although numerous theoretical studies are known for swimmers in the form of linear Felderhof (2006); Zargar, Najafi, and Miri (2009); Pickl et al. (2012); Taghiloo and Miri (2013); Pande and Smith (2015); Babel, Löwen, and Menzel (2016) or equilateral arrangements Chinomona et al. (2015), we are not aware of any simulations that include a rigorous description of the corresponding hydrodynamic flows, the time dependent deformation of the interface and thus the related capillary effects. It turns out to be rather challenging to analytically take the partial immersion of the swimmer particles in the liquid, particle rotations, inter-particle forces and their surface effects into account Chinomona et al. (2015). In the present study we aim at bridging this gap by employing the lattice Boltzmann (LB) method for the dynamics of the fluids coupled to a discrete element method for the particle dynamics. We perform extensive numerical simulations of triangular magnetocapillary swimmers consisting of three ferromagnetic beads. Our LB-simulations contain the full translational and rotational dynamics of the particles, the dynamics of the fluid flow and can provide insight into the fluid-particle interaction. Combining our simulations with a force driven bead spring model of a triangular swimmer, we explore the relation between propulsion velocity and the frequency of driving.

The remainder of this article is organized as follows: in the following section we present the numerical method which is followed by a systematic investigation of the equilibrium particle configuration in the presence of capillary and magnetic interactions. Then, we demonstrate our results of the dynamic behaviour of the triangular magnetocapillary swimmer. Subsequently, we analyze our simulation results with the help of a semianalytical model and summarize our findings.

II Simulation method and setup

For the simulation of the fluids we use a lattice Boltzmann (LB) method Benzi, Succi, and Vergassola (1992). The LB method allows a straightforward implementation of complex boundary conditions and due to the locality of the algorithm it is well suited for the implementation on parallel supercomputers. It is based on a discretized version of the Boltzmann equation

[TABLE]

which describes the time evolution of a single-particle distribution function $f^{c}_{i}(\bm{x},t)$ at time $t$ and position $\bm{x}$. $\bm{c}_{i}$ denotes the discrete velocity vector in the $i$th direction for fluid component $c=\{1,2\}$. Here, we use a so-called D3Q19 lattice with $i=1,\ldots,19$Qian, D’Humières, and Lallemand (1992). The left hand side of Eq. (1) describes the free streaming of fluid particles, while their collisions are modelled by a Bhatnagar-Gross-Krook (BGK) collision operator on the right hand side as Bhatnagar, Gross, and Krook (1954)

[TABLE]

Here, $f^{\mathrm{eq}}_{i}(\rho^{c}(\bm{x},t),\bm{u}^{c}(\bm{x},t))$ is a third-order equilibrium distribution function, and macroscopic densities and velocities are given by $\rho^{c}(\bm{x},t)=\rho_{0}\sum_{i}f^{c}_{i}(\bm{x},t)$ and $\bm{u}^{c}(\bm{x},t)=\sum_{i}f^{c}_{i}(\bm{x},t)\bm{c}_{i}/\rho^{c}(\bm{x},t)$, respectively ($\rho_{0}$ is a reference density). $\tau{{}_{\mathrm{c}}}$ is the relaxation rate of component $c$, which determines the relaxation of $f^{c}_{i}(\bm{x},t)$ towards the equilibrium. We employ a three-dimensional lattice with the cell size $\Delta x$ in space and the time $t$ is discretized with $\Delta t$-steps. The speed of sound $c{{}_{\mathrm{s}}}=1/\sqrt{3}\Delta x/\Delta t$ depends on the choice of the lattice geometry and allows one to obtain the kinematic $\nu^{c}=c^{2}_{\mathrm{s}}\Delta t(\tau^{c}/\Delta t-1/2)$ or the dynamic $\eta^{c}=\nu^{c}\rho^{c}$ fluid viscosities. For simplicity, we set $\Delta x=\Delta t=\rho_{0}=\tau^{c}=1$ in the remainder of this paper.

Capillary interactions can be modeled when more than one fluid species is present. In this case a mean-field interaction force between several fluid components is calculated according to the pseudopotential method of Shan and Chen asShan and Chen (1993); Liu et al. (2016)

[TABLE]

Here, $c$ and $c^{\prime}$ refer to different fluid components, $\bm{x}^{\prime}$ denotes the nearest neighbours of the lattice site $\bm{x}$ and $g_{cc^{\prime}}$ describes a coupling constant determining the surface tension. $\psi^{c}(\bm{x},t)$ has the functional form $\psi^{c}(\bm{x},t)\equiv\psi^{c}(\rho^{c}(\bm{x},t))=1-\mathrm{e}^{-\rho^{c}(\bm{x},t)}$. The force (3) is applied to the fluid component $c$ by adding a shift $\Delta\bm{u}^{c}(\bm{x},t)=\tau^{c}\bm{F}^{c}_{\mathrm{C}}(\bm{x},t)/\rho^{c}(\bm{x},t)$ to the velocity $\bm{u}^{c}(\bm{x},t)$ in the equilibrium distribution. The method is a diffuse interface method, with an interface width of typically $5$ lattice sites depending weakly on the coupling strength Frijters, Günther, and Harting (2012); Krüger et al. (2013). In the binary fluid system we refer to the fluids as “red” (r) and “blue” (b) Jansen and Harting (2011).

Rigid particles are simulated by solving Newton’s equations of motion for translational and rotational degrees of freedom by means of a leap-frog algorithm. The particles are discretized on the lattice and are coupled to both fluid species by means of a modified bounce-back boundary condition for multiple fluid components Ladd (1994); Ladd and Verberg (2001); Jansen and Harting (2011); Günther, Frijters, and Harting (2014).

When a lattice site $\bm{x}$ is occupied by the surface of a particle, the following equation is applied to its neighbouring fluid lattice site $(\bm{x}+\bm{c}_{i}$)

[TABLE]

where $C$ is a linear function of the local velocity of particle surface Jansen and Harting (2011) and $\bar{i}$ is defined such that $\bm{c}_{i}=-\bm{c}_{\bar{i}}$. In order to conserve the total momentum of the system, an additional force $\bm{F}_{p}$ and torque $\bm{D}_{p}$ is applied on the particle Jansen and Harting (2011) to compensate for the momentum change of the fluid caused by Eq. (4)

[TABLE]

in which $\bm{r}(t)$ is a vector directed from the particle center to the lattice site of reflection. While the particle moves, the configuration of lattice sites occupied by the particle changes. For newly occupied sites, the fluid on that site is removed and its momentum is added to the particle through a force Jansen and Harting (2011)

[TABLE]

Vacated lattice sites need to be filled with fluid. In the case of two fluid components, we define an average density Jansen and Harting (2011); Frijters, Günther, and Harting (2012)

[TABLE]

where $N_{FN}$ is the number of neighbouring fluid sites with coordinates $\bm{x}_{i_{FN}}=\bm{x}+\bm{c}_{i_{FN}}$. The fluid on the vacated site is then initialised with distribution functions Jansen and Harting (2011); Frijters, Günther, and Harting (2012)

[TABLE]

where $\bm{u}_{\mathrm{surface}}(\bm{x},t)$ is the local velocity of the particle surface. $\rho^{c}_{\mathrm{new}}(\bm{x},t)$ corresponds to $\bar{\rho}^{c}(\bm{x},t)$ plus a small correction term to account for local density gradients imposed by the Shan-Chen forcingJansen and Harting (2011); Frijters, Günther, and Harting (2012).

When two particles approach each other very closely so that no lattice node is available between their surfaces to describe the hydrodynamic forces, a lubrication correction can be applied Ladd and Verberg (2001). However, we always assure a sufficient resolution and the particles never get closer than a few $\Delta x$ in the current paper. Therefore, such a correction is not required here.

The ferromagnetic particles are considered as rigid spheres with fixed orientation of the magnetic moment with respect to the particle. When applying a static magnetic field $B_{\mathrm{y}}$ along the positive $y$-direction (Fig. 2), repulsive dipolar forces are induced. The magnetic repulsion is balanced by an attractive capillary force which is due to the interface deformation caused by the gravity-induced immersion of the particles. This leads to a stable arrangement of the beads after a certain relaxation time. In analogy with the experiments on magnetocapillary swimmers Lumay et al. (2013); Grosjean et al. (2015) we choose the amplitude of the time-dependent magnetic field to be approximately three times lower than that of the static field to treat it as a modulation. The field $\bm{B}(t)=B_{0\mathrm{x}}\cos\omega t\bm{e}_{\mathrm{x}}$ drives the system out of the local equilibrium leading to a directed collective motion of the beads.

In our implementation we assume a homogeneous external magnetic field $\bm{B}$ and prescribe a magnetic moment $\bm{\mu}_{i}=\chi V\bm{B}/\mu_{0}$ to each particle $i$, where $\chi$ is the particle susceptibility, $V$ is its volume and $\mu_{0}=4\pi\times 10^{-7}$ (in lattice units) corresponds to the magnetic permeability of vacuum. The emergent magnetic dipole-dipole interaction between a pair of particles is then

[TABLE]

Here, $r_{ij}\equiv\left|\left|\bm{r}_{ij}\right|\right|\equiv\left|\left|\bm{r}_{i}-\bm{r}_{j}\right|\right|$ is the distance between the centres of two spheres $i,j$ located at $\bm{r}_{i}$ and $\bm{r}_{j}$, respectively, and $\bm{e}_{ij}$ = $(\bm{r}_{i}-\bm{r}_{j})/\left|\left|\bm{r}_{i}-\bm{r}_{j}\right|\right|$. The magnetic field generated by the magnetic moment $\bm{\mu}_{j}$ at the location of another particle $i$ is

[TABLE]

The resulting magnetic force acting on the $i$th particle is then $\bm{F}_{i}=-\bm{\nabla}\left(-\bm{\mu}_{i}\cdot(\bm{B}_{i}+\bm{B})\right)$, or more explicitly

[TABLE]

Similarly, the magnetic torque acting on the particle $i$ is $\bm{T}_{i}=\left[\bm{\mu}_{i}\times(\bm{B}_{i}+\bm{B})\right]$, or explicitly

[TABLE]

The total force and the total torque for each particle include a summation of expressions (12) and (13) over index $j$. We note that the homogeneous external magnetic field $\bm{B}$ does not produce a net force on the particles and serves exclusively to change the direction of the particle magnetic moment. The magnetic forces and the torques given by Eqs. (12) and (13) are chosen such that they balance capillary forces in the equilibrium. The method with implemented magnetic interactions has already been benchmarked and successfully applied for simulations of magnetocapillary phenomena Xie et al. (2015); Xie, Davies, and Harting (2016, 2017).

We consider a simulation box consisting of $128^{3}$ cubic cells containing two equally-sized fluid lamellae (Fig. 2). Walls with mid-grid bounce back boundary conditions are placed parallel to the interface at $y=0$ and $y=128$, whereas for the boundaries in all other directions we assume periodicity. All particles have equal radii $R=5$ which assures that they are small on the scale of the simulation box, but sufficiently large with respect to the interface thicknessXie et al. (2015); Xie, Davies, and Harting (2016, 2017). The beads are initially placed near the interface between the red and blue fluids with densities $\rho_{r}=\rho_{b}=0.7$, deforming thus the interface of the fluids due to the presence of a body force mimicking gravity acting on the particles. This leads to attractive capillary forces between the beads. The fluid-fluid coupling constant is chosen $g_{cc^{\prime}}=0.1$ giving rise to the numerical surface tension $\gamma\approx 0.04$ in LB-units (see Ref. Frijters, Günther, and Harting (2012)), which assures a well-defined interface profile on the one hand, and the numerical stability of the method on the other hand. In all simulations presented here, the two fluids are considered identical, i.e. with equal viscosity and density. Due to the limitations of the pseudopotential lattice Boltzmann method, we are not able to reach realistic values for the surface tension and density ratios to mimic the air and water phases of the experiment. Even though we keep the dimensionless parameters such as the Bond number close to the experimental values wherever possible, we do expect in particular the limitation in the choice of surface tensions to have an impact on the swimmer propulsion. The ratio of surface tension forces and magnetic forces for parameters used in the experiment is several orders of magnitude higher than in the simulations. In addition, even density ratios of the order of 10 or 20 comprise the numerical stability resulting in a smaller available range of driving frequencies and Bond numbers. The long-range nature of capillary interactions would require system sizes of at least $10$ times the size of the swimmer in order to remove finite size effects from the simulations Xie, Davies, and Harting (2016); C. Kunert, J. Harting and O.I. Vinogradova (lack); C. Kunert and J. Harting (lack). However, due to the inherently slow motion of the swimmer, our simulations require millions of time steps and thus the computational time required would be prohibitive. When checking the impact of the finite size effects on the maximum swimmer speed, we do find only a factor two increase in the maximum swimming speed even for ten times larger systems and we do not expect the qualitative behaviour of the swimmer to change substantially. Thus, we expect quantitative differences in the propulsion of the experimental swimmer as compared to our simulations and therefore refrain from a quantitative comparison. For a quantitative comparison, an improved multiphase model would be required.

III Equilibrium position of magnetic particles at a fluid-fluid interface

III.1 Vertical equilibrium for a single particle

Prior to the simulation of multiple particles, we focus first on the equilibrium properties of a single particle at the interface. In the case of a single bead the homogeneous magnetic field exerts no magnetic force. The particle reaches its equilibrium and floats when the gravitational force $\bm{F}{{}_{\mathrm{g}}}$ is balanced by the sum of the surface tension force $\bm{F}{{}_{\mathrm{st}}}$ (or its vertical projection) and the force related to the hydrostatic pressure $\bm{F}{{}_{\mathrm{hp}}}$ (Fig. 3a), i.e. $F_{\mathrm{g}}^{\mathrm{y}}=F_{\mathrm{st}}^{\mathrm{y}}+F_{\mathrm{hp}}^{\mathrm{y}}$. More explicitly, the force balance perpendicular to the interface reads Binks and Horozov (2006); Vella and Mahadevan (2005)

[TABLE]

where $g$ denotes the gravitational constant and $\rho{{}_{\mathrm{p}}},\rho{{}_{\mathrm{fl1}}},\rho{{}_{\mathrm{fl2}}}$ are the densities of the particle and the fluids 1 and 2, respectively. The contact angle $\theta{{}_{\mathrm{c}}}$, which is fixed in the present simulations ($\theta{{}_{\mathrm{c}}}=\pi/2$), is the angle between the tangent to the fluid interface and the perpendicular from the particle center of mass, both in the contact point. $\varphi{{}_{\mathrm{c}}}$ is defined as the angle between the vertical line and the line from the center of mass to the contact point. Finally, the angle $\psi$ is the difference $\psi=\theta{{}_{\mathrm{c}}}-\varphi{{}_{\mathrm{c}}}$. $h{{}_{\mathrm{c}}}$ is the immersion depth of the particle and $h{{}_{\mathrm{cl}}}$ is the height of the perturbed liquid relative to the contact line (Fig. 3a).

Based on Eq. (14) it is straightforward to introduce a dimensionless quantity, the Bond number, which relates the downward- and upward forces

[TABLE]

where it is assumed that the interface is flat ($h{{}_{\mathrm{cl}}}\rightarrow 0$) and the particle is almost completely immersed in the fluid ($\varphi{{}_{\mathrm{c}}}\approx 0$). Given the parameters of the experimental situation from Refs. Lumay et al. (2013); Grosjean et al. (2015), i.e. particles with a radius of $R^{\mathrm{exp}}=0.25\cdot 10^{-3}$ m, a density of $\rho{{}_{\mathrm{p}}}=7.8\cdot 10^{3}$ kg/m3 ($\rho{{}_{\mathrm{g}}}\approx 0$) and the surface tension of water $\gamma=73$ mN/m, we obtain for the experimental Bond numbers $Bo^{\mathrm{exp}}\approx 0.06$.

In the present LB-simulations for equal densities of the red and blue fluids $\rho{{}_{\mathrm{fl1}}}=\rho{{}_{\mathrm{fl2}}}=0.7$ one obtains $Bo^{\mathrm{LB}}=(\rho{{}_{\mathrm{p}}}-\rho{{}_{\mathrm{fl2}}})gR^{2}/\gamma$. The dependence of $h{{}_{\mathrm{c}}}(Bo)$ is shown in Fig. 3b. It is linear up to the point where the particle detaches from the interface. The LB-simulations are compared with the solution of the equation for the angle $\varphi{{}_{\mathrm{c}}}$ obtained from Eq. (14) for $h{{}_{\mathrm{cl}}}\approx 0$ as

[TABLE]

where $k_{0}=\frac{(\rho{{}_{\mathrm{fl1}}}-\rho{{}_{\mathrm{fl2}}})}{(\rho{{}_{\mathrm{p}}}-\rho{{}_{\mathrm{fl2}}})}$ (note that at the interface $\rho{{}_{\mathrm{fl1}}}-\rho{{}_{\mathrm{fl2}}}=0.04\rho_{0}$) and a solution of the form $h{{}_{\mathrm{c}}}\approx 2R\cos\varphi{{}_{\mathrm{c}}}$ is found. As inferred from Fig. 3b the LB-simulations are in good agreement with the solution (16).

In the general case the Bond number for the detachment of the particle from the interface is given by $Bo_{\mathrm{crit}}\approx 3/2\sin^{2}(\theta{{}_{\mathrm{c}}}/2)|_{\theta{{}_{\mathrm{c}}}=\pi/2}\approx 0.75$Vella, Lee, and Kim (2006). Corrected by the factor due to the presence of the other fluid, the critical Bond number in the simulations is in line with the analytical prediction $Bo^{\mathrm{LB}}_{\mathrm{crit}}\approx 0.21$.

III.2 Equilibrium of two particles

When two ferromagnetic particles are placed at the fluid/fluid interface at zero magnetic field, they experience a capillary attraction due to the deformation of the interface. For a pair of small particles the capillary potential reads $U^{\mathrm{cap}}=-C^{\mathrm{cap}}K_{0}\left(r/l{{}_{\mathrm{c}}}\right)$ (Ref. Vella and Mahadevan (2005)), where $r$ is the distance between the centers of the particles, $K_{0}(x)$ is the modified Bessel function of the second kind, $l{{}_{\mathrm{c}}}=\sqrt{\frac{\gamma}{(\rho{{}_{\mathrm{fl1}}}-\rho{{}_{\mathrm{fl2}}})g}}$ is the capillary length and $C^{\mathrm{cap}}=2\pi\gamma R^{2}(Bo^{\mathrm{LB}})^{2}\sigma^{2}$, with $\sigma=\frac{1}{3}(\frac{\rho{{}_{\mathrm{p}}}}{\rho{{}_{\mathrm{p}}}-\rho{{}_{\mathrm{fl2}}}}+1)$ Vella and Mahadevan (2005). Since the capillary length is dependent on the Bond number, i.e. $l{{}_{\mathrm{c}}}=R\sqrt{\frac{(\rho{{}_{\mathrm{p}}}-\rho{{}_{\mathrm{fl2}}})}{(\rho{{}_{\mathrm{fl1}}}-\rho{{}_{\mathrm{fl2}}})}\frac{1}{Bo^{\mathrm{LB}}}}$, we estimate the minimal capillary length for a high $Bo^{\mathrm{LB}}=0.16$ in the simulations as $\frac{l{{}_{\mathrm{c}}}}{2R}\approx 7$. In contrast, for a low value of $Bo^{\mathrm{LB}}=0.016$ we obtain for $\frac{l{{}_{\mathrm{c}}}}{2R}\approx 22$.

In the situation when the inter-particle distance $r$ is much shorter than the capillary length $r\ll l{{}_{\mathrm{c}}}$, the capillary potential can be well approximated by $U^{\mathrm{cap}}\approx C^{\mathrm{cap}}\ln\left(r/l{{}_{\mathrm{c}}}\right)$. For reaching an equilibrium a uniform static magnetic field aligned perpendicularly to the interface is required, since it turns the dipole-dipole interaction to the repulsive one (Eq. (10), $U_{12}(\uparrow\uparrow)=\mu_{0}\mu_{i}\mu_{j}/(4\pi r^{3}_{ij})$), which opposes the attractive capillary interaction for reaching a finite equilibrium distance. The total potential thus has the form

[TABLE]

where the term $(r-a)$ accounts for the finite size of the beads, $a=2R$ is the minimal distance between the centres of the beads and $C^{\mathrm{mag}}=\mu_{0}\mu_{i}\mu_{j}/(4\pi)$. Potential (17) is schematically illustrated in Fig. 4 and clarifies the role of the first and second derivatives.

After minimization of Eq. (17) one obtains for the equilibrium inter-particle distance

[TABLE]

where we note that the ratio $\left(\frac{3\mu_{0}\mu^{2}_{\mathrm{S}}}{4\pi R^{4}}\right)/(2\pi\gamma R)$ is the ratio of magnetic to capillary forces at equilibrium and should be of order one; $\mu{{}_{\mathrm{S}}}$ denotes here the saturated magnetic moment. The equilibrium position for inter-particle distances well below the capillary length thus scales as

[TABLE]

When $r$ approaches the capillary length ($r\approx l{{}_{\mathrm{c}}}$), the logarithmic approximation for the Bessel function no longer holds. A function which approximates the Bessel function in this range very well is $\sqrt{\pi/2}\mathrm{e}^{-r/l{{}_{\mathrm{c}}}}/\sqrt{r/l{{}_{\mathrm{c}}}}$ (Ref. I.N. Bronshtein, K. A. Semendyayev (lack)). However, an analytical handling of such a function is complex. An alternative approach is to use an empirical expression $U^{\mathrm{cap}}\approx-0.45\frac{C^{\mathrm{cap}}l{{}_{\mathrm{c}}}}{r-a}$ which fits numerically the Bessel function in this range well and allows to write the total potential as

[TABLE]

Minimization of expression (20) under similar assumptions as of Eq. (18) yields for the equilibrium $L_{12}$

[TABLE]

Fig. 5 demonstrates three regimes of the equilibrium $L_{12}$-dependence on the strength of the static magnetic field or the reduced magnetic moment $m_{y}=\mu{{}_{\mathrm{y}}}/\mu{{}_{\mathrm{S}}}$. The three areas correspond to experimental observations shown in Fig. 3 of Ref. Vandewalle et al. (2012). Below $m{{}_{\mathrm{y}}}\approx 0.05$ the equilibrium inter-particle distance becomes "locked". This phenomenon is due to capillary bridges Kralchevsky and K. (2001); Myshkis and V.G. (1987), where a small amount of liquid causes sticking of particles together because of the minimized liquid interface. In the range $0.05<m{{}_{\mathrm{y}}}<0.2$ we observe different growth dependencies for low and high Bond numbers, since the particles are already separated and $L_{12}$ is well below $l{{}_{\mathrm{c}}}$ for $Bo^{\mathrm{LB}}=0.016$ (Eq. (19)) or approximately around $l{{}_{\mathrm{c}}}$ for $Bo^{\mathrm{LB}}=0.16$ (Eq. 21). Finally, $L_{12}$ grows linearly for high magnetic moments which is related to the proximity of $L_{12}$ to the capillary length, as given by Eq. (21).

III.3 Equilibrium of three particles

At last we investigate the case of three particles. All three particles are placed onto the interface simultaneously and a vertical static magnetic field (no time-dependent field so far) is applied from the very beginning. The initial inter-particle distance is set to $\frac{L_{ij}}{2R}(t=0)=2.6$. From symmetry considerations we expect an equilateral triangle in the plane of the interface for the equilibrium configuration. As evidenced by Fig. 6, the length of the sides of the triangle $L_{ij}$ remains equal irrespective of the induced magnetic moment.

An analytical treatment of the problem is more involved, since it requires taking into account the asymmetry of the capillary potential and the finite size of all three particles. Nevertheless, one can clearly observe the "locked" range until $m{{}_{\mathrm{y}}}\approx 0.05$ and the range of the linear growth for $L_{ij}$. As in the case of two particles in Fig. 5, the present simulations fully recover the measurements demonstrated in Fig. 6 or Ref. Vandewalle et al. (2012) Noteworthy is also the fact that the values of the magnetic moment for which the transition from the “locked” to the linear regimes takes place, remain quantitatively unchanged for two and three beads at high $Bo=0.16$ (cf. Figs. 5 and 6). The latter we explain as a compensation of the amount to which the capillary interactions grow by the increased repulsive magnetic interaction all due to the presence of the third particle.

IV Dynamics of a magnetocapillary swimmer

IV.1 Directed motion

So far, only static configurations formed by one, two or three beads were considered. Now we intend to investigate the dynamical behaviour of the three-particle ensemble that performs directed motion powered by an external time-dependent magnetic field. The ensemble is termed a magnetocapillary swimmer.

The simplest stable system which can break the symmetry required by the scallop theorem is a three-beads swimmer having a triangular shape (Fig. 2)Purcell (1977). It takes approximately $30000$ time steps for the $x$-coordinate of all beads to reach the stable configuration (cf. onset of the trajectories in Fig. 7). This time is chosen rather empirically after examining the relaxation dynamics of particles in the LB-simulations. It assures a full vertical and horizontal equilibrium for the particles and the interface. Only after this initial relaxation the oscillating time-dependent in-plane magnetic field $\bm{B}(t)=B_{0\mathrm{x}}\cos\omega t\bm{e}_{\mathrm{x}}$ is applied. For the amplitudes $|B(t)|/|B|\approx 0.36$ we choose a ratio which is close to the one used in the experiments Lumay et al. (2013); Grosjean et al. (2015). In order to stay well below the critical Bond number for the detachment of the particle from the interface and to still assure a significant interface deformation we keep $Bo=0.16$.

For characterization of the swimmer motion we introduce an average velocity of the center of mass in the non-moving frame defined by

[TABLE]

where $n$ counts the number of the periods of the oscillating field and $\bm{r}_{i}$ stands for the coordinates of the beads’ centers.

Fig. 7 demonstrates the motion of the swimmer in the case where the direction of the oscillating field $\bm{B}(t)$ is fixed along the $X$-axis with the fixed given frequency. In this configuration the swimmer propagates itself with minor deviations of all beads from their equilibrium positions and shows a linear motion along the $z$-axis. Careful analysis of all coordinate components reveals that not only do the $Z$-coordinates demonstrate visible oscillations, but also the $Y$-components, i.e. the vertical coordinates, oscillate with amplitudes close to $0.08R$. Finally, since the AC magnetic field is applied along the line connecting particles 1 and 2, and since particle 3 is symmetric to both the field and those particles, one can observe stronger oscillations of the $Z$-coordinates for beads 1 and 2.

Fig. 7 does not answer the question how the particles propagate during the stroke cycle. This issue is clarified in Fig. 8, where the propagation of the orientation vector (director) $\bm{n}_{i}$, coinciding with the magnetization orientation for each bead is demonstrated. As expected, the particles do not spin until the external oscillating magnetic field starts driving them (after $t=30000$), while upon switching the $\bm{B}(t)$-field on, all three particles start oscillating. Note, that the particles do not fully rotate around their own axes and the plane of the director oscillations remains always perpendicular to the net swimmer motion. Although the amplitude of the oscillating field is of the order of one third of the static field amplitude $|B(t)|/|B|\approx 0.36$, the magnetic forces and torques acting on the beads scale as $|F^{\mathrm{magn}}|\sim|B^{2}|$, resulting therefore in $|F^{\mathrm{magn}}(t)|/|F^{\mathrm{magn}}|\approx 0.13$. This explains why the change in the rotational orientation of all the beads is so weak ($n{{}_{\mathrm{y}}}>n{{}_{\mathrm{x}}},n{{}_{\mathrm{z}}}$).

We should point to one significant difference in the propagation mechanism of the swimmer in the experimentsVandewalle et al. (2012); Lumay et al. (2013); Grosjean et al. (2015, 2016); Lagubeau et al. (2016) and the implementation presented above. In the aforementioned experiments a set of three magnetic fields was applied: a strong static field perpendicularly to the interface to induce the magnetic repulsion between the beads, an in-plane oscillating magnetic field to modify the local equilibrium and a weak static in-plane field applied under an angle to the oscillating one to select a unique oscillation modeGrosjean et al. (2015). The action of the three fields led to a strong rotation of the particles in the plane of the interface. The magnetic forces and the torques implemented using Eqs. (12) and (13), respectively, are not supposed to result in any in-plane rotations of the particles and demonstrate small out-of plane oscillations of the beads (Fig. 8). For this reason we expect different speeds of propulsion within the present simulations and the experiments.

IV.2 Frequency dependence of the swimmer motion

The most striking feature of our magnetocapillary swimmer is its response to changes in frequency of the external driving force. While not all cyclic driving protocols yield self-propulsion, if we fix the directions of both the external static and the time-dependent field as shown in Fig. 2, the average speed $\left<v_{\mathrm{z}}\right>$ experiences a sharp peak showing that there is a driving for which the swimmer is most efficient given its geometry. However, the effect of the geometry is non-trivial (Fig. 9e) as the average speed of swimming at the optimum driving frequency first increases (for $L_{ij}<3\times 2R$) and then decays (for $L_{ij}>3\times 2R$).

A similar non-monotonous frequency response was previously reported in the experimentally realized linear magnetocapillary swimmer Grosjean et al. (2016), where, due to the non-Stokesian nature of the device, the swimming speed demonstrates a true resonance at the characteristic frequency of the harmonic oscillator $\tau_{\mathrm{HO}}$ with

[TABLE]

Here, $m_{\mathrm{R}}=m/2$ and $m$ is the mass of the bead and $k$ the curvature of the potential between the beads Lagubeau et al. (2016). Nonetheless, in this case the resonant frequency monotonously decreases with increasing $L_{ij}$. In the case of a triangular geometry, the velocity-frequency relation could not be determined experimentally, while the modelling showed a more complex behaviour as a consequence of the observed inertial periodic rotation of the beads M. Hubert (lack). However, while an interesting direction to pursue, the peak observed in LB simulations (Fig. 9) cannot be associated with $1/\tau_{\mathrm{HO}}$, and actually appears at up to an order of magnitude smaller frequencies.

Notably, an optimal velocity as a function of frequency was also found in triangular bead-spring swimmers in the pure Stokes regime Rizvi, Farutin, and Misbah (2018). Calculation of the swimming velocity within a perturbative approach to the leading order in $R/L$, provides the dependency of $v\sim L^{-2}$. This dependency, as well as the existence of the peak velocity in the frequency domain seems to be insensitive to the geometry of the swimmer and the details of the driving protocol. Namely, this behaviour was also found for the linear swimmer, using the same methodology Pande and Smith (2015); Felderhof (2006); Pande et al. (2017a); Rizvi, Farutin, and Misbah (2018). However, at the same level of theory, the resonant frequency turns out to be independent of the swimmer size Pande and Smith (2015); Pande et al. (2017a); Rizvi, Farutin, and Misbah (2018), despite the fact that our simulations suggest a more complex behaviour (Fig. 9e).

Interestingly, when all orders in $R/L$ are taken into account, for a linear swimmer Felderhof (2006), an increase of the optimum frequency $\omega{{}_{\mathrm{St}}}$ with the inter-particle distance $L$, and $C_{1}$, $C_{2}$, and $C_{3}$ being functions of system parametersFelderhof (2006)

[TABLE]

has been reported. This result, obtained in the linear arrangement of beads, seems to capture the general trends as observed in the simulations.

To check if such a behavior is expected in a triangular geometry, we numerically calculate the swimming velocity of an externally driven, triangular bead-spring device (Appendix Fig. 10), by this avoiding the truncation of the perturbation series in the orders of $R/L$ (see Appendix for details). The model recovers the $\sim L^{-2}$ dependence of the maximum speed (Appendix Fig. 11) Rizvi, Farutin, and Misbah (2018); Pande and Smith (2015); Felderhof (2006); Pande et al. (2017a) as well as the dependence of the optimal frequency as a function of the increasing distance between the beads given by Eq. (24) ($C_{1},C_{2},C_{3}$ are here treated as fitting parameters).

Using dimensional analysis, however, one can show that the emergent velocity of the swimmer in all Stokes perturbative models depends only on $R/L$, $A/(kR)$, and $\omega\tau_{St}$, where $A$ denotes the amplitude of the applied forces. Here, $\tau_{\mathrm{St}}$ is the characteristic time scale of relaxation of the swimmer spring Pande et al. (2017a), relative to the Stokes drag

[TABLE]

In the absence of inertia, $\tau_{\mathrm{St}}$ sets the internal time scale of the swimmer. Unfortunately, $\tau_{\mathrm{St}}$ is not recovered in LB simulations as the time scale is characteristic for the frequency associated with the peak speed, and appears at up to an order of magnitude larger frequencies.

In fact, we find this peak to be characterized by the time scale of bead coasting through the fluid $\tau_{\mathrm{cs}}$

[TABLE]

This is a characteristic time for a spherical particle placed in a bulk fluid needed to come to equilibrium, whose equation of motion is governed by Newton’s law with a Stokes drag, i.e. $mdv/dt=-6\pi\eta Rv$. As it is clear from Eq. (26), the relaxation time changes either upon variation of the bead radius or the fluid viscosity, hence the properties of the fluid. Indeed, scaling frequencies using $1/\tau_{\mathrm{cs}}$ allows us to collapse all the peaks in Fig. 9a, suggesting that the optimal swimming speed in our LB simulations arises from integrating coasting in most appropriate manner. Consequently, the average maximum speeds are reached for moderate $L_{ij}$ ($\approx 3.2\times 2R$) and are of the order of $\approx 0.001\cdot 2R$ per period of the oscillating field as seen in Fig. 9a. This result is remarkable because it shows the emergence of inertial effects which cannot be captured by a simple Stokes dynamics.

The Stokes model (Appendix A) is, however, capable of reproducing several central results of LB-simulations:

1. The swimmer propagates perpendicularly to the excitation (Figs. 9b-d, and 9e);
2. The average speed of swimming decreases with increasing inter-particle distance in $L>3\times 2R$, which is in agreement with the far-field description of hydrodynamics in the model. For small $L_{ij}$, the LB simulations behave differently due to short-range capillary interactions between the spheres which cannot be described by the bead-spring model. Notably, higher order treatments of the hydrodynamic interactions do not change the trend in the frequency response.
3. The peak in the velocity-frequency response is sensitive to the inter-particle distance.
4. The characteristic frequency grows with increasing $L$. As for the behaviour of the average speed of swimming, the LB-simulations are expected to behave differently for small $L$.

Besides the dependency of speed dependence on the frequency (Fig. 9a), several additional effects cannot be captured by our Stokes modes. Namely, at small $L_{ij}$ ($\approx 2.3\times 2R$) the swimmer moves efficiently for a wide range of frequencies and the main peak is broadened by the relative proximity of the particles.

The motion at low frequencies is, furthermore, dominated by the sizable deformations of the interface and the capillary potential has a much more complex structure than described by Eqs. (17) and (20), making the determination of the spring constant even more challenging, given the current setup of the simulation. Importantly, the deformations of the interface introduce an additional time scale $\tau_{\mathrm{int}}$

[TABLE]

Interestingly, using this time scale equally well collapses the peaks shown in Fig. 9a, suggesting that this could be the relevant process dominating the dynamics. However, our test simulations with a significantly larger box do not reproduce this result, while the coasting time scale reappears as the relevant one in the system.

V Summary

In this study we performed extensive lattice Boltzmann simulations of magnetocapillary swimmers consisting of three ferromagnetic beads trapped at a fluid/fluid interface. The simulations are inspired by the experiments reported in Refs. Vandewalle et al. (2012); Lumay et al. (2013); Grosjean et al. (2015), where ferromagnetic beads placed at the water/air interface were driven by magnetic fields and reached propulsion speeds of hundreds of $\mu$m/s.

To match the corresponding regimes for the capillary potentials of the experiments we studied equilibrium properties of a single (Fig. 3), two (Fig. 5) and three particles (Fig. 6) at the interface. The resulting dependence of the center-to-center distance for two and three particles on the strength of the magnetic moment unveiled three main regimes that are in agreement with the experimental observations (Figs. 3 and 6 of Ref. Vandewalle et al. (2012)).

Adding a smaller oscillating magnetic field drives the system out of the local equilibrium and results in a directed motion of the three-bead swimmer (Fig. 7). By fixing the direction of the oscillating field we demonstrated (Fig. 9a) that the average speed of the swimmer nonmonotonously depends on the frequency of the $\bm{B}(t)$-field and reaches its maximum at approximately $0.001$ $2R/\tau^{\mathrm{LB}}_{\mathrm{v}}$ for frequencies in the vicinity of the inverse coasting time defined by Eq. (26). The emergence of this time scale is interesting as the swimmer operates in near the Stokes regime, yet the inertial effects start to appear. Consequently the time scale of the swimmer characteristic for $Re=0$Pande and Smith (2015) is not determined as the relevant one, although certain elements of the Stoksian dynamics are recovered. Furthermore, the motion observed in LB simulations is significantly different to experiments with magnetocapillary swimmers (Refs. Vandewalle et al. (2012); Lumay et al. (2013); Grosjean et al. (2015)), since in these systems the motion is dominated by inertia and a classical resonance is observed, at least in the linear geometry. The understanding of the emergent time-scale $\tau_{cs}$ requires thus further theoretical investigations which we hope to undertake in future.

The LB-simulations performed here also point out to the importance of simultaneous optimization of the swimmer’s geometry and driving which is a very important property for applications: swimmers with different inter-particle distances ($L_{ij}$) might be controlled independently by application of AC-fields with distinct frequencies. In the context of the current work, however, besides identifying these interesting phenomena, we clearly establish lattice Boltzmann simulations as a technique for studying microswimming on liquid interfaces, which hitherto was not possible.

V.1 Appendix: Triangular bead-spring swimmer in the Stokes regime

We assume three beads of equal radii $R$ pairwise connected by harmonic springs of spring constant $k$ and equilibrium length $L$ (Fig. 10), such that the potential energy of a spring connecting two beads at positions $\bm{r}_{i}$ and $\bm{r}_{j}$ reads

[TABLE]

In addition to the spring constant $k$, external forces denoted by $\bm{F}_{i}(t)$ act on each bead $i$, satisfying the condition that (a) the sum over all forces $\bm{F}{{}_{\mathrm{tot}}}(t)=\sum_{i=1}^{3}\bm{F}_{i}(t)=0$ as well as (b) the resulting torque $\bm{T}{{}_{\mathrm{tot}}}(t)=\sum_{i=1}^{3}\bm{r}_{i}\times\bm{F}_{i}(t)=0$ both vanish at each time step. Note that the origin of the coordinate system relative to which the torque is calculated is irrelevant when (a) holds. For the calculations below, we specifically set the swimmer in the $XZ$-plane, and exert a sinusoidal force specified by $\bm{F}_{1}(t)=A\sin(\omega t)\bm{e}{{}_{\mathrm{x}}}$ on the first bead. Consequently, forces on the other two beads are set by the constraints (a) and (b). Notably, this protocol allows only for pure translations of the swimmer’s center of mass in the steady state.

In the regime of the bead radius $R\ll L$, the hydrodynamic interaction between the beads is given by the Oseen tensor Oseen (1910)

[TABLE]

and the equation of motion for each bead is governed by the Stokes law, such that the system is described by the set of equations of motion describing the time evolution of the position $\bm{r}_{i}$ of each bead

[TABLE]

where $\nu_{j}=\frac{1}{6\pi\eta R}$ and $\bm{G}$ is the spring force $\bm{G}(\bm{r}):=-\bm{\nabla}\phi(\bm{r})$. The first (mobility) term of Eq. (30) describes the motion due to external and spring forces directly acting on the beads, while the second accounts for the hydrodynamic interactions between the beads.

The swimming velocity emerges from inserting the solution of Eq. (30) into Eq. (22). With the above-prescribed driving, the swimmer’s motion is found always to be associated with a translation in the $Z$-direction without a net rotation.

Acknowledgements.

This work was financially supported by the DFG Priority Programme SPP 1726 “Microswimmers – From Single Particle Motion to Collective Behaviour” and the Cluster of Excellence “Engineering of Advanced Materials” (project EXC 315). We further acknowledge the Jülich Supercomputing Centre (JSC) and the High Performance Computing Center Stuttgart (HLRS) for the allocation of computing time.