Tuning a magnetic field to generate spinning ferrofluid droplets with controllable speed via nonlinear periodic interfacial waves

Recently, there has been significant interest in the physics of active and responsive fluids [1, 2]. For example, swimming bacteria can take a suspension of microscopic gears out of equilibrium and extract rectified (useful) work out of an otherwise random system [3]. One promising approach to creating active fluids with controllable properties and behaviors is by suspending many mechanical microswimmers made from shape-programmable materials [4] and actuating them with an external magnetic field [5, 6]. This actuation mechanism is particularly enticing for biological applications due to the safe operation of magnetic fields in the medical setting (for, e.g., targeted therapies and drug delivery in vivo) [7]. Even simpler than a suspension of magnetically-responsive mechanical microswimmers is a suspension of ferrofluid droplets, which can also respond to an external magnetic field [8, 9]. Ferrofluids are colloidal dispersion of ferromagnetic nanoparticles in a carrier liquid, such as water, which can be immiscible when placed in another liquid. However, the ferrofluid droplet’s interface motion and response to different types of external magnetic fields is not well understood. Previous work has addressed the linear stability of such fluid–fluid interfaces [10, 11], including stationary shapes [12], but not a droplet’s nonlinear dynamics or controllable motion. Guided by the well-established ability of nonlinearity to “arrest” long-wave instabilities [13], we demonstrate, using theory and nonlinear simulation, that it is possible to “grow” linearly unstable ferrofluid interfaces into well-defined permanent shapes. These permanent shapes, which cannot be further deformed without changing the forcing of the system, can then be considers as solitary waves, in the sense of a “localized wave that propagates along one space direction only, with undeformed shape” [14, p. 11]. Importantly, unlike previous work discussing traveling waves on a ferrofluid interface in a Cartesian configuration [15], we analyze the fully nonlinear dynamics of these waves in a novel configuration, and thus ensure they satisfy the solitary wave definition above. Specifically, we show that the resulting coherent droplet shapes are reproducible and controllable via an external magnetic field. These droplets can be set into rotational motion with velocities predictable by the proposed theory, leading to the possibility of an externally-actuated active fluid suspension.

We study the dynamics of an initially circular ferrofluid droplet (radius $R$) confined in a Hele-Shaw cell with gap thickness $b$ and surrounded by air (negligible viscosity), as shown in Fig. 1, because “[i]f any [ferro]fluid mechanics problem is likely to be accessible to theory and to direct comparison of theory and experiment it should be this one” [16]. Both fluids are considered incompressible. We propose to apply the radially-varying external magnetic field

A long wire through the origin, carrying an electric current $I$, produces the azimuthal component $\bm{\mathrm{H}}_{a}$. Anti-Helmholtz coils produce the radial component $\bm{\mathrm{H}}_{r}$, where $H_{0}$ is a constant and $L$ is a length scale [17, 12]. The combined magnetic field $\bm{\mathrm{H}}=\bm{\mathrm{H}}_{a}+\bm{\mathrm{H}}_{r}$ forms an angle with the initially undisturbed interface [18]. The droplet experiences a body force $\propto|\bm{\mathrm{M}}|\bm{\mathrm{\nabla}}|\bm{\mathrm{H}}|$, where $\bm{\mathrm{M}}$ is the magnetization. To study shape dynamics, we assume the ferrofluid is uniformly magnetized, $\bm{\mathrm{M}}=\chi\bm{\mathrm{H}}$, where $\chi$ is its constant magnetic susceptibility. So, $\bm{\mathrm{\nabla}}|\bm{\mathrm{H}}|\neq\bm{\mathrm{0}}$ is the main contribution to the body force, and the demagnetizing field is negligible, as shown in previous work [18, 19, 17, 20].

Enforcing no-slip on the confining boundaries and neglecting inertial terms, the confined flow is governed by a modified Darcy’s law [17] with gap-averaged velocity:

where $p$ is the pressure in the droplet, $\eta$ is the ferrofluid’s viscosity, $\Psi=\mu_{0}\chi|\bm{\mathrm{H}}|^{2}/2$ is a scalar potential accounting for the magnetic body force, and $\mu_{0}$ is the free-space permeability. Here, $\bm{\mathrm{v}}$ is the velocity field of the “inner” ferrofluid, while the viscosity of the “outer” fluid is considered negligible (i.e., it is considered inviscid), so the flow exterior to the droplet is neglected. The resulting model is thus, essentially, a one-phase model.

At the boundary of the droplet, the pressure is given by a modified Young–Laplace law [8, 9]:

where $\tau$ is the constant surface tension, and $\kappa$ is the curvature of the droplet shape. The second term on the right-hand side of Eq. (3) is the magnetic normal traction [8, 9], where $\bm{\mathrm{\hat{n}}}$ denotes the outward unit normal vector at the interface. This contribution breaks the symmetry of the initial droplet interface, due to the projection of $\bm{\mathrm{M}}$ onto $\bm{\mathrm{\hat{n}}}$, and causes the droplet to rotate. The kinematic boundary condition

requires that the droplet boundary is a material surface.

We employ the weakly nonlinear approach [21] previously adapted to ferrofluid interfacial dynamics (e.g., [18, 17, 12]). The droplet interface is written as $h(\theta,t)=R+\xi(\theta,t)$, where

represents the perturbation of the initially circular interface, with complex Fourier amplitudes $\xi_{k}(t)\in\mathbb{C}$ and azimuthal wavenumbers $k\in\mathbb{Z}$. The velocity potential $\phi=p-\Psi$ is then expanded into a Fourier series as

and $\phi_{k}$ is expressed in terms of $\xi_{k}$ through the kinematic boundary condition (4). Substituting Eqs. (5), (6) and (3) into Eq. (2), keeping only terms up to second order in $\xi$, we find the dimensionless equations of motion ($k\neq 0$):

The mode-coupling functions in Eq. (7) are given by

where $\operatorname{sgn}(x)=x/|x|$ for $x\neq 0$ and $\operatorname{sgn}(0)=0$.

From mass conservation, $\xi_{0}=-\sum_{k>0}|\xi_{k}|^{2}/R$ $\forall t\geq 0$. Here,

denotes the (complex) linear growth rate, and

are the magnetic Bond numbers quantifying the ratio of azimuthal and radial magnetic forces to the capillary force, respectively. Terms multiplied by $\chi$ arise from the magnetic normal stress. The time and length scales used in the nondimensionalization are $12\eta L^{3}/\tau b^{2}$ and $L$, respectively.

First, consider Eq. (7), neglecting quadratic terms in $\xi_{k}$, then $\Re[\Lambda(k)]=\lambda(k)$ governs the exponential growth or decay of infinitesimal perturbations. For $\lambda(k)>0$, the interface is unstable. Specifically, Eq. (9) indicates that the radial magnetic field term $\propto(1+\chi)\mathrm{N_{Br}}$ is destabilizing, while the azimuthal term $\propto\mathrm{N_{Ba}}$ and surface tension are stabilizing. The most unstable mode $k_{m}$ solves $d\lambda(k)/dk=0$:

This wavenumber characterizes the dominant $\lfloor k_{m}\rfloor$-fold symmetry of a pattern. Note that the normal stress from the azimuthal magnetic field does not contribute to the linear dynamics.

The phase velocity of each mode,

in the linear regime, is set by $\Im[\Lambda(k)]$. A periodic shape on $[0,2\pi]$ forms a closed curve, meaning wave propagation is manifested as rotation of the droplet. Motion is caused by the magnetic normal stresses arising from the combined magnetic field. Intuitively, from vector projection, we observe that only the combined azimuthal and radial magnetic field can break the symmetry and cause a force imbalance leading to motion. This linear analysis indicates that perturbations of the droplet interface can propagate (and, since $v_{p}=v_{p}(k)$, they also experience dispersion). Such wavepackets will either decay or blow-up exponentially according to the sign of $\lambda(k)$. However, this is not the whole story, and nonlinearly stable traveling shapes exist, as we now show.

To demonstrate the possibility of nonlinear traveling waves in this system, we numerically solve the weakly nonlinear mode-coupling equations (7) for five modes (i.e., $k,2k,\ldots,5k$). The fundamental mode $k=7$ is chosen to allow propagating solutions over a wider swath of the $(\mathrm{N_{Ba}},\mathrm{N_{Br}},k_{m})$ space (compared to choosing $k<7$), while only requiring modest spatial resolution for simulations (compared to $k>7$). We verified that the amplitudes $|c_{n}|$ and phases $\angle{[c_{n}]}$ of modes saturate at late times, leading to permanent propagating profiles with $\xi_{nk}(t)=c_{n}e^{in\omega(k)t}$ (see below).

Next, we perform fully nonlinear simulations to validate the weakly nonlinear predictions. The vortex sheet method is a standard sharp-interface technique for simulating dynamics of Hele-Shaw flows [22]. It is based on a boundary integral formulation in which the interface is formally replaced by a generalized vortex sheet [23] with a distribution of vortex strengths $\gamma(s,t)$, where $s$ is the arclength coordinate. We adapt this approach to handle ferrofluids under imposed magnetic fields. First, we express the velocity of the interface solely in terms of the interface position. To do so, it is convenient to identify the position vector in $\mathbb{R}^{2}$ with a scalar $z(s,t)\in\mathbb{C}$ (^∗ denotes complex conjugate) [23, 24, 25]. Second, to advance the interface, we solve the dimensionless equations

iteratively for the velocity $z_{t}$, where $(\cdot)_{t}\equiv\partial(\cdot)/\partial t$, $(\cdot)_{s}\equiv\partial(\cdot)/\partial s$, $i=\sqrt{-1}$, and $\mathcal{P}$ represents principal value integration. Here, $\Psi=\mathrm{N_{Ba}}r^{-2}+\mathrm{N_{Br}}r^{2}$, and $(\bm{\mathrm{M}}\cdot\bm{\mathrm{\hat{n}}})^{2}=\chi\left[\sqrt{\mathrm{N_{Ba}}}r^{-1}(\bm{\mathrm{\hat{e}}}_{\theta}\cdot\bm{\mathrm{\hat{n}}})+\sqrt{\mathrm{N_{Br}}}r(\bm{\mathrm{\hat{e}}}_{r}\cdot\bm{\mathrm{\hat{n}}})\right]^{2}$ is the dimensionless magnetic normal stress. Time advancement is accomplished by the Crank–Nicolson scheme. The spatial discretization is implemented on an array of Lagrangian points ($N=1024$) with uniform $\Delta s$; see Appendix B for further details, including algorithm flowchart and grid convergence study.

The evolution of perturbed harmonic modes $\xi_{k}$ under the fully nonlinear simulation and the weakly nonlinear approximation are shown in Fig. 2(a). Starting from small initial values ($\xi_{k}|_{k=7}=0.002$, $\xi_{nk}=0$ for $n>1$) with $\mathrm{N_{Ba}},\mathrm{N_{Br}},R,\chi$ set so that the most unstable mode is equal to the fundamental mode ($k_{m}=k=7$), they saturate at late times. The perturbed circular interface grows exponentially in the linear regime and then matches the weakly nonlinear approximation at intermediate times ($t\in[0,t_{w}]$). The nonlinear simulations take longer to saturate ($t\in[t_{w},t_{e}]$) and do so at higher final amplitudes compared to the weakly nonlinear result. The time-domain evolution is also shown in Fig. 2(b), evolving from a nearly flat (unwound circular) interface into a permanent propagating profile ¹¹1See Supplemental Material at [URL will be inserted by publisher] for a video of this process..

The rotating droplet, shown in Fig. 2(c), has a polygonal shape with the symmetry set by the fundamental mode, $k=7$. The fully nonlinear profile has a sharper peak compared to the weakly nonlinear approximation, which is otherwise in good agreement. The key discovery of the present work is the stable rotating shape, which we now seek to analyze as a nonlinear wave phenomenon [14].

A deficiency of linear and weakly nonlinear analyses is that they do not provide sufficient conditions for stability. Linearly stable base states can be nonlinearly unstable [27], and vice versa. Importantly, however, our nonlinear traveling wave solution is a local attractor (following the terminology from [28]); see Fig. 2(d).

Shapes in a neighborhood of the propagating profile, subject to small ($\xi_{k,2k}/\xi_{k,2k}^{f}\ll 1$) or intermediate ($\xi_{k,2k}/\xi_{k,2k}^{f}=\mathcal{O}(1)$) initial perturbations, converge to it. Larger perturbations (shaded region) lead to nonlinear instability of the weakly nonlinearly stable profiles; “fingers” continue to rotate and grow without bound under the effect of the radial magnetic field $\propto\mathrm{N_{Br}}$, which increases with distance to the center of the droplet. Convergence to the attractor is sensitive to the initial amplitude of the first harmonic mode $\xi_{k}$. For the chosen parameters, $\lambda(k)>0$ and $\lambda(2k)<0$: high wavenumbers decay and the fundamental wavenumber grow in the linear stage. Consequently, for low $\xi_{k}/\xi_{k}^{f}$ and high $\xi_{2k}/\xi_{2k}^{f}$, the low wavenumber modes grow and saturate, as high wavenumber modes decay exponentially in the linear regime. With higher initial $\xi_{k}/\xi_{k}^{f}$, the perturbed droplet will not go through the linear regime, and the amplitudes of both modes will rapidly increase to create a skewed shape, with multivalued $h(\theta,t)$, for which harmonic modes can no longer be defined. Note that Fig. 2(d) is a projection in the $(\xi_{k},\xi_{2k})$ plane, where the initial values of $\xi_{3k}$, $\xi_{4k}$, $\xi_{5k}$ are set as the final amplitudes (and phases) from the weakly nonlinear equations. A fast Fourier transform was used to decompose the nonlinear profile into normal modes that we plot in this figure. Note that even though Fig. 2(a) indicates $\xi_{3k}$ makes a non-trivial contribution to the final shape, while $\xi_{4k},\xi_{5k}$ play a smaller role, the projection is sufficient to conclude that the propagating wave profile is an attractor.

A permanent traveling wave profile has $\xi(\theta,t)=\Xi(k\theta-\omega t)$, and $v_{f}=\omega/k$ is its propagation velocity. Expressing the modes’ complex amplitudes as $\xi_{nk}(t)=c_{n}e^{-in\omega(k)t}$, with constant $c_{n}\in\mathbb{C}$ that account for their relative phases, we have $v_{p}(k,t)=n\omega(k)/nk=v_{f}$. The mean $v_{p}$ of the first five harmonics is used to calculate $v_{f}^{F}$ for the fully nonlinear simulation and also $v_{f}^{W}$ for the weakly nonlinear approximation. Meanwhile, $v_{f}^{L}=v_{p}$ as given by Eq. (12).

For a quantitative comparison, three sets of parameters are considered, fixing $\chi=1$. Two sets (i) and (ii) are for $k_{m}=7$, and the variation of $\mathrm{N_{Br}}$ is according to Eq. (11). A third set (iii) explores the effect of $k_{m}$ under the same linear propagation velocity $v_{f}^{L}$. Figure 3(a) compares the final propagating velocity predictions. Both $v_{f}^{L}$ and $v_{f}^{W}$ are in relatively good agreement with $v_{f}^{F}$ for small velocities. When $\mathrm{N_{Ba}}\to 0$ (the magnetic field becomes radial), only a stationary (non-rotating) droplet ($v_{f}=0$) exists [12]. For higher $v_{f}$, the larger deviation in the predictions highlights the importance of nonlinearity. Nevertheless, the linear and weakly nonlinear results follow a similar trend. Importantly, $v_{f}^{L}$ and $v_{f}^{W}$ help identify the key control factors: the coupled magnetic field strength $\sqrt{\mathrm{N_{Ba}}\mathrm{N_{Br}}}$ and the radius of the initial droplet $R$. The salient physics uncovered is that the propagating velocity can be non-invasively tuned.

The most unstable mode $k_{m}$ sets the propagating profile, which has a sharper peak for higher $k_{m}$ [Fig. 3(c)]. To quantify the shape change, we introduce the skewness $Sk(t)=\langle\xi^{3}\rangle/\langle\xi^{2}\rangle^{3/2}$, which is used to define the vertical asymmetry of nonlinear surface water waves [29, 30]; $Sk>0$ corresponds to narrow crests and flat troughs. Here, $\left\langle\,\cdot\,\right\rangle=\frac{1}{2}\pi\int_{0}^{2\pi}(\,\cdot\,)\,d\theta$. Figure 3(b) shows that $Sk$ (for the fully nonlinear propagation) increases with $k_{m}$, as expected from the sharper peaks in Fig. 3(c). This observation also explains why $v_{f}^{L}$ becomes a worse approximation of $v_{f}^{F}$ as $k_{m}$ increases [inset of Fig. 3(a)]: smoother peaks (lower $k_{m}$) are better captured by the linear theory based on harmonic modes.

Figure 3(c) reveals that the wave profile for $k_{m}=5$ ($\mathrm{N_{Ba}}=1.9$, $\mathrm{N_{Br}}=19.5$) is more asymmetric than the one for $k_{m}=9$ ($\mathrm{N_{Ba}}=0.60$, $\mathrm{N_{Br}}=60.8$). Under a purely radial magnetic field ($\mathrm{N_{Ba}}=0$), the stationary shape has azimuthal symmetry [12]. For the combined magnetic field, on the other hand, the dimensionless governing Eq. (2) and pressure boundary condition in Eq. (3) can be rewritten as

where $\bm{\mathrm{\hat{e}}}_{\theta}\cdot\bm{\mathrm{\hat{n}}}=-h_{\theta}/(h^{2}+h_{\theta}^{2})$, $\bm{\mathrm{\hat{e}}}_{r}\cdot\bm{\mathrm{\hat{n}}}=h/(h^{2}+h_{\theta}^{2})$, and $h_{\theta}=\partial h/\partial\theta$. The magnetic scalar potential in Eq. (14) results from the body force, and the terms pre-multiplied by $\chi$ in Eq. (15) represent the magnetic normal stress. For a droplet with symmetric azimuthal perturbation, the body force alone cannot break the symmetry. Therefore, the asymmetry of shapes discussed is to be attributed to the magnetic normal stress.

This observation can be intuitively understood by considering one wavelength of a symmetric waveform. The first three terms on the right-hand side of Eq. (15) are equal on both sides of the peak, while the fourth term changes at the peak due to the sign of $h_{\theta}$, which requires different curvatures on either side of the peak to remain balanced. Therefore, $\sqrt{\mathrm{N_{Ba}}\mathrm{N_{Br}}}$ can be taken as the measure of the coupling effect between the magnetic field components.

To further understand the asymmetry of propagating shapes induced by the combined magnetic field, we extend the parameters of case (i) to a new case (iv): $k_{m}=7$, $R=1$ and $\mathrm{N_{Br}}$ varying according to $\mathrm{N_{Ba}}$ (see Fig. 4 caption). To quantify the fore-aft asymmetry of the shape, we introduce $As(t)=\langle\mathcal{H}[\xi]^{3}\rangle/\langle\xi^{2}\rangle^{3/2}$ [29, 30]; $\mathcal{H}[\,\cdot\,]$ is the Hilbert transform. For $As>0$, waves tilt “forward” (i.e., counter-clockwise).

Figure 4(a) shows $As(t)$ for different $\sqrt{\mathrm{N_{Ba}}\mathrm{N_{Br}}}$, which quantifies the coupled field effect, starting with small symmetric perturbations. For a stable case, $As(t)$ reaches a maximum value ($t\approx t_{1}$) during the initial unstable weakly nonlinear growth (dark shadow region), and asymptotes to a value close to zero ($t\geq t_{6}$). The differences in the final propagating profile (under the same $k_{m}$) shown in Fig. 4(b) are hard to capture, which is consistent with the observation in Fig. 3(b). For the unstable cases, “wave breaking” occurs, which is highlighted by a change of sign of $As$. Also, now, $As(t)$ no longer saturates at late $t$. Instead $As(t)$ crosses zero (at $t\gtrsim t_{3}$) and approaches a singularity. This unstable example is shown in the second row of Fig. 4(c). As its amplitude first grows, the wave tilts forward ($t=t_{2}$), but nonlinear effects restore its symmetry ($t=t_{3}$). Subsequently, the wave tilts backwards ($t=t_{4},t_{5}$) and its amplitude continues to grow ($t=t_{6}$). The calculation of $As$ then fails because $\mathcal{H}$ requires the perturbation $\xi(\theta,t)$ to be single-valued in $\theta$. The distorted wave has a wider base and evolves into long unstable fingers.

Note that $\mathrm{N_{Ba}}/\mathrm{N_{Br}}$ also increases with $\sqrt{\mathrm{N_{Ba}}\mathrm{N_{Br}}}$ for our choices of $\mathrm{N_{Ba}}$ and $\mathrm{N_{Br}}$. Equation (9) shows that the radial magnetic field is destabilizing, while surface tension ($k>2$ here) and the azimuthal field are stabilizing. However, the nonlinear simulations indicate that, for the same $k_{m}$, increasing $\mathrm{N_{Ba}}/\mathrm{N_{Br}}$ can induce instability because it engenders a larger $v_{f}$ (and $As$), leading to a global bifurcation with Fig. 2(d) as one stable slice. This result has an analogy to solitary waves in equations of the Kortweg–de Vries (KdV) type. Specifically, initial perturbations grow, deforming a shape until nonlinearity is balanced by dispersion, when a permanent wave emerges [31]. However, depending on the form of the nonlinearity, not all such permanent waves are stable attractors, and conditions must be placed on the wave speed [32].

This study demonstrates how a perturbed circular ferrofluid droplet can evolve into a nonlinearly stable rotating shape. The most unstable mode sets how perturbations evolve into a permanent profile (and its skewness and asymmetry). Weakly nonlinear theory, in hand with fully nonlinear simulations, revealed permanent rotating shapes (traveling waves) with predictable propagation velocity. We showed how the coupling of the magnetic field components modifies the asymmetry and the nonlinear instability.

Although the manipulation of the linear growth rate of interfacial perturbations in Hele-Shaw cells is well studied [33], including extensions based on the weakly nonlinear expansion from Eq. (7) [34], the control of the dynamic, fully nonlinear, patterns is not. Our approach harnesses the magnitude and the direction of coupled magnetic fields to generate ferrofluid droplets, with well-characterized shapes and rotational speeds, by purely external means.

Open questions remain: e.g., which fundamental modes evolve into propagating shapes? Work on the stationary problem [17, 35] gives a hint, however, for a propagating shape the Birkhoff integral equation [36] must be solved, making an extension of [17, 35] challenging. Interestingly, our simulations also reveal that patterns predicted as stable by weakly nonlinear analysis can be unstable. In Appendix A, we provide an example showing that perturbations with $k=4$ will not evolve into either a stationary or a propagating shape (although both are predicted to exist by weakly nonlinear analysis).

Additionally, does this system accommodate more than one propagating wave? If so, do such waves keep their shapes upon collision, as with soliton interactions [31, 37]? Previous studies derived KdV equations for unidirectional small-amplitude, long-wavelength disturbances on fluid–fluid interfaces in Hele-Shaw [38] and axisymmetric ferrofluid configurations [19, 39], demonstrating the celebrated “$\operatorname{sech}^{2}$” solitary wave. Instead, in our study without such restrictions, we discovered periodic traveling nonlinear waves, which are akin to the cnoidal solutions of periodic KdV, i.e., the fundamental nonlinear modes (“soliton basis states”) [40]. Additionally, we observed wave breaking [Fig. 4(c,bottom)].

Finally, it would be of interest to verify the proposed shape manipulation strategies by laboratory experiments. Previous theoretical studies [17, 41, 42, 43] suggest that many exact stationary droplet shapes are unstable, thus their relevance to experimental studies is limited. On the other hand, the nonlinear simulations in our study, showing stable rotation, pave the way for future experimental realizations.