Controlling stability and transport of magnetic microswimmers by an external field

As a microswimmer propels itself through a fluid, it generates a long-range disturbance. This perturbation propagates through the fluid and influences the motion of other swimmers. Self-propulsion in conjunction with fluid-mediated interactions in active suspensions give rise to a wealth of collective phenomena that are very distinct from those found in passive systems at equilibrium [1–5]. Some examples include hydrodynamic instabilities that lead to spatio-temporal pattern formation [6–8], active turbulence [9–12] and unusual rheological properties [13–16]. Moreover, microswimmers exhibit new motility patterns in response to external stimuli such as chemical signals [12,17,18], light [19,20], gravitational [21–25] and magnetic fields [26–30]. The control and regulation of collective motion of microswimmers via an external field offers a promising route for their exploitation in high-tech applications such as micro-scale cargo transport, targeted drug delivery, and microfluidic devices [31–34].

Presently, a theoretical understanding of collective behaviour and transport of microswimmers in an external field is largely missing. Here, we put forward a kinetic theory for active suspensions that extends the previous kinetic models [6,7] to include the effects of an external field. Our model is applicable to any active suspension driven by an external aligning torque. Examples include magnetotactic bacteria (MTB) carrying an intrinsic weak magnetic moment [35–39] and weakly magnetised artificial swimmers [40–49] in an external magnetic field or bottom-heavy swimmers in a gravitational field [22]. MTB driven by a sufficiently strong magnetic field exhibit particularly intriguing patterns of collective behavior such as band formation [26,27] and pearling instability under flow [28]. Thus, we focus on the dynamics of active magnetic suspensions in a uniform magnetic field.

We investigate the instabilities and transport of dilute suspensions of spherical magnetic swimmers in an external field combining linear stability analysis and 3D numerical simulations. At low magnetic fields, a homogeneous weakly polarised state is stable, akin to an isotropic suspension of spherical swimmers. However, for sufficiently high activity strengths and moderately strong magnetic fields, a homogeneous polar phase becomes unstable for both pushers and pullers. These instabilities significantly reduce the polarisation of the swimmers and lead to a decrease in their mean transport speed. In the unstable regime, we observe a rich phenomenology of pattern formation by varying the magnetic field and activity strengths. Representative examples of pattern development for pushers and pullers are shown in fig. 1. Notably, pushers and pullers exhibit distinct instability patterns that result from the divergence of bend and splay fluctuations [50–52], respectively. The 3D visualisation of the density field in fig. 1(c) and (d) (see also the corresponding supplementary videos pusher.mp4 and puller.mp4) clearly shows that the pushers concentrate in band-like structures perpendicular to the magnetic field that migrate in the field direction, whereas pullers form lane-like patterns parallel to the field. Our results for the pushers are remarkably similar to the observed magnetotactic bands reported for spherical MTB [26,27]. The instability of the polar state induced by an external field shares similarities with the instability of aligned swimmers with nematic or polar interactions [6,7]. However, for an externally induced polar state such an instability disappears by a further increase of the magnetic field strength. To our knowledge such a re-entrant hydrodynamic stability has not been previously reported in active systems and calls for further experimental investigations.

Zoom In Zoom Out Reset image size

We consider a dilute suspension of N spherical magnetic microswimmers with a hydrodynamic radius a immersed in a fluid of volume V at a number density $\varrho=\frac{N}{V}$. We assume that the self-propulsion is generated by a force-free mechanism of hydrodynamic origin such that the far field flow of a swimmer is well represented by that of a point-force dipole with an effective dipolar strength S_eff [7,53–55]. S_eff depends on the geometrical parameters of the model swimmer [42,43,45,55], for instance on a and the flagellum length ℓ [55]. Each swimmer carries a weak magnetic dipole moment configuration of an ensemble of the $\mu_m$ along its body axis $\mathbf{\hat{n}}$ and has a self-propulsion velocity $v_{\mathrm{sp}}\mathbf{\hat{n}}$. The suspension is exposed to a uniform magnetic field B that exerts an aligning torque on each swimmer. The magnetic moment values of MTB are of the order of $\mu_m \approx 1 \times10^{-16}\text{JT}^{-1}$ [36,38,56,57] and their size $a \approx 1\ \mu \text{m}$. For dilute suspensions with inter-particle distances $r \gtrsim 3 a$, their dipole-dipole interactions are small compared to the thermal energy scale and we can neglect their effect. Instead, we focus on the interplay between the hydrodynamic interactions and the aligning torque.

For sufficiently low ϱ, the mean-field configuration of an ensemble of the swimmers at a time t can be described by the probability density $\frac{1}{V} \Psi (\textbf{x}, \mathbf{\hat{n}}, t)$ of finding a particle with the center-of-mass position x and the unit orientation $\mathbf{\hat{n}}$. $\Psi$ is normalized such that $\frac{1}{V} \int \text{d}\mathbf{x} \int \text{d}\mathbf{\hat{n}} \Psi = 1$. The time evolution of Ψ is governed by a Smoluchowski-equation of the form

$$\begin{equation} \partial_{t} \Psi + \nabla \cdot \mathbf{J}_{\mathrm{tr}}[\Psi] + \nabla_{\mathrm{S}} \cdot \mathbf{J}_{\mathrm{rot}}[\Psi] =\mathbb{D} \Psi, \end{equation} \tag{ 1 }$$

where $\nabla_{\mathrm{S}} = (\mathds{1} - \mathbf{\hat{n}} \mathbf{\hat{n}}) \cdot \nabla_{\mathbf{\hat{n}}}$ denotes the angular gradient operator; J_tr and J_rot are the translational and rotational drift currents. $\mathbb{D} = D_t \nabla^2 + D_r \nabla_{\text{S}}^{2}$, with $\nabla_{\text{S}}^{2}$ as the Laplace operator on a unit sphere, accounts for the evolution of Ψ resulting from the translational and rotational diffusive currents. D_t and D_r represent the effective long-time translational and rotational diffusion coefficients that can be of thermal or biological origin, e.g., due to tumbling of bacteria. $\mathbf{J}_{\mathrm{tr}} \equiv \Psi \mathbf{v_{\mathrm{tr}}}$ describes the translational current stemming from the self-propulsion of a swimmer and its convection in the local flow u,

$$\begin{equation} \mathbf{v}_{\mathrm{tr}} = v_{\mathrm{sp}} \mathbf{\hat{n}} + \mathbf{u}[\Psi]. \end{equation} \tag{ 2 }$$

The rotational current $\mathbf{J}_{\mathrm{rot}} \equiv \Psi \mathbf{v_{\mathrm{rot}}}$ incorporates contributions from the rotational velocities resulting from the torque generated by the aligning magnetic field and the local flow vorticity $W_{jl} \equiv 1/2 ( \partial_j u_l - \partial_l u_j )$ according to Jeffery's equation [58,59]:

$$\begin{equation} \mathbf{v}_{\mathrm{rot}} = (\mathds{1} - \mathbf{\hat{n}} \mathbf{\hat{n}}) \cdot \left( \mu_m / \xi_{\mathrm{R}} \mathbf{B} - \mathbf{W}[\mathbf{u}] \cdot \mathbf{\hat{n}} \right), \end{equation} \tag{ 3 }$$

where $\xi_{\mathrm{R}}$ is the rotational friction coefficient.

The flow field $\mathbf{u}[\Psi]$ in the low Reynolds number limit is captured by the incompressible Stokes equation

$$\begin{equation} \eta \Delta \mathbf{u} - \nabla P = - \nabla \cdot \mathbf{\Sigma}[\Psi], \quad \nabla \cdot \mathbf{u} = 0 \end{equation} \tag{ 4 }$$

in which P denotes the isotropic pressure and η the viscosity of the suspending fluid and $\nabla \cdot \mathbf{\Sigma} = \partial_{x_i} \Sigma_{ij} \hat{\mathbf{e}}_j$. The flow is determined by the state of the system encoded by Ψ via a mean-field stress profile $\mathbf{\Sigma}[\Psi]$. It includes two contributions: an active stress $\mathbf{\Sigma}_a$, generated by the self-propulsion of force-free dipolar swimmers [53,54], and an antisymmetric magnetic stress $\mathbf{\Sigma}_m$, caused by reorientation of swimmers in the magnetic field. The active stress is proportional to the angular expectation value of the nematic order tensor $\mathbf{\Sigma}_a(\mathbf{x}) = \Sigma_a \int \text{d}\mathbf{\hat{n}} \Psi \, (\mathbf{\hat{n}} \mathbf{\hat{n}} - \frac{1}{3} \mathds{1})$ [7,60]. The strength of the active stress is given by $\Sigma_a = \pm \varrho S_{\mathrm{eff}}$. The sign of $\Sigma_a$ determines the nature of the swimmers, being a puller $\Sigma_a > 0$ or a pusher $\Sigma_a < 0$. The magnetic stress is given by $\mathbf{\Sigma}_m(\mathbf{x}) = \frac{\Sigma_m}{2} \int \text{d}\mathbf{\hat{n}} \Psi \, ( \mathbf{\hat{n}} \mathbf{\hat{B}} - \mathbf{\hat{B}} \mathbf{\hat{n}}),$ in which $\mathbf{\hat{B}} = \mathbf{B} / B$ and $\Sigma_m = \varrho \mu_m B$ [61]. Note that the symmetric part of the stress is zero for spherical particles [61].

To facilitate the analysis of our model, we render the equations dimensionless, using the following characteristic velocity, length, and time scales: $u_c=v_{\mathrm{sp}}$, $x_c=\varrho^{-1/3}$ (average inter-particle distance) and $t_c=x_c/u_c$. These scaling choices leave the distribution function unchanged: $\Psi(\mathbf{x}, \mathbf{\hat{n}}, t) \equiv \Psi_{\mathrm{scaled}}(\mathbf{x}/x_c, \mathbf{\hat{n}}, t/t_c)$. The corresponding dimensionless model parameters are the magnetic field strength $b = t_c \mu_m B / \xi_{\mathrm{R}}$, the rotational and translational diffusion coefficients $d_r = D_r \varrho^{-1/3} v_{\mathrm{sp}}^{-1}$ and $d_t = D_t \varrho^{1/3} v_{\mathrm{sp}}^{-1}$, the active stress amplitude $\sigma_a = t_c \eta^{-1} \Sigma_a = \varrho^{2/3} S_{\mathrm{eff}} v_{\mathrm{sp}}^{-1} \eta^{-1}$ and the magnetic stress amplitude $\sigma_m = t_c \eta^{-1} \Sigma_m = \varrho^{2/3} \mu_m B v_{\mathrm{sp}}^{-1} \eta^{-1}$.

Let us first consider a solution of eq. (1) satisfying, $\partial_t \Psi_0 = 0$ and $\nabla \Psi_0 = 0$. $\Psi_0$ is given by

$$\begin{equation} \Psi_0(\mathbf{\hat{n}}) = \frac{\alpha}{4 \pi \sinh \alpha} e^{\alpha \mathbf{\hat{n}} \cdot \mathbf{\hat{B}}}. \end{equation} \tag{ 5 }$$

in which $\alpha = b/d_r \equiv \mu_m B/ (\xi_{\mathrm{R}} \, D_r)$ and it is identical to the steady state solutions obtained in [16,28,39]. We call α the alignment parameter as it is equal to the ratio of two characteristic reorientation times; $\alpha \equiv \tau_r / \tau_m$. $\tau_r = D_r^{-1}$ represents the average decorrelation time of the particle from its initial orientation and $\tau_m = \xi_{\mathrm{R}} / \mu_m B$ describes the typical time a non-diffusive particle needs to align itself with the magnetic field. Hence, the degree of alignment is determined by the competition between the aligning magnetic torque and the randomizing rotational diffusion. The $\Psi_0(\mathbf{\hat{n}})$ corresponds to a homogeneous polar state with a mean polarization vector $\mathbf{p}_0 \equiv \int \text{d}\mathbf{\hat{n}}\, \Psi_0(\mathbf{\hat{n}})=p_0(\alpha) \mathbf{\hat{B}}$ in which

$$\begin{equation} p_0(\alpha) = \coth \alpha - 1 / \alpha, \end{equation} \tag{ 6 }$$

is known as the Langevin function in the context of paramagnetism. Note that a full alignment is only achieved for $\alpha \gg 1$.

We investigate the linear stability of the homogeneous polar state by considering a small perturbation of the form $\Psi_0 + \varepsilon \, \tilde{\Psi}(\mathbf{k}, \mathbf{\hat{n}}) e^{{i} \mathbf{k} \cdot \mathbf{x} + \lambda t}$. The equation of motion linearised in $0 < \varepsilon \ll 1$ can be expressed as an eigenvalue problem of the form $\mathbb{L} \tilde{\Psi} = \lambda \tilde{\Psi}$, where $\mathbb{L}(\mathbf{k},\mathbf{\hat{n}},\mathbf{b},\Psi_0)$ is a linear differentio-integro-operator (see the Supplementary Material Supplementarymaterial.pdf (SM)). We solve the eigenvalue problem in the basis of spherical harmonics $Y_l^h$ numerically by truncating the matrix $\langle Y_l^h | \mathds{L} | Y_j^m\rangle$ at sufficiently large number of modes such that the convergence of the dominant eigenvalues are ensured. The external field breaks the rotational symmetry. Hence, the stability depends on the direction $\hat{\mathbf{k}}$ of the perturbation wave vector with respect to the magnetic field direction, which can be characterized by a single angle $\Theta_{\mathrm{B}}\equiv \cos^{-1}(\hat{\mathbf{k}} \cdot \hat{\mathbf{B}})$.

We first examine the stability of swimmers with moderate values of activity and magnetic field strength leading to $\sigma_a = \pm 3/2$ and $\alpha=4$. The remaining parameters are chosen to be comparable to those of MTB [28] and they are fixed to: $\sigma_m=0.01 \, \alpha \, d_r$, $d_r = 0.05$, and $d_t = 0.003$. Figure 2 shows the real part of the eigenvalue with the largest magnitude $\text{Re} \lambda_{\max}(k)$, the so-called maximum growth rate, as a function of $k=|\mathbf{k}|$ at various perturbation angles $\Theta_{\mathrm{B}}$. For both puller and pusher swimmers, long-wavelength perturbations dominate the instabilities and destabilize the homogeneous polar state defined by p₀. For pushers, fluctuations of the linearised equations in the direction of magnetic field grow fastest whereas for pullers both perturbation directions parallel and perpendicular to B predominate. Thus, we expect pushers and pullers to exhibit distinct instability patterns as confirmed by the non-linear dynamics simulations; see fig. 1.

Zoom In Zoom Out Reset image size

Next, we present the stability phase diagram in which we vary the strengths of both activity $S_{\mathrm{eff}} \propto \sigma_a$ and magnetic field $B \propto \alpha$. Figure 3 depicts the stability diagram for $\Psi_0(\alpha)$ in the $(\sigma_a,\,\alpha)$-plane. The magnetic stress is varied concomitant with α as $\sigma_m = 0.01\, \alpha d_r$; the remaining parameters are kept constant at the values given in the caption. With this choice of parameters the diagram should represent an experimentally accessible range. From the linear stability analysis, we determine the border lines that separate the stable from the unstable regions. Note that the steady state $\Psi_0(\alpha)$ and its polarization depend on α. $p_0(\alpha)$ given by eq. (6) is shown in the right panel of fig. 3. At low values of α where polarization is weak, $p_0(\alpha) \lesssim 0.3$, $\Psi_0(\alpha)$ remains stable. At moderately strong fields and for $\sigma_a \gtrsim 1$, the hydrodynamic interactions become amplified as a result of the increased polarization and destabilize the steady state. Strikingly, at stronger magnetic fields the hydrodynamic instabilities can be overcome. At such strong fields, the randomizing effect of the rotational diffusion can be ignored. The magnetic torque $\propto \alpha$ dominates over the hydrodynamic torque and the steady state becomes stable again. This is a consequence of the fact that for a fixed active stress amplitude the hydrodynamic stress and the resulting torque have an upper bound (perfect polarization), whereas we can increase the external torque by increasing the magnetic field. To examine the validity of these predictions, we study the dynamics of swimmers by non-linear simulations.

Zoom In Zoom Out Reset image size

We perform non-linear simulations of the kinetic model in three dimensions to study the long-time dynamics and pattern formation resulting from the instabilities. To solve the Smoluchowski equation, eq. (1), with periodic boundary conditions, we use a hybrid stochastic particle based sampling method to obtain $\Psi(\mathbf{x},\mathbf{\hat{n}},t)$ and a spectral method to solve for the flow field $\mathbf{u}(\mathbf{x})$. For further details see the SM. Probing the stability of $\Psi_0(\alpha)$ for different activity and magnetic field strengths by simulations, we find excellent agreement with the predictions of the linear stability analysis, as demonstrated in fig. 3. $\Psi_0(\alpha)$ is stable for the parameter points denoted by discs, whereas $(\sigma_a,\,\alpha)$ values for which Ψ evolves towards an inhomogeneous time-dependent density profile are depicted by stars. The angular distribution of these unstable states appears to converge towards a steady state, whereas their density fields evolve towards dynamic spatial patterns. In the unstable regime, density and polarization gradients generate a flow with an inhomogeneous vorticity field that is coupled to the swimmers orientations and rotates them away from the $\mathbf{\hat{B}}$ direction. To quantify these disorienting effects, we measure the time-averaged global polarization defined as $\bar{\mathbf{p}}_t = \frac{1}{N_t} \sum_{j=0}^{N_t} \|\langle\mathbf{\hat{n}} \rangle(t_0 + j \Delta t)\|$, where $\Delta t = 40 \, \delta t$ (with the simulation time step $\delta t = 0.2 \,$ in units of t_c) is about an order of magnitude larger than the reorientation time $\tau_m$. The braces ${\langle \bullet \rangle} \equiv \frac{1}{V} \int \text{d}\mathbf{x} \int \text{d}\mathbf{\hat{n}}\, \Psi(\mathbf{x}, \mathbf{\hat{n}}, t) \bullet$ define the expectation value. t₀ marks a relaxation time after which $\bar{p}_t$ is nearly time-independent despite exhibiting non-stationary patterns (see fig. 3 in the SM). Choosing $N_t=20$ allowed us to obtain sufficient statistics. We find that $\mathbf{\hat{B}} \cdot \mathbf{\bar{p}}_t \approx \bar{p}_t$ and $\bar{p}_t$ is independent of the system size for $L\gtrsim 50$ (see SM). Figures 4(c) and (d) present the $\bar{p}_t$ as a function of α for pushers and pullers at different activity strengths $\sigma_a$. For moderate $\sigma_a$ and α values falling in the unstable regime, we observe a significant reduction in the polarization compared to $p_0(\alpha)$ (eq. (6)), shown by the dashed line. The decrease in the mean polarization is stronger for larger activity strengths. Stronger magnetic fields drive the system back into the stable regime and $\bar{p}_t$ agrees with $p_0(\alpha)$ in those regions.

Zoom In Zoom Out Reset image size

The mean polarization governs the mean transport speed $\bar{v}$ in the direction of magnetic field that additionally includes a contribution from the convective flow component along $\mathbf{\hat{B}}$:

$$\begin{equation} \bar{v} = \mathbf{\hat{B}} \cdot \left< v_{\mathrm{sp}} \mathbf{\hat{n}} + \mathbf{u} \right> = (v_{\mathrm{sp}} \bar{\mathbf{p}}_t + \bar{\mathbf{u}}_t) \cdot \mathbf{\hat{B}}. \end{equation} \tag{ 7 }$$

For an efficient transport in the direction of the magnetic field, a high polarization of swimmers parallel to B is desirable that can be achieved by increasing the field strength. To evaluate the contribution of $\bar{\mathbf{u}}_t \cdot \mathbf{\hat{B}}$ to the transport, we calculate the space and time-averaged flow velocity as $\mathbf{\bar{u}}_t = \frac{1}{N_t} \sum_{j=0}^{N_t} \langle \mathbf{u} \rangle(t_0 + j \Delta t)$.

Figures 4(e) and (f) show $\mathbf{\bar{u}}_t \cdot \mathbf{\hat{B}}$ vs. $\alpha \propto B$ for pushers and pullers at different values of $\sigma_a$ that is almost independent of box size for $L\ge 50$ (see SM). The mean flow velocity created by pushers has a vanishing component along $\mathbf{\hat{B}}$. Hence, their transport speed is determined by the mean polarization whereas for the pullers the contribution of convective flow to the transport is significant. This dissimilarity originates from distinct nature of instabilities that prevail the pushers and pullers and distort their polarization field $\mathbf{p}(\mathbf{x},t) = \frac{1}{V \rho(\mathbf{x},t)} \int \text{d}\hat{\mathbf{n}} \; \Psi(\mathbf{x},\hat{\mathbf{n}},t) \, \hat{\mathbf{n}}$, in which $\rho(\mathbf{x},t) = \frac{1}{V} \int \text{d}\hat{\mathbf{n}} \; \Psi(\mathbf{x},\hat{\mathbf{n}},t)$. For pushers, bending deformations in their polarization field produce alternating shear flow layers perpendicular to $\hat{\mathbf{B}}$ as schematically drawn in fig. 4(a). Because of imbalance of the magnetic and the flow-induced torques, the bending fluctuations are further amplified and decrease the mean polarization. These distortions also increase the density where $\nabla \cdot \mathbf{p} < 0$. As a result, pushers form dense layers perpendicular to $\hat{\mathbf{B}}$ that migrate parallel to $\hat{\mathbf{B}}$, see fig. 1. Similarly for pullers, splay deformations generate alternating pillar-like flow regions along $\hat{\mathbf{B}}$ as shown in fig. 4(b). Denser regions of pullers where $\nabla \cdot \mathbf{p} < 0$ coincide with regions carrying a flow anti-parallel to $\hat{\mathbf{B}}$. They result in a net convection anti-parallel the magnetic field (fig. 4(f)) and reduce the mean transport speed.

The mean polarization of magnetic swimmers increases continuously with the field at a pace that depends on their activity strength. For $B \lesssim 0.1\ \text{mT}$, a weakly polar state is always stable. At moderate fields, for sufficiently strong activities where $\bar{p} \gtrsim 0.3$, pushers and pullers exhibit distinct bend and splay instability patterns and propose a pragmatic approach for distinguishing them in experiments. These instabilities highlight the significance of hydrodynamic interactions that hinder the directed transport of swimmers. The reduced transport speed results from the coupling of density and flow field to the distortions of the polarization field. Interestingly, at field strengths beyond an activity-dependent value ($B \gtrsim 1\ \text{mT}$ for $\Sigma_a \approx 0.02\ \text{Jm}^{-3}$ and $B \gtrsim 3 \ \text{mT}$ for $\Sigma_a \approx 0.04\ \text{Jm}^{-3}$) a homogeneous polar state becomes stable. We defer a classification of patterns as a function of activity and magnetic field strengths to a future work.

Finally, we note that our results are valid in the limit of negligible magnetic interactions that are of relevance to dilute suspensions of swimmers with weak magnetic dipole moments such as magnetotactic bacteria. For synthetic magnetic microswimmers with larger magnetic dipole moments or dense suspensions, the magnetic dipolar interactions alone can lead to clustering instabilities [62] and the interplay between long-range magnetic and hydrodynamic interactions on development of instabilities deserves to be explored. Moreover, clarifying the role of swimmer-swimmer correlations [63], and near-field hydrodynamic interactions in more concentrated suspensions merits further investigations.

We thank Eric Clément, M. Cristina Marchetti and Friederike Schmid for fruitful discussions and Tapan Adhyapak for a critical reading of the manuscript. We acknowledge the financial support from the German Research Foundation (http://www.dfg.de) within SFB TRR 146 (https://trr146.de). The simulations were performed using the MOGON II computing cluster. This research was supported in part by the National Science Foundation under Grant No. NSF PHY17-48958.